CN103226193B - 基于gpu架构的自适应递归最小二乘杂波对消方法 - Google Patents

Abstract

本发明公开了一种基于GPU架构的自适应递归最小二乘杂波对消方法，主要解决基于CPU架构的自适应杂波对消处理速度受限于硬件平台和对消处理难于实现并行的问题。其实现步骤为：在CPU主机端对杂波对消的参数进行初始化，并把用于杂波对消的数据传输至图形处理器GPU显存中；将递归最小二乘RLS杂波对消处理中涉及到的矩阵相乘、矩形相减、规约求和等运算，调度到GPU硬件平台上并行执行，得到杂波对消后的数据；将对消处理后的数据传输至主机端内存中，由CPU负责后续的数据处理工作。本发明具有硬件平台搭建简单、稳定性好、加速性能好、扩展性强、易于开发的优点，可用于缩短实际工程中杂波对消处理的时间。

Description

基于GPU架构的自适应递归最小二乘杂波对消方法

技术领域

本发明属于雷达技术领域，主要涉及杂波对消的实现方法，具体来说就是一种基于GPU架构的自适应递归最小二乘杂波对消方法，可在数学仿真软件MATLAB和微软提供的软件集成开发平台Visual studio上，快速进行杂波对消。

技术背景

杂波对消广泛应用于各种信号处理系统中，特别是外辐射雷达信号处理系统。杂波抑制一般位于信号处理的前端，是整个信号处理中最重要的模块，杂波抑制的好坏直接影响着后续信号处理的性能高低。杂波对消处理的目的主要是解决强直达波、多径杂波干扰的问题，通过主、辅两个通道的信号相关来对消主通道中掺杂的杂波，理想的情况下，主通道中的杂波可以被完全抑制掉。

杂波对消主要有最小均方误差LMS、直接矩阵求逆SMI和递归最小二乘RLS三种算法，其中，递归最小二乘RLS的杂波对消性能最好，精度最高，但运算量最大，处理时间也最长，这是RLS主要的缺点，如何提高RLS的处理速度是RLS算法存在的一个关键问题。现有技术中，RLS杂波对消是基于CPU架构实现的，杂波抑制速度受到了硬件的限制，不能满足各种信号处理系统的需求。

发明内容

本发明的目的在于针对上述基于CPU架构的RLS杂波对消的不足，提出一种基于GPU架构的自适应递归最小二乘杂波对消方法，以在不损失处理精度的情况下，加速杂波对消的处理速度，满足各种信号处理系统的需求。

本发明的技术思路是：使用图形处理器GPU代替CPU，利用GPU强大的浮点运算能力，对杂波对消中涉及到的矩阵运算进行并行加速，缩短杂波对消处理的时间。

为实现上述目的，本发明方法包括如下步骤：

(1)在CPU主机端对递归最小二乘RLS滤波器的参数进行初始化：

(1a)将雷达目标信号、杂波信号和噪声信号作为自适应递归最小二乘RLS滤波器主通道的接收信号sig；将杂波信号和噪声信号作为自适应递归最小二乘RLS滤波器辅助通道的接收信号ref，这两种信号sig、ref的数据长度S根据接收信号的采样率设定，并根据杂波强度设定杂波对消阶数L；

(1b)在CPU主机端分配杂波对消所需的内存、显存资源；

(1c)设置递归最小二乘RLS的遗忘因子λ，0＜λ＜1；

(1d)设递归最小二乘RLS杂波对消起始时刻k为1；

(1e)设杂波对消权向量w的初始值为：w_L×1＝[00...00]_L×1；

(1f)设递归最小二乘RLS杂波对消起始时刻误差信号ε＝sig；

(1g)在CPU主机端构造杂波对消起始时刻的协方差矩阵：P_L×L＝E_L×L/δ，其中，E为L阶单位方阵，δ为常数，取值为0.05；

(1h)在CPU主机端，把主通道信号sig、辅助通道信号ref的数据加载到CPU内存中；

(2)在CPU主机端启动杂波对消工作线程；

(3)选择能支持计算统一设备架构CUDA的图形处理器GPU，并将该图形处理器GPU与杂波对消线程进行连接；

(4)把CPU主机端内存中的主通道信号sig、辅助通道信号ref、误差信号ε、协方差矩阵P_L×L、杂波对消初始权向量w_L×1传输至图形处理器GPU的设备端显存中；

(5)在图形处理器GPU的设备端进行递归最小二乘RLS杂波对消处理：

(5a)在图形处理器GPU设备端上，获取k时刻的主通道信号和辅助通道信号：

(5a1)在图形处理器GPU的设备端上，加载k时刻的主通道信号sig_k+L-1；

(5a2)在图形处理器GPU的设备端上，加载k时刻辅助通道信号ref_k+L-1，并向前滑动L个时间单元，记为参考输入信号x_k，即：

x_k＝ref[kk+1...k+L-1]_L×1；

(5b)在图形处理器GPU的设备端调用内核函数计算中间矩阵： 其中代表k时刻的协方差矩阵，H代表共轭转置；

(5c)在图形处理器GPU设备端调用内核函数，更新k时刻的误差信号ε_k+L-1和k+1时刻的杂波对消权向量w_L×1；

(5d)在图形处理器GPU设备端调用内核函数，更新k+1时刻的协方差矩阵P_L×L；

(5e)在CPU主机端更新时刻值：k＝k+1；

(6)重复迭代步骤(5a)～(5e)共S-L+1次，在图形处理器GPU设备端上，得到杂波对消后的误差信号ε；

(7)把杂波对消后的误差信号ε传输至CPU主机端的内存并保存；

(8)在CPU主机端上释放进行杂波对消所分配的内存、显存资源。

本发明与现有技术相比具有以下优点：

1)本发明由于采用了基于浮点运算能力大的GPU架构系统，与CPU相比，数据扩展性强。

2)本发明将递归最小二乘RLS杂波对消处理中涉及到的矩阵运算，调度到图形处理器GPU上处理，提高了杂波对消的并行度，从而提高了杂波对消处理速度。

3)本发明中的硬件平台易于搭建，与专业的信号处理板相比，降低了开发费用，缩短了开发周期。

附图说明

图1是本发明应用的自适应滤波器原理图；

图2是本发明的流程图；

图3是在对消阶数L为64时，用本发明对仿真数据处理与CPU处理的时间对比图；

图4是在对消阶数L为128时，用本发明对仿真数据处理与CPU处理的时间对比图；

图5是在对消阶数L为256时，用本发明对仿真数据处理与CPU处理的时间对比图；

具体实施方式

参照图1，本发明应用的自适应滤波器是基于双通道实现的，其中一个通道为主通道，接收的信号包括目标信号、直达波信号、多径杂波信号和噪声信号，另一个为辅助通道，接收的信号包括直达波信号和噪声信号。主通道和辅助通道中的杂波信号必须相关才可以进行杂波对消，在对消时，用主通道的信号减去不同延时参考信号的加权和，就可以把主通道中的杂波干扰滤除掉，得到相对纯净的目标信号。

参照图2，本发明基于GPU架构的自适应递归最小二乘杂波对消方法，其实现步骤如下：

步骤1：在CPU主机端对递归最小二乘RLS滤波器的参数进行初始化。

(1a)将雷达目标信号、杂波信号和噪声信号作为自适应递归最小二乘RLS滤波器主通道的接收信号sig；将杂波信号和噪声信号作为自适应递归最小二乘RLS滤波器辅助通道的接收信号ref，这两种信号sig、ref的数据长度S根据接收信号的采样率设定，并根据杂波强度设定杂波对消阶数L；

(1b)在CPU主机端上调用Windows应用编程接口API中的malloc()函数和统一设备计算架构CUDA中的cudaMalloc()函数分配杂波对消所需的内存、显存资源；

(1c)设置递归最小二乘RLS的遗忘因子λ，0＜λ＜1；

(1d)设递归最小二乘RLS杂波对消起始时刻k为1；

(1e)设杂波对消权向量w的初始值为：w_L×1＝[00...00]_L×1；

(1f)设递归最小二乘RLS杂波对消起始时刻误差信号ε＝sig；

(1g)在CPU主机端构造杂波对消起始时刻的协方差矩阵：P_L×L＝E_L×L/δ，其中，E为L阶单位方阵，δ为常数，取值为0.05；

(1h)在CPU主机端，把主通道信号sig、辅助通道信号ref的数据加载到CPU内存中。

步骤2：在CPU主机端调用Windows应用编程接口API中的CreateThread()创建一个工作线程，杂波对消处理的所有操作均在此工作线程中进行。

步骤3：选择能支持计算统一设备架构CUDA的图形处理器GPU，并调用计算统一设备架构CUDA中的cudaSetDevice()函数将该图形处理器GPU与杂波对消线程进行连接，此后杂波对消工作线程就可以对图形处理器GPU进行调度。

步骤4：在CPU主机端调用计算统一设备架构CUDA中的同步数据拷贝函数cudaMemcpy()把CPU内存中的主通道信号sig、辅助通道信号ref、误差信号ε、协方差矩阵P_L×L、初始杂波对消权向量w_L×1传输至图形处理器GPU的设备端显存中。

步骤5：在图形处理器GPU的设备端进行递归最小二乘RLS杂波对消处理。

(5a)在图形处理器GPU设备端上，获取k时刻的主通道信号和辅助通道信号：

(5a1)在图形处理器GPU的设备端上，加载k时刻的主通道信号sig_k+L-1；

(5a2)在图形处理器GPU的设备端上，加载k时刻辅助通道信号ref_k+L-1，并向前滑动L个时间单元，记为参考输入信号x_k，即：

x_k＝ref[kk+1...k+L-1]_L×1；

(5b)在图形处理器GPU的设备端调用内核函数计算中间矩阵： 其中代表k时刻的协方差矩阵，H代表共轭转置：

(5b1)设内核函数线程格grid为(L/D,1)，线程块block为(L,1)，其中，(L/D,1)表示内核函数的线程格grid由L/D个线程块block组成，(L,1)表示内核函数的每一个线程块block由L个图形处理器GPU线程组成，D代表一个线程块block所能计算出M_L×1、N_1×L的数据点数，受图形处理器GPU存储资源限制，D取值为4，L/D个线程块block计算出M_L×1、N_1×L的所有L点数据；

(5b2)根据每个线程块block的索引号，在内核函数的每个线程块block中，同时加载协方差矩阵P_L×L中各自对应的D行数据，分别与当前k时刻的参考输入信号x_k进行规约求和运算，得到D点输出数据，并写入到M_L×1的对应位置；

(5b3)根据每个线程块block的索引号，在内核函数的每个线程块block中，同时加载协方差矩阵P_L×L中各自对应的D列数据，分别与k时刻的参考输入信号x_k进行规约求和运算，得到D点输出数据，并写入到N_1×L的对应位置；

(5c)在图形处理器GPU设备端调用内核函数，更新k时刻的误差信号ε_k+L-1和k+1时刻的杂波对消权向量w_L×1：

(5c1)设内核函数线程格grid为(1,1)，线程块block为(L,1)，其中(1,1)表示内核函数由一个线程块block组成，(L,1)表示线程块block由L个图形处理器GPU线程组成；

(5c2)规约计算k时刻自适应滤波器的输出值：H代表共轭转置；

(5c3)更新k时刻的误差信号为：ε_k+L-1＝sig_k+L-1-T；

(5c4)规约计算中间变量

(5c5)并行计算增益向量K_L×1的L点数据：K_L×1＝M_L×1/(λ+F)；

(5c6)并行更新k+1时刻杂波对消权向量w_L×1的L点权系数为：

$w_{L\×1}^{(k+1)}=w_{L\×1}^{\left(k\right)}+K_{L\×1}*{\mathrm{\ϵ}}_{k+L-1}^H;$

(5d)在图形处理器GPU设备端调用内核函数，更新k+1时刻的协方差矩阵P_L×L：

(5d1)设内核函数线程格grid为((L+15)/16,(L+15)/16)，线程块block为(16,16)，其中((L+15)/16,(L+15)/16)表示内核函数的线程格grid由((L+15)/16)*((L+15)/16)个线程块block组成，(16,16)表示内核函数的一个线程块block由16*16个图形处理器GPU线程组成；

(5d2)计算中间矩阵：U_L×L＝K_L×1*N_1×L；

(5d3)并行更新k+1时刻协方差矩阵P_L×L的L*L点数据：

(5e)在CPU主机端更新时刻值：k＝k+1。

步骤6：重复迭代步骤(5a)～(5e)共S-L+1次，在图形处理器GPU设备端上，得到杂波对消后的误差信号ε。

步骤7：在CPU主机端调用统一计算设备架构CUDA中的同步数据拷贝函数cudaMemcpy()把杂波对消后的误差信号ε传输至CPU主机端的内存并保存。

步骤8：在CPU主机端上调用Windows应用编程接口API中的free()函数和统一设备计算架构CUDA中的cudaFree()函数，释放进行杂波对消所分配的内存、显存资源。

本发明的效果可通过以下仿真进一步说明：

1）实验条件：

在实验中，硬件平台选用HP Z820工作站，GPU显卡型号为NVIDA Telsa C2075，Intel Xeon多核处理器，Win7系统，软件平台为Visual Stdio2008+CUDA4.0和MATLAB2009b。

实验采用一帧外辐射源雷达仿真数据，其中主通道中包括两个目标信号、直达波信号、多径信号和噪声信号，辅助通道中包括直达波信号和噪声信号，另外，主通道和辅助通道的数据量均为200000点复数数据。

2）实验内容及结果：

采用本发明提出的方法，在对消阶数L为64、128、256的情况下，采用单精度浮点型float，对一帧外辐射源雷达的仿真数据进行杂波对消并保存，滤除主通道中含有的直达波信号和多径信号。实验中，λ＝0.95、δ＝0.05。在同样的条件下，分别进行15次重复实验，记录每次实验的运行时间，并与CPU处理的时间进行对比，实验结果如图3、图4、图5。其中，图3为对消阶数L为64时GPU处理与CPU处理时间对比图，图4为对消阶数L为128时GPU处理与CPU处理时间对比图，图5为对消阶数L为256时GPU处理与CPU处理时间对比图。

由图3、图4、图5可见，使用本发明提出的方法进行杂波对消的时间均显著小于MATLAB处理的时间，当L=64时，平均加速比为19倍，当L=128时，平均加速比为38倍，当L=256时，平均加速比为85倍，而且处理时间很稳定，验证了本发明提出的方法具有良好的稳定性和加速性能。

Claims (4)

1.一种基于GPU架构的自适应递归最小二乘杂波对消方法，包括如下步骤：

(1)在CPU主机端对递归最小二乘RLS滤波器的参数进行初始化：

(1a)将雷达目标信号、杂波信号和噪声信号作为自适应递归最小二乘RLS滤波器主通道的接收信号sig；将杂波信号和噪声信号作为自适应递归最小二乘RLS滤波器辅助通道的接收信号ref，这两种信号sig、ref的数据长度S根据接收信号的采样率设定，并根据杂波强度设定杂波对消阶数L；

(1b)在CPU主机端分配杂波对消所需的内存、显存资源；

(1c)设置递归最小二乘RLS的遗忘因子λ，0＜λ＜1；

(1d)设递归最小二乘RLS杂波对消起始时刻k为1；

(1e)设杂波对消权向量w的初始值为：w_L×1＝[00…00]_L×1；

(1f)设递归最小二乘RLS杂波对消起始时刻误差信号ε＝sig；

(1g)在CPU主机端构造杂波对消起始时刻的协方差矩阵：P_L×L＝E_L×L/δ，其中，E为L阶单位方阵，δ为常数，取值为0.05；

(1h)在CPU主机端，把主通道信号sig、辅助通道信号ref的数据加载到CPU内存中；

(2)在CPU主机端启动杂波对消工作线程；

(3)选择能支持计算统一设备架构CUDA的图形处理器GPU，并将该图形处理器GPU与杂波对消线程进行连接；

(4)把CPU主机端内存中的主通道信号sig、辅助通道信号ref、误差信号ε、协方差矩阵P_L×L、杂波对消初始权向量w_L×1传输至图形处理器GPU的设备端显存中；

(5)在图形处理器GPU的设备端进行递归最小二乘RLS杂波对消处理：

(5a)在图形处理器GPU设备端上，获取k时刻的主通道信号和辅助通道信号：

(5a1)在图形处理器GPU的设备端上，加载k时刻的主通道信号sig_k+L-1；

(5a2)在图形处理器GPU的设备端上，加载k时刻辅助通道信号ref_k+L-1，并向前滑动L个时间单元，记为参考输入信号x_k，即：

x_k＝ref[kk+1…k+L-1]_L×1；

(5b)在图形处理器GPU的设备端调用内核函数计算中间矩阵： 其中代表k时刻的协方差矩阵，H代表共轭转置；

(5c)在图形处理器GPU设备端调用内核函数，更新k时刻的误差信号ε_k+L-1和k+1时刻的杂波对消权向量w_L×1；

(5d)在图形处理器GPU设备端调用内核函数，更新k+1时刻的协方差矩阵P’_L×L；

(5e)在CPU主机端更新时刻值：k＝k+1；

(6)重复迭代步骤(5a)～(5e)共S-L+1次，在图形处理器GPU设备端上，得到杂波对消后的误差信号ε’；

(7)把杂波对消后的误差信号ε’传输至CPU主机端的内存并保存；

(8)在CPU主机端上释放进行杂波对消所分配的内存、显存资源。

2.根据权利要求1所述的基于GPU架构的自适应递归最小二乘杂波对消方法，其特征在于步骤(5b)所述的在图形处理器GPU的设备端调用内核函数计算中间矩阵： $M_{L\×1}=P_{L\×L}^{\left(k\right)}*x_k,N_{1\×L}=x_k^H*P_{L\×L}^{\left(k\right)},$ 按如下步骤进行：

(5b1)设内核函数线程格grid为(L/D,1)，线程块block为(L,1)，其中，(L/D,1)表示内核函数的线程格grid由L/D个线程块block组成，(L,1)表示内核函数的每一个线程块block由L个图形处理器GPU线程组成，D代表一个线程块block所能计算出M_L×1、N_1×L的数据点数，取值为4，L/D个线程块block计算出M_L×1、N_1×L的所有L点数据；

(5b2)在内核函数的每个线程块block中同时加载协方差矩阵P_L×L的D行数据，分别与当前k时刻的参考输入信号x_k进行规约求和运算，得到D点输出数据，并写入到M_L×1的对应位置；

(5b3)在内核函数的每个线程块block中同时加载协方差矩阵P_L×L的D列数据，分别与k时刻的参考输入信号x_k进行规约求和运算，得到D点输出数据，并写入到N_1×L的对应位置。

3.根据权利要求1所述的基于GPU架构的自适应递归最小二乘杂波对消方法，其特征在于步骤(5c)所述的在图形处理器GPU设备端调用内核函数，更新k时刻的误差信号ε_k+L-1和k+1时刻的杂波对消权向量w_L×1，按如下步骤进行：

(5c1)设内核函数线程格grid为(1,1)，线程块block为(L,1)，其中(1,1)表示内核函数由一个线程块block组成，(L,1)表示线程块block由L个图形处理器GPU线程组成；

(5c2)规约计算k时刻自适应滤波器的输出值：H代表共轭转置；

(5c3)更新k时刻的误差信号为：ε_k+L-1＝sig_k+L-1-T；

(5c4)规约计算中间变量

(5c5)并行计算增益向量K_L×1的L点数据：K_L×1＝M_L×1/(λ+F)，F为中间变量；

(5c6)并行更新k+1时刻杂波对消权向量w_L×1的L点权系数为：

$w_{L\×1}^{(k+1)}=w_{L\×1}^{\left(k\right)}+K_{L\×1}*{\mathrm{\ϵ}}_{k+L-1}^H.$

4.根据权利要求1所述的基于GPU架构的自适应递归最小二乘杂波对消方法，其特征在于步骤(5d)所述的在图形处理器GPU设备端调用内核函数，更新k+1时刻的协方差矩阵P_L×L，按如下步骤进行：

(5d1)设内核函数线程格grid为((L+15)/16,(L+15)/16)，线程块block为(16,16)，其中((L+15)/16,(L+15)/16)表示内核函数的线程格grid由((L+15)/16)*((L+15)/16)个线程块block组成，(16,16)表示内核函数的一个线程块block由16*16个图形处理器GPU线程组成；

(5d2)计算中间矩阵U_L×L＝K_L×1*N_1×L，K_L×1为并行计算增益向量；

(5d3)并行更新k+1时刻协方差矩阵P_L×L的L*L点数据：

$P_{L\×L}^{(k+1)}=(P_{L\×L}^{\left(k\right)}-U_{L\×L})/\mathrm{\λ}.$