基于近似计算的正交匹配追踪重构方法及系统与流程

本发明涉及信号处理领域，具体涉及一种基于近似计算的正交匹配追踪重构方法及系统。

背景技术：

压缩感知(compressedsensing，cs)是近十多年来信号处理领域最著名的研究进展之一，主要由压缩采样和重构两个步骤组成。压缩感知中的采样过程相对简单，但重构非常复杂，所以对于重构算法的研究是压缩感知中非常重要的一方面。由于算法软件实现速度慢，功耗大，不能达到实时重建的效果，部分研究者开始利用fpga(fieldprogrammablegatearray)技术对算法进行加速设计。正交匹配追踪(orthogonalmatchingpursuit，omp)算法是应用非常广泛的重构算法之一，且相对简单，便于硬件实现。但是，omp算法依然包含大量的内积、比较和矩阵求逆运算。特别是最小二乘(leastsquares，ls)计算，其涉及到的矩阵求逆运算的计算复杂度高且不利于硬件实现。

技术实现要素：

有鉴于此，本发明的目的在于提供一种基于近似计算的正交匹配追踪重构方法及系统，能够通过减少最小二乘运算的次数，进而减少矩阵求逆的次数，有效降低算法的计算复杂度，提高信号重构效率。

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

一种基于近似计算的正交匹配追踪重构方法，包括以下步骤：

步骤s1：采集原始信号，并进行预处理；

步骤s2:将预处理后的输入信号，基于压缩感知理论，映射成低维的测量信号；

步骤s3:采用基于近似计算的正交匹配追踪重构算法，对测量信号进行重构，进一步得到的重构信号。

进一步的，所述原始信号是由两个不同频率的正弦波相加得到的。

进一步的，所述步骤s2具体为:

设有原始信号f∈rⁿ和测量矩阵φ∈r^m×n，原始信号f通过稀疏基ψ∈r^n×n变换成一个k稀疏信号x∈rⁿ，则压缩感知中的采样过程用下式表示：

y＝φf＝φψx＝ax(1)

其中a＝φψ被称作cs矩阵,φ是随机解调架构矩阵,ψ是傅里叶变换基。

进一步的，所述步骤s3具体为：

步骤s31：残差rt-1与cs矩阵a的列向量φj进行内积运算，找到内积值最大的列的列索引λt，并根据最大列索引λt更新原子集

步骤s32:利用近似计算方法替代最小二乘计算

步骤s33:更新残差更新迭代次数t＝t+1，判断若迭代次数t＜k，则执行步骤s31，否则执行下一步；

步骤s34:输出最终信号估计值

进一步的，所述步骤s32具体为：

步骤(1)：计算残差rt-1与cs矩阵a的列向量φj内积中的最大值并记录所对应的列索引，即残差rt-1初始化为观测向量y；

步骤(2)：根据最大列索引λt更新原子集at是根据每次迭代的最大列索引λt从cs矩阵a中提取对应的最大列组成的新的矩阵。

进一步的，所述步骤s32具体为：

对于最小二乘解，暂不提取最大列进行计算：

将a＝φψ代入(2)式得：

对于(3)式中的ψ^-1用ψ^t去近似替代，(ψ^t)^-1ψ^t的结果是单位矩阵直接化简掉，于是(3)式变成了下式：

由(4)式求出的是一个1024*1维的列向量，求出中对应最大列序号位置上的值，进一步化简

表示由

中对应最大列序号位置上的值组成的向量；表示ψ^t中对应最大列序号的列组成的矩阵；

用来表示近似计算的最小二乘解：

进一步的，所述步骤s3在最后一次迭代求解最小二乘时，采用改进的cholesky分解算法进行矩阵求逆的运算：

将矩阵进行如下分解：

c＝ldl^t(7)

其中l是对角线元素为1的下三角矩阵，d是对角矩阵，分别按以下公式进行计算：

然后c^-1按以下公式计算获得：

c^-1＝(l^-1)^td^-1l^-1(10)

其中矩阵d^-1可以通过将矩阵d对角线上的元素求倒数获得，而矩阵l^-1可以由下式获得:

在第k次迭代时用正常的最小二乘步骤进行求解：

一种基于近似计算的正交匹配追踪重构系统，包括依次连接的计算模块、存储模块和控制模块；所述计算模块完成算法迭代过程中的所有运算；所述存储模块用来存储cs矩阵，稀疏基，观测向量和所有的中间矩阵，向量和数据；所述控制模块控制算法的各个状态的转换以及从存储模块中读取相应的数据到计算模块中和将结果写回到存储模块中。

进一步的，所述计算模块包括矩阵内积单元，近似计算单元，残差计算单元，c矩阵求逆单元和最小二乘计算单元。

进一步的，所述矩阵内积单元由256个复数乘法器构成的并行乘法器和255个复数加法器构成的加法器树和一个比较器组成。

本发明与现有技术相比具有以下有益效果：

本发明能够通过减少最小二乘运算的次数，进而减少矩阵求逆的次数，有效降低算法的计算复杂度，提高信号重构效率，更适合于硬件加速实现，达到实时重建的效果。

附图说明

图1是本发明一实施例中信号重构方法流程图；

图2是本发明一实施例中整体硬件架构图；

图3是本发明一实施例中矩阵内积单元硬件架构图；

图4是本发明一实施例中矩阵分解元素计算硬件架构图；

图5是本发明一实施例中牛顿-拉夫森求倒数硬件架构图

图6是本发明一实施例中系统状态转换图；

图7是本发明一实施例中重构结果图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

请参照图1，本发明提供一种基于近似计算的正交匹配追踪重构方法，包括以下步骤：

步骤s1：采集原始信号，并进行预处理；

步骤s2:将预处理后的输入信号，基于压缩感知理论，映射成低维的测量信号；

设有原始信号f∈rⁿ和测量矩阵φ∈r^m×n，原始信号f通过稀疏基ψ∈r^n×n变换成一个k稀疏信号x∈rⁿ，则压缩感知中的采样过程用下式表示：

y＝φf＝φψx＝ax(1)

其中a＝φψ被称作cs矩阵,φ是随机解调架构矩阵,ψ是傅里叶变换基。

步骤s3:采用基于近似计算的正交匹配追踪重构算法，对测量信号进行重构，进一步得到的重构信号。

在本实施例中，所述步骤s3具体为：

步骤s31：残差rt-1与cs矩阵a的列向量φj进行内积运算，找到内积值最大的列的列索引λt，并根据最大列索引λt更新原子集

步骤(1)：计算残差rt-1与cs矩阵a的列向量φj内积中的最大值并记录所对应的列索引，即残差rt-1初始化为观测向量y；

步骤(2)：根据最大列索引λt更新原子集at是根据每次迭代的最大列索引λt从cs矩阵a中提取对应的最大列组成的新的矩阵。

步骤s32:利用近似计算方法替代最小二乘计算

对于最小二乘解，暂不提取最大列进行计算：

将a＝φψ代入(2)式得：

对于(3)式中的ψ^-1用ψ^t去近似替代，(ψ^t)^-1ψ^t的结果是单位矩阵直接化简掉，于是(3)式变成了下式：

由(4)式求出的是一个1024*1维的列向量，求出中对应最大列序号位置上的值，进一步化简

表示由中对应最大列序号位置上的值组成的向量；表示ψ^t中对应最大列序号的列组成的矩阵；

用来表示近似计算的最小二乘解：

步骤s33:更新残差更新迭代次数t＝t+1，判断若迭代次数t＜k，则执行步骤s31，否则执行下一步；

步骤s34:输出最终信号估计值

优选的，在本实施例中，所述步骤s3在最后一次迭代求解最小二乘时，采用改进的cholesky分解算法进行矩阵求逆的运算：

将矩阵进行如下分解：

c＝ldl^t(7)

其中l是对角线元素为1的下三角矩阵，d是对角矩阵，分别按以下公式进行计算：

然后c^-1按以下公式计算获得：

c^-1＝(l^-1)^td^-1l^-1(10)

其中矩阵d^-1可以通过将矩阵d对角线上的元素求倒数获得，而矩阵l^-1可以由下式获得:

在第k次迭代时用正常的最小二乘步骤进行求解：

在本实施例中，如图2所示，提供一种基于近似计算的正交匹配追踪重构系统，包括依次连接的计算模块、存储模块和控制模块；所述计算模块完成算法迭代过程中的所有运算；所述存储模块用来存储cs矩阵，稀疏基，观测向量和所有的中间矩阵，向量和数据；所述控制模块控制算法的各个状态的转换以及从存储模块中读取相应的数据到计算模块中和将结果写回到存储模块中。

优选的，在本实施例中，所述计算模块包括矩阵内积单元，近似计算单元，残差计算单元，c矩阵求逆单元和最小二乘计算单元。

参考图3，在本实施例中，所述矩阵内积单元由256个复数乘法器构成的并行乘法器和255个复数加法器构成的加法器树和一个比较器组成。一个复数乘法器通过调用4个dsp资源实现。加法器通过调用slices资源实现。每个周期从存储模块中读取a的一列和y或(r)进行内积操作，结果送入到比较器内进行比较，并记录较大列的列序号。最终在比较完所有列的内积后，输出内积值最大的列的列序号。近似计算单元，残差计算单元和最小二乘计算单元的大部分运算都可以复用这里的并行乘法器和加法器树结构实现。

在本实施例中，c矩阵求逆单元具体的硬件实现架构如图4所示；因为矩阵c是正定对称矩阵，所以c^-1也是对称矩阵，所以只需要求出c^-1对角线上的元素和下三角部分的元素就可以了。首先需要根据c矩阵用公式(8)和(9)求出l矩阵和d矩阵的元素。由公式(8)和(9)可知l矩阵和d矩阵的元素具有相关性，所以只能按照一定的顺序去计算。然后在求矩阵d^-1的元素时会涉及到除法，本文选用的是比较经典的牛顿-拉夫森方法去求倒数。这种方法是找到一个在零点有解x＝1/m的函数f(x)，一般这个函数就是下式：

f(x)＝1/x-m(13)

这个函数在零点的解就是m的倒数。牛顿-拉夫森方法给出了迭代公式：

本实施例中使用初始值x0＝3-2m，然后利用公式(14)进行迭代计算，xi的值就会收敛于m的倒数。用这种方法将除法运算转化成了乘法和减法运算。具体的硬件实现架构如图5所示。

然后利用公式(11)求l^-1矩阵的元素，同样可以利用图5中的架构实现。最后可以根据公式(10)计算矩阵c^-1。

在本实施例中，存储模块主要用来存储算法计算所需要的数据，包括cs矩阵，观测向量，稀疏基，残差等和从计算模块写回的各种数据。各个数据的存储具体情况如表1所示。为了能够在一个周期读出矩阵一列的数据并行的送入到并行乘法器中，本发明将矩阵一列的数据存在bram同一个地址中。其中cs矩阵a用单独的一块bram存储，地址是0到1023，傅里叶变换基ψ用单独的一块bram存储，地址是0到1023，观测向量y和残差r存储到相同的一块bram中，观测向量存到0地址，残差r按迭代顺序存在y向量后面。

表1内存资源使用情况

在本实施例中，控制模块主要通过状态机完成对各个运算子模块工作条件的控制。状态机的状态转换图如6所示，各状态说明如表2所示。每个子模块都有一个控制信号，在子模块需要工作的时候状态机将该子模块的控制信号置1，当任务完成时，子模块反馈给控制模块一个信号，然后状态机将该子模块的控制信号置0，进入下一个状态。当迭代次数t＜k时，状态机在state0到state4之间循环，当t＝k时，从state4跳到state5然后经过state6，state7和state8结束循环，系统复位到state0。

表2状态机的各状态说明

实施例1：

本实施例中，输入原始信号f由四个频率分别为50hz，70hz，90hz，110hz，幅值分别0.3，0.4，0.5，0.6的正弦波相加得到，采样序列1:1024，采样频率为1024；2)在maltab中运行acomp程序，获得cs矩阵a，傅里叶变换基ψ，观测向量y的原始数据，并对数据进行定点量化处理；3)将得到的数据导入coe文件，在xilinx公司的开发工具vivado2016.2中进行仿真测试；4)将fpga中仿真得到的结果导入matlab中，在matlab中进行反傅里叶变换得到重构信号，并与原始信号进行对比，验证是否重构成功。

acomp算法在fpga中输出原始信号在频域的值，是放大了2^9倍的定点数据，分别为169-140i，167+139i，148-105i，150+103i，62-138i,59+140i,35+103i,36-102i对应的在matlab中仿真的数据是176-140i，181+139i，163-111i，158+111i,67-135i,63+136i,34+107i,35-111i。图7是k＝8时的仿真结果，分别为原始信号和根据acomp算法在fpga中的输出结果重构的信号,重构信噪比为22.8537db。与现有技术中的重构速度进行对比，在相同的133.33mhz频率下，当k＝8时本发明重构时间为66.904us，现有技术为68.55us，当k＝36时本发明重构时间为296.14us，现有技术为327us。可见本发明有效的提高了加速效果，并且在稀疏度较高时加速效果更加明显。

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。