CN104133200A - 基于fpga的正交匹配追踪方法 - Google Patents

Abstract

本发明公开了一种基于FPGA的正交匹配追踪方法，主要解决现有雷达系统处理复杂算法能力有限及处理速度慢的问题，其实现方法是：根据雷达回波数据得到观测向量和字典矩阵；初始化误差余量为观测向量并对其补零至2048长度后调用FFT核对其做逆傅里叶变换；对FFT核输出的数据取模的平方并找到最大值的位置，即为误差余量相关性最大的列位置索引；根据找到的列位置索引更新列向量集合；对列向量集合进行施密特正交化处理更新误差余量；对误差余量取二范数并判断该二范数是否小于1，若满足此条件则用共轭梯度算法恢复出原始信号。本发明具有速度更快，资源占用更低的优点，可用于对基于傅里叶基的大规模字典矩阵和高稀疏度信号进行重建。

Description

基于FPGA的正交匹配追踪方法

技术领域

本发明属于数字信号处理技术领域，特别涉及利用观测信号进行完整信号重建算法的加速硬件实现方法，可用于合成孔径雷达/逆合成孔径雷达、步进频率雷达成像。

背景技术

雷达进行目标检测时，由于受到复杂电磁环境、雷达系统本身的缺陷、目标运动姿态变化等影响，可能会导致回波数据稀疏采样，若采用传统的傅里叶分析方法，会出现目标点散焦和图像模糊等现象，因此如何实现稀疏信号目标检测和成像非常重要。2006年由D.Donoho、E.Candes、T.Tao等人提出的一种全新的信息获取指导理论，即压缩感知CS，该理论指出：对稀疏或可压缩的信号通过远低于奈奎斯特标准的采样率进行数据采样，通过解一个非线性最优化问题，即可精确地恢复出原始信号。该理论一经提出，就在信号/图像处理、医疗成像、模式识别、光学/雷达成像、信息论、无线通信等领域受到高度关注。2007年，R.Baraniuk等人首次正式讨论了压缩感知原理在雷达成像中的应用，他们通过理论分析和数值仿真证明了压缩感知在雷达成像领域中的可行性。压缩感知理论的建立和快速发展，为稀疏信号重构技术的实现提供了理论依据。

稀疏信号重构需要通过求解l₁范数优化问题得到，l₁范数优化算法虽然计算精度较高，但是算法的计算复杂度极高，不利于实时处理与移动平台的实现。运算效率较高的求解方法为贪婪算法，主要包括匹配追踪算法MP、正交匹配追踪算法OMP、分布正交匹配追踪算法StOMP、正则化匹配追踪算法ROMP等。其中，OMP算法具有极高的运算效率，基于硬件平台的算法实现具有较高的应用价值，但对现有的OMP算法进行FPGA实现时仍存在以下缺点：

1)用FPGA实现OMP算法时，会大量涉及到矩阵—向量相乘、向量—向量相乘操作，而随着字典矩阵规模的增加，矩阵—向量相乘所占用的时钟周期和系统资源会急剧增加，从而对求解速度的提升和硬件资源的优化造成限制。

2)在求解最小二乘部分，现有的求解方法多是通过LU分解、QR分解等对矩阵进行求逆来解最小二乘问题，这种方法的缺点是当矩阵规模增大时，计算复杂度会大大增加，相应的计算时间会增加，进而降低了对信号的处理速度。

3)现有技术所使用的数据类型多为定点型的实数，而雷达接收到的数据则是含有相位信息的复数，算法的实用性不强。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于FPGA的正交匹配追踪方法，以对基于稀疏傅里叶基的OMP算法进行优化，在满足复杂运算的同时减少系统资源和计算延时，提高信号重构速度和实用性。

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

(1)根据雷达回波数据，利用Matlab软件产生一个2048×2048大小的傅里叶基矩阵Φ和512×2048大小的随机矩阵A，用该随即矩阵A分别对傅里叶基矩阵Φ和雷达回波数据进行抽取生成512×2048大小的字典矩阵Ψ和512长度的测量向量y，并将该字典矩阵Ψ和测量向量y中的数据量化成32bit单精度浮点型数据后存储到ROM中；

(2)初始化迭代次数t为1，初始化误差余量R₀为测量向量y，初始化列向量集合和位置索引集合Λ₀为空集；

(3)计算字典矩阵Ψ与误差余量R_t-1的相关性：即按随机抽取的先验知识对512长度的误差余量R_t-1补零生成2048长度的向量再对该向量进行IFFT运算；

(4)对步骤(3)的计算结果取模的平方得到2048长度的实向量h_t，求出h_t中的最大元素所处的位置λ_t，该位置λ_t即为字典矩阵Ψ中与误差余量R_t-1相关性最大的列所处的位置，再利用该位置λ_t更新列向量集合和位置索引集合Λ_t；

(5)对列向量集合进行施密特正交化处理生成正交矩阵q_t，并用正交矩阵q_t更新误差余量：其中是正交矩阵q_t的共轭转置矩阵；

(6)判断误差余量R_t的2范数是否小于1，若满足此收敛条件，则执行步骤(7)，否则将迭代次数t加1后重复步骤(3)～(5)；

(7)根据观测向量y和最新的列向量集合利用共轭梯度算法对雷达目标信号x进行恢复与重建：

7.1)对列向量集合进行共轭转置得到矩阵将该共轭转置矩阵分别与列向量集合和观测向量y相乘得到矩阵C和迭代向量b；

7.2)初始化迭代次数k为1，初始化残差向量r⁰和搜索方向d⁰为迭代向量b，初始化待恢复信号x⁰为零向量；

7.3)将矩阵C和搜索方向d^k-1相乘得到过渡向量q^k，同时将残差向量r^k-1进行共轭转置得到共轭转置向量(r^k-1)^H，将该共轭转置向量(r^k-1)^H与残差向量r^k-1相乘得到旧残差因子δ_o；

7.4)利用过渡向量q^k、搜索方向d^k-1和旧残差因子δ_o计算搜索步长α：α＝δ_o/[(d^k-1)^H·q^k]，其中(d^k-1)^H是搜索方向d^k-1的共轭转置；

7.5)利用搜索步长α、搜索方向d^k-1更新待恢复信号：x^k＝x^k-1+αd^k-1，其中x^k-1是第k-1次迭代得到的待恢复信号，同时利用搜索步长α和过渡向量q^k更新残差向量r^k：r^k＝r^k-1-αq^k，其中r^k-1是第k-1次迭代得到的残差向量；

7.6)将更新后的残差向量r^k进行共轭转置得到行残差向量(r^k)^H，再将该行残差向量(r^k)^H与残差向量r^k相乘得到新残差因子δ_n；

7.7)利用新残差因子δ_n、旧残差因子δ_o和更新后的残差向量r^k更新搜索方向d^k：d^k＝r^k+βd^k-1，以在下次迭代时更新搜索步长α和待恢复信号x^k，其中d^k-1是第k-1次迭代得到的搜索方向；

7.8)判断迭代次数k是否等于12，若k＝12，则迭代结束，输出待恢复信号x^k，否则将迭代次数k加1后重复步骤7.3)～7.7)。

本发明具有如下优点：

1.处理速度快，资源和功耗小

本发明基于FPGA实现，由于FPGA有丰富的内部资源，在数据处理过程中可调用FPGA内部的IP核实现数据的复乘、FFT运算等功能；且FPGA并行处理度高，可大幅提升数据的处理速度；另外FPGA还具有灵活的算法适应性和较低的功耗等优点，用FPGA作为稀疏信号重构算法的实现平台可减少系统资源和功耗。

2.优化资源，进一步提高了处理速度

在现有的OMP算法中，当字典矩阵规模较大时，用字典矩阵与误差余量相乘求相关性将耗费大量的时间和资源，本发明通过对误差余量进行补零并对补零后的向量进行IFFT运算的方法来计算字典矩阵与误差余量的相关性，减小了资源的占用，进一步提高了处理速度。

3.适用于大稀疏度信号的重建

本发明用共轭梯度算法CG对雷达目标信号进行恢复与重建，由于该CG算法的数据一致性很好，计算比较密集、规整，恢复大稀疏度信号时，CG算法与QR分解法相比具有更快的速度和更少的资源占用。

4.适用性强

现有OMP算法的FPGA实现大多采用的是定点型的实数，而雷达回波信号是带有相位信息的复数，使得算法在实际应用中受到限制。本发明采用单精度浮点型的复数，在雷达成像中具有更强的适用性。

附图说明

图1是本发明的实现流程图；

图2是在占用资源相同的情况下，用本发明和矩阵-向量相乘法求字典矩阵与误差余量相关性的时间对比图；

图3是用本发明和QR分解法恢复目标信号的时间对比图；

图4是用本发明和QR分解法恢复目标信号占用的DSP48E1s资源对比图。

具体实施方案

参照附图1，本发明的具体实施步骤如下：

步骤1.根据雷达回波数据产生观测向量和字典矩阵。

1.1)根据雷达回波数据利用Matlab软件产生一个2048×2048大小的傅里叶基矩阵Φ和512×2048大小的随机矩阵A，用该随机矩阵A分别与傅里叶基矩阵Φ和雷达回波数据相乘生成512×2048大小的字典矩阵Ψ和512长度的测量向量y：Ψ＝A·Φ，y＝A·x，其中x是雷达回波数据，包含2048个复数；再把字典矩阵Ψ和测量向量y中的数据量化成32bit单精度浮点型数据后存储到ROM中；

1.2)初始化误差余量R₀为观测向量y，并将误差余量R₀存储在RAM中，初始化列向量集合和位置索引集合Λ₀为空集：Λ₀＝φ，其中φ表示空集，初始化迭代次数t为1。

步骤2.计算字典矩阵与误差余量的相关性。

2.1)按随机抽取的先验知识对512长度的误差余量R_t-1补零生成2048长度的插值向量

2.2)调用FPGA内部的FFT核，对其进行配置，即将控制FFT/IFFT功能选择的fwd_inv信号设为0以实现IFFT功能，将数据类型选择为浮点型，数据输出方式选择自然顺序输出；

2.3)用配置好的FFT核对插值向量作IFFT运算，得到相关性向量：其中相关性向量g_t是长度为2048的复向量。

步骤3.寻找字典矩阵Ψ中与误差余量R_t-1相关性最大的列，并更新位置索引集合和列向量集合。

3.1)调用浮点乘法器IP核和浮点加法器IP核对相关性向量g_t中的各元素取模的平方，得到一个2048长度的实相关向量：h_t＝|g_t(i)|²，i＝1,2,...,2048，其中g_t(i)是相关性向量的g_t的第i个元素；

3.2)利用冒泡法找到实相关向量h_t中的最大元素在向量中的位置，记为λ_t，该位置λ_t即为字典矩阵Ψ中与误差余量R_t-1相关性最大的列在字典矩阵Ψ中所处的位置；

3.3)根据最新找到的位置λ_t更新位置索引集合Λ_t和列向量集合Λ_t＝Λ_t-1∪λ_t，其中Λ_t-1是第t-1次迭代后的位置索引集合，是第t-1次迭代得到的列向量集合，是字典矩阵Ψ中第λ_t列对应的列向量；并将得到的列向量集合存储在RAM中。

步骤4.更新误差余量R_t。

4.1)对列向量集合进行施密特正交化处理：

4.1a)进行第n次迭代时，列向量集合中共有n个长度为512的列向量，即其中φ_i是列向量集合的第i个列向量；

4.1b)初始化正交基矩阵η＝[η₁,η₂,...,η_i,...,η_n]，其中第i个正交基向量η_i是正交基矩阵η的第i个列向量，初始化第1个正交基向量η₁为列向量集合中的第1个列向量：η₁＝φ₁；

4.1c)调用FPGA内部的乘累加IP核和减法器IP核，按如下公式计算第i个正交基向量η_i：

${\mathrm{\&eta;}}_i={\mathrm{\&phi;}}_i-{\mathrm{\&Sigma;}}_{j=1}^{i-1}<{\mathrm{\&phi;}}_i,{\mathrm{\&kappa;}}_j>\&CenterDot;{\mathrm{\&kappa;}}_j,i=2,...,n$

其中，其中φ_i是列向量集合的第i个列向量，κ_j是施密特正交基向量，由以下公式计算：

κ_j＝η_j/||η_j||₂，j＝1,2,...,n-1

其中，||η_j||₂是第j个正交基向量η_j的2范数；

4.1d)将n个正交基向量η₁～η_n组合得到正交基矩阵η＝[η₁,η₂,...,η_i,...,η_n]，对该正交基矩阵η中的n个正交基向量η₁～η_n进行归一化处理得到正交矩阵q_t，该正交矩阵q_t即为对列向量集合进行施密特正交化处理得到的矩阵；

4.2)利用正交矩阵q_t更新误差余量R_t

4.2a)将正交矩阵q_t进行共轭转置得到转置正交矩阵调用FPGA内部的乘累加IP核将该转置正交矩阵与误差余量R_t-1相乘得到归一化误差矩阵：

4.2b)调用FPGA内部的乘累加IP核将正交矩阵q_t与归一化误差矩阵l_t进行乘累加操作得到减量误差向量：γ_t＝q_t·l_t，再调用两个浮点减法IP核将误差余量R_t-1与减量误差向量γ_t相减得到更新后的误差余量R_t。

步骤5.判断是否满足收敛条件

将更新后的误差余量R_t取共轭转置得到共轭误差向量调用乘累加核将该共轭误差向量与误差向量R_t进行乘累加操作，得到收敛因子σ，再判断该收敛因子σ是否小于1，若σ＜1，则进行步骤6，否则将迭代次数t加1后返回步骤2。

步骤6.用共轭梯度法对雷达目标信号进行恢复与重建。

6.1)对列向量集合进行共轭转置得到矩阵调用FPGA内的乘累加IP核将该共轭转置矩阵分别与列向量集合和观测向量y相乘，得到矩阵A和迭代向量b： $A={\widetilde{\mathrm{\&Phi;}}}_t^H\&CenterDot;{\widetilde{\mathrm{\&Phi;}}}_t,b={\widetilde{\mathrm{\&Phi;}}}_t^H\&CenterDot;y;$

6.2)初始化迭代次数k为1，初始化残差向量r⁰和搜索方向d⁰为迭代向量b，初始化待恢复信号x⁰为零向量；

6.3)调用乘累加IP核将矩阵C和搜索方向d^k-1相乘得到过渡向量q^k:q^k＝A·d^k-1，同时将残差向量r^k-1进行共轭转置得到共轭转置向量(r^k-1)^H，再调用乘累加IP核将该共轭转置向量(r^k-1)^H与残差向量r^k-1相乘得到旧残差因子：δ_o＝(r^k-1)^H·r^k-1；

6.4)调用乘累加IP核和浮点除法器IP核按照下面的公式来计算搜索步长α：

α＝δ_o/[(d^k-1)^H·q^k]

其中，δ_o是旧残差因子，(d^k-1)^H是搜索方向d^k-1的共轭转置，q^k是过渡向量；

6.5)调用浮点乘法器IP核和浮点加法器IP核，根据搜索步长α、搜索方向d^k-1和过渡向量q^k更新待恢复信号x^k和残差向量r^k：

x^k＝x^k-1+αd^k-1，r^k＝r^k-1-αq^k

其中x^k-1是第k-1次迭代得到的待恢复信号，r^k-1是第k-1次迭代得到的残差向量；

6.6)将更新后的残差向量r^k进行共轭转置得到行残差向量(r^k)^H，调用乘累加IP核将该行残差向量(r^k)^H与残差向量r^k相乘得到新残差因子：δ_n＝(r^k)^H·r^k；

6.7)调用浮点乘法器核和浮点加法器核，根据新残差因子δ_n、旧残差因子δ_o和更新后的残差向量r^k更新搜索方向d^k：

d^k＝r^k+βd^k-1

其中d^k-1是第k-1次迭代得到的搜索方向，更新搜索方向d^k是为了在下次迭代时更新搜索步长α和待恢复信号x^k；

6.8)判断迭代次数k是否等于12，若k＝12，则迭代结束，输出待恢复信号x^k，否则将迭代次数k加1后重复步骤6.3)～6.7)。

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

仿真一：

对用本发明和现有的矩阵—向量相乘方法求字典矩阵与误差余量相关性仿真在占用资源相同情况的时间，结果如图2所示：图2中，横坐标是测量向量长度，纵坐标是时间，其单位是μs。由图2可见，当测量向量长度较小时，矩阵—向量乘方法的速度稍快于本发明，但测量向量长度增大时，矩阵—向量乘方法的耗时远大于本发明，可见测量向量长度越大，本发明与矩阵—向量乘相比在速度方面的优优越性越高。

仿真二：

对用本发明和现有QR分解法仿真恢复原始信号在不同稀疏度下的时间，结果如图3所示：图3中，横轴是稀疏度，纵轴是时间，其单位是μs，从图3可见，当稀疏度k＜6时，用QR分解方法比本发明用时少，当稀疏度k≥6时，QR分解方法的用时大于本发明，且稀疏度越大，本发明的优势越明显。

仿真三：

对用本发明和现有QR分解法恢复原始信号仿真在不同稀疏度下的DSP48E1s资源占用情况，仿真结果如图4所示：图4中，横轴是稀疏度，纵轴是DSP48E1s资源占用量。从图4可见，当稀疏度k＜2时，QR分解法占用的DSP48E1s资源少于本发明，当稀疏度k≥2时，QR分解法占用的DSP48E1s资源多于本发明，且稀疏度越大，两者占用的DSP48E1s资源对比越明显，越能体现本发明的优势。

Claims (1)

1.一种基于FPGA的正交匹配追踪方法，包括如下步骤：

(1)根据雷达回波数据，利用Matlab软件产生一个2048×2048大小的傅里叶基矩阵Φ和512×2048大小的随机矩阵A，用该随即矩阵A分别对傅里叶基矩阵Φ和雷达回波数据进行抽取生成512×2048大小的字典矩阵Ψ和512长度的测量向量y，并将该字典矩阵Ψ和测量向量y中的数据量化成32bit单精度浮点型数据后存储到ROM中；

(2)初始化迭代次数t为1，初始化误差余量R₀为测量向量y，初始化列向量集合和位置索引集合Λ₀为空集；

(3)计算字典矩阵Ψ与误差余量R_t-1的相关性：即按随机抽取的先验知识对512长度的误差余量R_t-1补零生成2048长度的插值向量再对该插值向量进行IFFT运算得到相关性向量：

(4)对相关性向量g_t中各元素取模的平方得到2048长度的实相关向量h_t，求出实相关向量h_t中的最大元素所处的位置λ_t，该位置λ_t即为字典矩阵Ψ中与误差余量R_t-1相关性最大的列所处的位置，再利用该位置λ_t更新列向量集合和位置索引集合Λ_t；

(5)对列向量集合进行施密特正交化处理生成正交矩阵q_t，并用正交矩阵q_t更新误差余量：其中是正交矩阵q_t的共轭转置矩阵；

(6)判断误差余量R_t的2范数是否小于1，若满足此收敛条件，则执行步骤(7)，否则将迭代次数t加1后重复步骤(3)～(5)；

(7)根据观测向量y和最新的列向量集合利用共轭梯度算法对雷达目标信号x进行恢复与重建：

7.1)对列向量集合进行共轭转置得到矩阵将该共轭转置矩阵分别与列向量集合和观测向量y相乘得到矩阵C和迭代向量b；

7.2)初始化迭代次数k为1，初始化残差向量r⁰和搜索方向d⁰为迭代向量b，初始化待恢复信号x⁰为零向量；

7.3)将矩阵C和搜索方向d^k-1相乘得到过渡向量q^k，同时将残差向量r^k-1进行共轭转置得到共轭转置向量(r^k-1)^H，将该共轭转置向量(r^k-1)^H与残差向量r^k-1相乘得到旧残差因子δ_o；

7.4)利用过渡向量q^k、搜索方向d^k-1和旧残差因子δ_o计算搜索步长α：α＝δ_o/[(d^k-1)^H·q^k]，其中(d^k-1)^H是搜索方向d^k-1的共轭转置；

7.5)利用搜索步长α、搜索方向d^k-1更新待恢复信号：x^k＝x^k-1+αd^k-1，其中x^k-1是第k-1次迭代得到的待恢复信号，同时利用搜索步长α和过渡向量q^k更新残差向量r^k：r^k＝r^k-1-αq^k，其中r^k-1是第k-1次迭代得到的残差向量；

7.6)将更新后的残差向量r^k进行共轭转置得到行残差向量(r^k)^H，再将该行残差向量(r^k)^H与残差向量r^k相乘得到新残差因子δ_n；

7.7)利用新残差因子δ_n、旧残差因子δ_o和更新后的残差向量r^k更新搜索方向d^k：d^k＝r^k+βd^k-1，以在下次迭代时更新搜索步长α和待恢复信号x^k，其中d^k-1是第k-1次迭代得到的搜索方向；

7.8)判断迭代次数k是否等于12，若k＝12，则迭代结束，输出待恢复信号x^k，否则将迭代次数k加1后重复步骤7.3)～7.7)。