一种基于GPU的子空间快速求解方法技术

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

【技术实现步骤摘要】
一种基于GPU的子空间快速求解方法
本专利技术属于电子信息
，具体的说是一种基于GPU的子空间快速求解方法。
技术介绍
波达方向估计是阵列信号处理领域的重要研究方向，它在雷达、声纳、通信、导航等领域有着极为广阔的应用前景。基于子空间分解的测向算法是一种高分辨测向算法，包括多重信号分类算法(MUSIC，multiplesignalclassification)以及基于旋转不变技术的参数估计算法(ESPRIT，estimatingsignalparameterviarotationalinvariancetechniques)，能够实现对处于同一波束宽度内的多个信源的分离，但是这类算法需要对阵列信号协方差矩阵进行特征分解来获取子空间。对于单个信号源的情况，现有的幂法能够快速有效求解子空间。但是当信号源数量较多时，幂法不能直接给出用于构建子空间需要的全部特征向量，现有的方法一般采用Jacobi算法或者QR分解算法对阵列信号协方差矩阵进行特征分解来获取子空间。为了提高算法实时性，现有技术一般采用FPGA(FieldProgrammableGateArray)或者DSP(DigitalSignalProcessor)平台来加速子空间求解过程。文献1，徐德琛，刘志文，徐友根.某测向系统中MUSIC算法的FPGA实现[J].北京理工大学学报，2010.文献1提出在FPGA平台上采用round-robin-Jacobi算法求解子空间，在阵元数为12的情况下，子空间求解过程耗时约为0.2ms。然而，随着通信技术的发展，大规模天线阵列应用逐渐广泛，5G基站最大可配置的天线数量为1024根。在阵元数较多时，采用Jacobi算法或者QR分解算法求解子空间存在计算复杂度高、迭代时间长等问题，通过FPGA或者DSP平台加速子空间求解过程其加速效果仍不够，而且开发难度大，平台维护成本较高。近十年来，随着图形处理器(GPU，GraphicProcessUnit)的不断发展，其架构也在不断完善，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总线(peripheralcomponentinterconnectexpress)写入设备内存，采用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的具体方本文档来自技高网...

【技术保护点】
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计算阵列信号协方差矩阵

【技术特征摘要】
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即为...

【专利技术属性】
技术研发人员：常煜，万群，徐保根，何林立，王子强，
申请(专利权)人：电子科技大学，同方电子科技有限公司，
类型：发明
国别省市：四川,51