专利摘要

本发明属于电子信息技术领域，具体的说是一种基于GPU的子空间快速求解方法。本发明在幂法的基础上，利用阵列信号协方差矩阵的共轭对称性，提出使用Householder变换对阵列信号协方差矩阵进行降阶，得到一种快速子空间求解方法。同时利用GPU多核并行处理能力，进一步提高子空间求解过程实时性。

权利要求

1.一种基于GPU的子空间快速求解方法，其特征在于，包括以下步骤：

S1、初始化：

设置阵元数量为M，空间来波信号数量为K，误差范围eps＝1；

设定接收机分别在每个采样时刻对所接收的信号进行采样，并将数据写入主机内存，第n个采样时刻的采样序列表示为x(n)＝[x1(n),x2(n),...,xM(n)]T，T表示转置；

S2、将N个时刻的采样数据x(n),n＝1,2,..,N通过PCIe总线写入设备内存，采用GPU计算阵列信号协方差矩阵其中，*表示共轭转置，阵列信号协方差矩阵R是M×M的共轭对称矩阵；

S3、在GPU中，利用幂法直接计算阵列信号协方差矩阵R的主特征值λ1及主特征向量u1；

S4、判断空间来波信号数量K＝1是否成立，若是，则进入步骤S11，若否，则进入步骤S5；

S5、利用向量u1构造Householder变换矩阵H1，对阵列信号协方差矩阵R执行Householder变换，变换后的矩阵为

S6、选取矩阵B1的后M-1行和后M-1列构成矩阵R1，利用幂法计算矩阵R1的主特征值λ2及主特征向量v1，然后对向量v1进行升维，升维后的向量u2＝H1[0；v1]即是阵列信号协方差矩阵R第二大特征值对应的特征向量；

S7、判断空间来波信号数量K＝2是否成立，若是，则进入步骤S11，若否，则令k＝3，并进入步骤S8；

S8、利用向量vk-2构造Householder变换矩阵Hk-1，对矩阵Rk-2执行Householder变换，变换后的矩阵为

S9、选取矩阵Bk-1的后M-(k-1)行和M-(k-1)列构成矩阵Rk-1，利用幂法计算矩阵Rk-1的主特征向量vk-1，将向量vk-1升维为M维向量uk，即：

向量uk即为阵列信号协方差矩阵R第k大特征值对应的特征向量；

S10、若k等于K，进入步骤S11，否则，k＝k+1，回到步骤S8；

S11、输出阵列信号协方差矩阵R的K个大特征值对应的特征向量u1,u2,...,uK，即子空间求解结果。

2.根据权利要求1所述的一种基于GPU的子空间快速求解方法，其特征在于，步骤S3的具体方法为：

给定一个初始向量向量的第一个元素为1，其余元素为0；

重复执行直到时迭代结束，迭代结束后主特征向量其中eps为预先给定的误差范围；每轮迭代后对向量进行归一化，具体操作为

3.根据权利要求2所述的一种基于GPU的子空间快速求解方法，其特征在于，步骤S5的具体方法为：

给定向量e1，向量e1的第一个元素为u1的第一个元素，其余元素为0；

用向量e1构造Householder变换矩阵H1，其中：

对阵列信号协方差矩阵进行Householder变换，变换后的矩阵为矩阵B1的表现形式如下

4.根据权利要求3所述的一种基于GPU的子空间快速求解方法，其特征在于，步骤S6的具体方法为：

选取矩阵B1的后M-1行和后M-1列构成新的M-1维矩阵R1，矩阵R1包含矩阵R的特征值λ2,λ3,...,λM，将步骤3中矩阵R替换为R1，利用幂法计算矩阵R1主特征向量v1，对向量v1进行升维，升维后的向量即为阵列信号协方差矩阵R第二大特征值对应的特征向量。

5.根据权利要求4所述的一种基于GPU的子空间快速求解方法，其特征在于，步骤S8的具体方法为：

利用向量vk-2构造Householder变换矩阵Hk-1，其中：

向量ek-1第一个元素为vk-2第一个元素，其余元素为0。对矩阵Rk-2执行Householder变换，变换后的矩阵为

说明书

技术领域

本发明属于电子信息技术领域，具体的说是一种基于GPU的子空间快速求解方法。

背景技术

波达方向估计是阵列信号处理领域的重要研究方向，它在雷达、声纳、通信、导航等领域有着极为广阔的应用前景。基于子空间分解的测向算法是一种高分辨测向算法，包括多重信号分类算法(MUSIC，multiple signal classification)以及基于旋转不变技术的参数估计算法(ESPRIT，estimating signal parameter via rotational invariancetechniques)，能够实现对处于同一波束宽度内的多个信源的分离，但是这类算法需要对阵列信号协方差矩阵进行特征分解来获取子空间。对于单个信号源的情况，现有的幂法能够快速有效求解子空间。但是当信号源数量较多时，幂法不能直接给出用于构建子空间需要的全部特征向量，现有的方法一般采用Jacobi算法或者QR分解算法对阵列信号协方差矩阵进行特征分解来获取子空间。为了提高算法实时性，现有技术一般采用FPGA(FieldProgrammable Gate Array)或者DSP(Digital Signal Processor)平台来加速子空间求解过程。文献1，徐德琛，刘志文，徐友根.某测向系统中MUSIC算法的FPGA实现[J].北京理工大学学报，2010.文献1提出在FPGA平台上采用round-robin-Jacobi算法求解子空间，在阵元数为12的情况下，子空间求解过程耗时约为0.2ms。然而，随着通信技术的发展，大规模天线阵列应用逐渐广泛，5G基站最大可配置的天线数量为1024根。在阵元数较多时，采用Jacobi算法或者QR分解算法求解子空间存在计算复杂度高、迭代时间长等问题，通过FPGA或者DSP平台加速子空间求解过程其加速效果仍不够，而且开发难度大，平台维护成本较高。

近十年来，随着图形处理器(GPU，Graphic Process Unit)的不断发展，其架构也在不断完善，GPU的功能不再仅限于图形渲染，逐渐发展成为具备超高浮点运算能力的计算工具，为解决子空间求解实时性问题提供了有力工具。GPU的并行处理能力是由其片内大量的流处理器决定的，这种架构使得GPU在处理数据量大，数据之间依赖性较小的数据时相对于传统的CPU串行运算具备更大的运算优势。

发明内容

本发明提供一种基于GPU的信号子空间快速求解方法，能够有效解决信号子空间求解实时性差的问题。本发明在幂法的基础上，利用阵列信号协方差矩阵的共轭对称性，提出使用Householder变换对阵列信号协方差矩阵进行降阶，得到一种快速子空间求解方法。同时利用GPU多核并行处理能力，进一步提高子空间求解过程实时性。

为实现上述技术目的，本发明采用如下技术方案进行实现。

一种基于GPU的子空间快速求解方法，其特征在于，包括以下步骤：

S1、初始化：

设置阵元数量为M，空间来波信号数量为K，误差范围eps＝1；

设定接收机分别在每个采样时刻对所接收的信号进行采样，并将数据写入主机内存，第n个采样时刻的采样序列表示为x(n)＝[x1(n),x2(n),...,xM(n)]T，T表示转置；

S2、将N个时刻的采样数据x(n),n＝1,2,..,N通过PCIe总线(peripheralcomponent interconnect express)写入设备内存，采用GPU计算阵列信号协方差矩阵 其中，*表示共轭转置，阵列信号协方差矩阵R是M×M的共轭对称矩阵；

S3、在GPU中，利用幂法直接计算阵列信号协方差矩阵R的主特征值λ1及主特征向量u1；

S4、判断空间来波信号数量K＝1是否成立，若是，则进入步骤S11，若否，则进入步骤S5；

S5、利用向量u1构造Householder变换矩阵H1，对阵列信号协方差矩阵R执行Householder变换，变换后的矩阵为

S6、选取矩阵B1的后M-1行和后M-1列构成矩阵R1，利用幂法计算矩阵R1的主特征值λ2及主特征向量v1，然后对向量v1进行升维，升维后的向量u2＝H1[0；v1]即是阵列信号协方差矩阵R第二大特征值对应的特征向量；

S7、判断空间来波信号数量K＝2是否成立，若是，则进入步骤S11，若否，则令k＝3，并进入步骤S8；

S8、利用向量vk-2构造Householder变换矩阵Hk-1，对矩阵Rk-2执行Householder变换，变换后的矩阵为

S9、选取矩阵Bk-1的后M-(k-1)行和M-(k-1)列构成矩阵Rk-1，利用幂法计算矩阵Rk-1的主特征向量vk-1，将向量vk-1升维为M维向量uk，即：

向量uk即为阵列信号协方差矩阵R第k大特征值对应的特征向量；

S10、若k等于K，进入步骤S11，否则，k＝k+1，回到步骤S8；

S11、输出阵列信号协方差矩阵R的K个大特征值对应的特征向量u1,u2,...,uK，即子空间求解结果。

上述方案的步骤S2中，对于阵列信号协方差矩阵R的计算，即 每个元素的计算过程都是相对独立的，是一种典型的适合GPU并行实现的运算。

在GPU上计算阵列信号协方差矩阵，一般的实现方式是分配和协方差矩阵元素数量相等的线程，每个线程负责读取全局内存中矩阵x的一行和x*的一列进行乘累加。然而这种方法的运算效率受限于全局内存带宽，频繁访问全局内存会对整个运算过程造成较大的延迟，全局内存的带宽一般为100GB/s。

为了避免频繁访问全局内存，我们需要在存储器层面进行优化。因此本发明提出使用矩阵分块的方法，将矩阵相乘的过程分割成很多小矩阵的乘法。矩阵分块相乘运算如附图-3所示，例如计算R中R(1:16,1:16)的值等价于 其中R(1:16,1:16)表示矩阵R第1到15行，1到15列的元素。此时，可以将小矩阵加载到共享内存中(sharedmemory)，小矩阵的乘法运算就不需要再频繁访问全局内存，共享内存的访问延迟可以忽略。设置小矩阵的维度为k，经过分块优化之后，实际上需要的全局内存访问次数变为2M3/k，理论上数据存取次数变为原来的1/k。

当采样点数L较大时，(通常L>512，不同的GPU平台略有不同)，可采用并行归约法进行求和运算，首先求出相邻两个元素之和，然后在上一次计算结果的基础上对相邻两个和再进行求和运算，以此类推。归约运算如附图2所示。

可以看出，理论上这种基于GPU的归约求和只需要logL的时间即可完成CPU上L时间的运算，很大程度上降低了计算时间，同时，L越大，这种加速效果越好。

进一步的，步骤S3的具体方法为：

给定一个初始向量 向量 的第一个元素为1，其余元素为0；

重复执行 直到 时迭代结束，迭代结束后主特征向量 其中eps为预先给定的误差范围；每轮迭代后对向量 进行归一化，具体操作为

在步骤S3中，为满足合并访问，最大化利用GPU带宽，本发明提出将步骤2中得到的阵列信号协方差矩阵按列优先进行存储，然后在GPU上执行 通过多次调用kernel来达到迭代效果。

进一步的，步骤S5的具体方法为：

给定向量e1，向量e1的第一个元素为u1的第一个元素，其余元素为0；

用向量e1构造Householder变换矩阵H1， 其中：

对阵列信号协方差矩阵进行Householder变换，变换后的矩阵为 矩阵B1的表现形式如下

进一步的，步骤S6的具体方法为：

选取矩阵B1的后M-1行和后M-1列构成新的M-1维矩阵R1，矩阵R1包含矩阵R的特征值λ2,λ3,...,λM，将步骤3中矩阵R替换为R1，利用幂法计算矩阵R1主特征向量v1，对向量v1进行升维，升维后的向量 即为阵列信号协方差矩阵R第二大特征值对应的特征向量。

进一步的，步骤S8的具体方法为：

利用向量vk-2构造Householder变换矩阵Hk-1， 其中：

向量ek-1第一个元素为vk-2第一个元素，其余元素为0。对矩阵Rk-2执行Householder变换，变换后的矩阵为

本发明的有益效果为，本发明在幂法的基础上，针对阵列信号协方差矩阵的共轭对称性，提出一种快速子空间求解方法，同时利用GPU的多核并行处理能力，能够有效提高子空间求解实时性。

附图说明

图1为本发明一种基于GPU的信号子空间快速求解方法流程图。

图2为并行归约线程调度图。

图3为矩阵分块乘法示意图

图4为优化后GPU上存储器带宽示意图。

图5为GPU上资源利用率简图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述，以便本领域的技术人员更好地理解本发明。

实施例

参照图1，为本发明一种基于GPU的信号子空间快速求解方法实现流程。本例采用的是基于CPU+GPU异构计算平台。

本实施例以来波信号数量K＝3，阵元数M＝64,128,256,512为例来描述本发明在计算信号子空间时的效果。

本例的具体实施方式流程如下：

步骤1：步骤1：阵列天线接收到的数据经采样后首先存储在主机内存，首先通过cudaMemcpy()操作将采样数据从主机内存拷贝到设备内存，在GPU端计算阵列信号协方差矩阵R。

步骤2：为了避免频繁访问全局内存，我们采用矩阵分块乘法以及共享内存对内存访问进行优化。矩阵分块乘法如图-3所示。在CPU端申请二维block，block的大小为(16,16)，grid大小为((M+15)/16,(M+15)/16)。每个block负责计算阵列信号协方差矩阵中每个块矩阵，每个块矩阵大小为16*16。每个block通过不断将x的行块和x*的列块写入共享存储器，然后通过乘累加得到阵列信号协方差矩阵中某个块元素。图-5展示了这种设计模式下的GPU利用率，从图中可以看出GPU理论利用率达到100％，实际利用率达到99.19％。此外，图-4展示了存储器的访问带宽，共享存储器的带宽达到669.72GB/s。

步骤3：给定一个初始向量 在GPU端执行幂迭代过程 分配的线程数量与矩阵R维度一致。为了满足合并访问，在该步骤中，矩阵R按照列优先进行存储，每个线程读取矩阵R中的一列和向量 进行乘累加，得到向量 对应位置元素。在该步操作中，通过列优先存储进行优化，全局内存访问带宽能够达到57GB/s(极限带宽为80.2GB/s)。

向量乘累加是个串行执行的过程，为了提高并行度，采用归约运算计算两个向量的內积。

图2给出了并行归约线程调度图。假设我们需要对一个拥有8个元素的数组进行求和，第一次归约时S＝1，每个线程对数组中相邻两个位置元素进行求和，并将求和结果统一存储在左侧存储区。

第二次归约时S＝2，使用的线程数量为上一次执行归约运算的线程数量的一半，此时，线程索引为tid＝0,2,4,6的线程对相邻两个位置进行求和运算。

第三次归约S＝3，线程索引为tid＝0,4的线程访问相应位置元素进行求和运算。对于长度为L的数组，只需要logL次归约即可完成L次的乘加运算。

在步骤3中，通过多次调用kernel来完成迭代过程。迭代门限eps设置为1，迭代结束后向量 向量u1即为矩阵R最大特征值对应的特征向量。

步骤4：根据向量u1构造变换矩阵H1， 其中

向量e1第一个元素为u1第一个元素，其余元素为0。

为了避免设备内存和主机内存之间低效率的数据交换，步骤4也在GPU端进行计算。由于该步骤中运算量较小，对该步进行优化对整体加速效果意义很小，因此该步不做优化。

在GPU端对矩阵R执行Householder变换， 该步是基本的矩阵相乘运算，因此和步骤2中实施方式相同。

步骤5：选取矩阵B1的后M-1行和后M-1列构成矩阵R1，将步骤3中的矩阵R替换为R1，计算矩阵R1的主特征向量v1，然后对向量v1进行升维，升维后的向量u2＝H1[0；v1]即是阵列信号协方差矩阵R第二大特征值对应的特征向量。

步骤6：利用向量v1构造Householder变换矩阵H2， 其中

向量e2第一个元素等于v1的第一个元素，其余元素为0。对矩阵R1进行Householder变换，变换后的矩阵为 取矩阵B2的后M-2行和后M-2列构成新的矩阵R2，用R2替换步骤3中的矩阵R，得到矩阵R2的主特征向量v2。对向量v2进行升维，升维后的向量 即为矩阵R第三大特征值对应的特征向量。

将向量u1,u2,u3传输至主机内存，完成子空间求解。

为了验证本发明所提方法的有效性以及高效性，表1展示了子空间求解过程的加速效果，表-2给出了不同方法的子空间求解数值实验结果。采用统一硬件环境和开发平台进行对比，GPU型号为NVIDIA GeForce GTX 960M，CPU型号为Intel Core i5-4210。CUDA版本为CUDA7.5，编程环境为Visual Studio 2010，编程语言为CUDA C。

表1给出了在不同阵元数量的情况下本发明在计算信号子空间时的加速效果。

表1子空间求解加速效果。

本发明的基于GPU架构的信号子空间快速求解方法，通过提出一种快速、高并发的信号子空间求解方法，利用GPU强大的多核并行处理能力，能够有效解决信号子空间求解实时性差的问题。

表2在阵元数为64，空间来波信号数量为3的情况下，给出了使用本发明所提方法求解子空间的数值实验结果，同时给出了使用MATLAB eig()函数求解子空间的数值实验结果。为了进一步说明本发明所提方法在求解子空间时的有效性，我们给出了子空间向量之间的空间夹角余弦值，空间夹角余弦值的范围在[-1,1]之间，值越趋近于1，表示两个向量的方向越接近，越趋近于-1，表示两个向量的方向越相反，越接近于0，表示两个向量越近乎于正交。

cosθ1＝0.9999

cosθ2＝0.9993

cosθ3＝0.9977

其中θ1表示本发明所提方法最大特征值对应的子空间向量与MATLAB eig()函数得到的最大特征值对应的子空间向量之间的空间夹角。θ2表示本发明所提方法第二大特征值对应的子空间向量与MATLAB eig()函数得到的第二大特征值对应的子空间向量之间的空间夹角。θ3表示本发明所提方法第三大特征值对应的子空间向量与MATLAB eig()函数得到的第三大特征值对应的子空间向量之间的空间夹角。

表2本例所提方法的数值实验结果

可得出：本发明的基于GPU的子空间快速求解方法，在阵元数量为64～512时，能够将子空间求解的速度分别提高3.5～24.3倍。

一种基于GPU的子空间快速求解方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。