CN103178853B - 基于压缩感知的稀疏信号欠采样方法 - Google Patents

Abstract

基于压缩感知的稀疏信号欠采样方法，涉及一种稀疏信号欠采样方法。它是为了在保证信号恢复效果的情况下降低频域稀疏信号采样率。本发明采用FPGA产生触发信号和m序列，m序列调理后与被测稀疏信号混频，并经过低通滤波器滤波。数据采样模块检测到触发信号后对滤波后信号进行采样，存储采样数据。在信号重构时，需要先求得***的传递函数，以及和采样数据对应的m序列，然后应用OMP信号重构算法恢复出原信号。本发明适用于频域稀疏模拟信号的欠采样。

Description

基于压缩感知的稀疏信号欠采样方法

技术领域

本发明涉及一种稀疏信号欠采样方法。

背景技术

传统的信息采样过程必须遵循奈奎斯特采样定理，即采样速率至少要大于原信号最高频率的2倍，这样才能从采样得到的离散数据中不失真地恢复出原始信号。然而随着信息技术的发展，以奈奎斯特采样定理为基础的信号处理框架对前端ADC的采样速率和处理速度提出更高的要求，也给后端信息的传输、存储等环节带来巨大压力。解决这些压力常见的方案是信号压缩，但是，这种先采样后压缩的方法并没有减小前端ADC的压力，而且信号压缩意味着采样过程中有大量的冗余信息，浪费了大量的传感元、时间和存储空间等资源。

2004年由Candès和Donoho等人提出了压缩感知理论(CompressiveSensing，CS)。该理论能够将信号压缩和采样合并进行，即在信号获取的同时，就对数据进行适当地压缩，当信号具有稀疏性时就可以显著降低信号采样率。

在信号具有稀疏性的前提下，可以用一个与正交基矩阵Ψ不相关的观测矩阵Φ：M×N(M□N)(这里Φ的每一行可以看作是一个传感器，它与系数相乘，获取了信号的部分信息)，对信号x执行压缩观测：

y＝Φx(1)

就可以得到M个线性观测值(投影)y∈□^M。这些少量线性投影中则包含了重构信号x的足够信息。

然后通过特殊的算法从观测值y中得到原信号的重构信号

频域稀疏信号具有频谱稀疏的特性，即在频域上是稀疏的，这样正交基矩阵Ψ就可以取傅里叶变换矩阵。设x(t)是一模拟信号，对它进行傅里叶变换，即把它用N×1维的傅里叶变换基向量的线性组合来表示。则x(t)可以展开为：

$x\left(t\right)=\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\α}}_n{\mathrm{\ψ}}_n\left(t\right)---\left(2\right)$

其中：ψ_n(t)＝e^j2πt(n-1)/N。

傅里叶变换系数：

α_n＝<x(t),ψ_n(t)>＝ψ_n ^T(t)x(t)

这些系数组成系数向量α＝(α₁,α₂,…,α_n)^T，这就是信号的频谱，每个元素都是信号的一个频谱分量。如果将系数向量α中的元素按降序排列后，元素值迅速衰减，或者系数向量α中值较大的系数个数为K，并且K比N小很多，则表明该信号在频域是稀疏的，称x(t)是频域稀疏信号，其稀疏度为K。

对于频域稀疏信号采样，利用压缩感知理论可以突破乃奎斯特采样定律的限制，实现欠采样。压缩感知理论表示，频域稀疏信号的采样速率不再取决于信号的最高频率或带宽，而是取决于信号的稀疏度K。设信号稀疏度为K，当信号为离散信号时，由公式(3)可得到最少采样点数M的范围，其中N为原信号的长度，c为一个很小的常数。当信号为连续模拟信号时，N为原信号的乃奎斯特频率，由公式(3)可得到最低采样率M的范围。

M≥cKlog(N/K)(3)

发明内容

本发明是为了在保证信号恢复效果的情况下降低频域稀疏信号采样率，从而提供一种基于压缩感知的稀疏信号欠采样方法。

基于压缩感知的稀疏信号欠采样方法，它由以下步骤实现：

步骤一、采用嵌入在FPGA中m序列发生器产生m序列；并采用FPGA同步产生触发信号；

步骤二、将步骤一产生的m序列采用信号调理电路进行信号调理，获得调理后的m序列；

步骤三、将步骤二获得的调理后m序列与被测稀疏信号采用乘法器进行混频，获得混频信号；

步骤四、采用低通滤波器对步骤三获得的混频信号进行低通滤波，获得低通滤波后信号；

步骤五、采用步骤一中产生触发信号触发采样电路，并采用采样电路对步骤四获得的低通滤波后信号进行采样，获得采样结果；

步骤六、采用上位机对步骤五获得的采样结果进行信号重构，获得原始被测稀疏信号的重构结果。

步骤六中所述的采用上位机对步骤五获得的采样结果进行信号重构的具体方法是：

步骤A、根据公式：

H(s)＝H₁(s)G(s)

获得混频到采样之间的传递函数H(s)；

式中：H₁(s)为低通滤波器的传递函数，所述低通滤波器为二阶巴特沃斯模拟低通滤波器；

G(s)为混频到采样之间各级运算放大电路传递函数的连乘积；

H₁(s)是根据公式：

$H_1\left(s\right)=\frac{w_c^2}{s^2+\frac{w_c}{Q}s+w_c^2}$

获得的；

式中：w_c是低通滤波器的截止频率，Q是品质因数；s是符号变量；

G(s)是根据公式：

G(s)＝G₁(s)G₂(s)…G_n(s)

获得的；

式中：G₁(s)是第一个运算放大电路传递函数；G₂(s)是第二个运算放大电路的传递函数；G_n(s)是第n个运算放大电路的传递函数；n为正整数；

步骤B、采用双线性变换法对步骤A获得的混频到采样之间的传递函数H(s)进行离散化处理，获得数字滤波器传递函数H(z)，并求得数字滤波器传递函数H(z)的脉冲响应h(n)；

步骤C、对步骤B中的数字滤波器传递函数H(z)的脉冲响应h(n)进行翻转，获得翻转结果h_r(n)；

步骤D、将步骤一中的m序列p_c(n)和步骤C获得的翻转结果h_r(n)做卷积运算，获得观测矩阵Φ；

步骤E、根据公式：

Θ＝ΦΨ

获得感知矩阵Θ；

式中：Ψ为正交基矩阵，所述Ψ根据公式：

Ψ＝dftmtx(N)

获得；所述dftmtx为MATLAB软件的自带函数；

步骤F、根据步骤E获得的感知矩阵Θ通过正交匹配追踪算法实现信号重构。

基于压缩感知的稀疏信号欠采样方法的实现装置，它包括FPGA15、调理电路16、乘法器11、低通滤波器12、采样电路13和上位机14；

所述FPGA15内嵌有m序列发生器，所述m序列发生器用于产生m序列；

所述FPGA15的m序列输出端与信号调理电路16的信号输入端连接；所述信号调理电路16的信号输出端与乘法器11的一号信号输入端连接；所述乘法器11的二号信号输入端用于接收被测稀疏信号；所述乘法器11的信号输出端与低通滤波器12的信号输入端连接；所述低通滤波器12的信号输出端与采样电路13的信号输入端连接；所述采样电路13的信号输出端与上位机14的信号输入端连接。

本发明的方法利用频域稀疏信号频谱稀疏的特性，在保证信号恢复效果的前提下，大幅度降低了频域稀疏信号的采样率，且当信号所含频率分量的最大值越大时，该效果越明显，能够使实际采样率远小于信号的乃奎斯特频率。本发明的能够降低对AD器件的速率的要求，且采样数据量小，节省存储空间，有利于数据的传输。

附图说明

图1是本发明的***结构示意图；图2是FPGA中嵌入的m序列发生器和控制逻辑示意图；图3是具体实施方式五中减法电路的连接示意图；图4是具体实施方式五中AD633模拟乘法器的电路连接示意图；图5是具体实施方式五中MAX275内部结构和二阶巴特沃斯模拟低通滤波器电路连接示意图；图6是原始被测稀疏信号的数据点连线波形示意图；图7是含3kHz、4kHz两个频点的原始被测稀疏信号的频谱示意图；图8是重构后信号中数据点的连线波形示意图；图9是用随机解调***的数据重构含3kHz、4kHz两个频点的原始被测稀疏信号的实验结果图。

具体实施方式

具体实施方式一、基于压缩感知的稀疏信号欠采样方法，它由以下步骤实现：

步骤一、采用嵌入在FPGA中m序列发生器产生m序列；并采用FPGA同步产生触发信号；

步骤二、将步骤一产生的m序列采用信号调理电路进行信号调理，获得调理后的m序列；

步骤三、将步骤二获得的调理后m序列与被测稀疏信号采用乘法器进行混频，获得混频信号；

步骤四、采用低通滤波器对步骤三获得的混频信号进行低通滤波，获得低通滤波后信号；

步骤五、采用步骤一中产生触发信号触发采样电路，并采用采样电路对步骤四获得的低通滤波后信号进行采样，获得采样结果；

步骤六、采用上位机对步骤五获得的采样结果进行信号重构，获得原始被测稀疏信号的重构结果。

具体实施方式二、基于压缩感知的稀疏信号欠采样方法的实现装置，它包括FPGA15、调理电路16、乘法器11、低通滤波器12、采样电路13和上位机14；

所述FPGA15内嵌有m序列发生器，所述m序列发生器用于产生m序列；

所述FPGA15的m序列输出端与信号调理电路16的信号输入端连接；所述信号调理电路16的信号输出端与乘法器11的一号信号输入端连接；所述乘法器11的二号信号输入端用于接收被测稀疏信号；所述乘法器11的信号输出端与低通滤波器12的信号输入端连接；所述低通滤波器12的信号输出端与采样电路13的信号输入端连接；所述采样电路13的信号输出端与上位机14的信号输入端连接。

具体实施方式三、本具体实施方式与具体实施方式二所述的基于压缩感知的稀疏信号欠采样方法的实现装置的区别在于，FPGA15的型号为EP2C8Q208。

具体实施方式四、本具体实施方式与具体实施方式二所述的基于压缩感知的稀疏信号欠采样方法的实现装置的区别在于，乘法器11采用型号为AD633的芯片实现。

具体实施方式五、本具体实施方式与具体实施方式二所述的基于压缩感知的稀疏信号欠采样方法的实现装置的区别在于，低通滤波器12为采用型号为MAX275的芯片搭建的二阶巴特沃斯模拟低通滤波器。

工作原理：准备工作：

***运行前的准备工作是获得某一组初值下的m序列值。操作步骤如下：

1、给m序列发生器设定一组初值，并保持不变，这样每次重新产生m序列时都按同样的规律，保证每次的产生的序列都一样。

2、将均匀采样模块的数据输入端连接到FPGA的m序列输出端，采样模块设置为触发采样模式，触发控制端连接到FPGA的触发信号输出端。数据采集模块的采样率设置为和m序列的跳变频率(即m序列发生器的时钟频率)相同，例如，m序列的时钟是10kHz，则采样率也设为10kHz。

3、使控制触发模块产生触发信号，使m序列模块产生m序列，同时触发信号被数据采集模块检测到，开始对m序列进行采样。采样一定时间后得到这组初值下完整的m序列值，将其保存到一个文件中备用。

实际运行过程：

1、给m序列发生器设定一组与准备工作时相同的初值，并保持不变。

2、将被测信号连接到***的输入端，即乘法器的另一个输入端。

3、将均匀采样模块的数据输入端连接到滤波器的输出端。采样模块设置为触发采样模式，触发控制端还是连接到FPGA的触发信号输出端。采样率设置的只要比由公式(3)计算出来的值大即可。例如，被测信号的稀疏度为2，信号最高频率分量为4kHz，则计算得到采样率需大于6.6c。其中c是一个常数，和实际***本身的误差大小、准确度、工作环境等有关，需要实验来确定。一般情况下，采样率大一些重构效果更好。

4、控制触发模块产生触发信号，使m序列模块产生m序列，同时触发信号被数据采集模块检测到，开始对滤波后的信号进行采样。将采样数据传输到上位机中存储。

5、数据采集模块采集到的数据是被测信号和m序列经过乘法器混频，低通滤波器滤波以及调理后的结果。由于采样速率已知，通过计算即可得知这段采样时间里参与***作用的m序列的长度，从保存好的m序列文件中取出对应长度的一段即可。例如，采样速率为2kHz，截取0～0.1s的采样数据，即开始的200个数据点。若m序列的跳变频率是10kHz，则0～0.1s对应的数据量是1000，从之前保存好的m序列文件中截取最开始的1000个数据即可。在上位机中调用编好的信号重构算法，利用截取的采样数据和m序列值对信号进行重构，得到原信号的频点位置和对应的幅值。如果有必要，则将重构出来的数据经过D/A变为模拟信号。

***各模块详细说明：

FPGA模块用于产生伪随机序列，同时控制触发采样信号的产生。

本发明在具体应用时采用的是m序列这种伪随机序列，利用Altera公司的CycloneII系列EP2C8Q208这款FPGA产生。m序列是最长线性反馈移位寄存器序列的简称，是由带线性反馈的移位寄存器产生的序列，并且具有最长的周期，具有容易产生、规律性强等许多优良的特性，是一种重要的伪随机序列。

FPGA中还实现了触发采样信号的控制逻辑，用于同步产生触发采样信号，控制后端数据采集模块对末端输出信号进行采样。

图2是在FPGA中实现的一个10级m序列发生器和触发信号控制逻辑示意图。

m序列最初的幅值范围为0～3.3V，为了产生负的幅值，要用减法器进行调整，减法电路如图3所示，使其从原来的0～3.3V变为-1.6～+1.7V，同时也大大减小了m序列最初带有的较大毛刺，提高了质量。

乘法器：m序列经过调理之后与频域稀疏信号在模拟乘法器中相乘，实现混频。信号在时域相乘，在频域表现为频谱的卷积。由于m序列的频谱范围很广，通过卷积就对频域稀疏信号进行了调制，实现了扩频，将原信号的频谱涂抹到整个频率轴，频谱上每个位置处都含有原信号的所有信息，并且，由于m序列的参与，相当于对原信号进行了信号加扰，频谱上每一点处的信息都具有了独特的标识。

本发明在具体应用时采用的是AD633模拟乘法器芯片，电路如图5所示。AD633内部对相乘后的信号有1/10的衰减，信号幅值太小会给后面的分析处理带来不便，为了弥补这一衰减，对乘法器的输出信号进行了放大。运放电路可使用固定增益放大器设计，以要保证对不同频率的信号分量增益都相同。

低通滤波器：这部分的关键是尽量使实际滤波器的特性与理论计算的传递函数相符。

本发明具体实现时采用有源低通滤波器，有源滤波器可采用放大器搭建，也可以使用集成的滤波器芯片来设计。采用集成的滤波器芯片一般可以获得指标更好、参数更精确的滤波器，且***电路简单，设计方便，在该***中选用了MAX275有源滤波器芯片搭建二阶巴特沃斯模拟低通滤波器，如图5所示。

采样电路：用于对滤波后的信号进行低速采样并存储。

该部分的关键是保证采样开始时刻与m序列经过***到达采样模块的时刻同步，即采样的起始时刻要合适。采样触发信号由FPGA产生，输入给采样模块，只有收到触发采样信号后，采样模块才开始进行数据采集和存储。

上位机14：信号重构需要先对观测过程进行建模，求得各个模块的传递函数，从而求得对应的观测矩阵Φ，然后求得感知矩阵Θ。再利用Θ以及低速均匀采样得到的数据通过信号重构算法恢复原信号。

具体步骤为：

1、通过计算推导或实际采样的方式得到一段采样时间内的数据和参与***作用的m序列p_c(n)；由于采样速率已知，通过计算即可得知这段采样时间里参与***作用的m序列的长度，从保存好的m序列文件中取出对应长度的一段。例如，采样速率为2kHz，截取0～0.1s的采样数据，即开始的200个数据点。若m序列的跳变频率是10kHz，则0～0.1s对应的数据量是1000，从之前保存好的m序列文件中截取最开始的1000个数据即可。

2、求混频模块之后的各级调理电路的传递函数。

混频模块、滤波器模块后有放大调理模块，计算出这些模块的传递函数G₁(s)，G₂(s)…；然后，求这些调理模块的总传递函数G(s)，G(s)为各级传递函数的连乘积：

G(s)＝G₁(s)G₂(s)…(4)

如果这些放大调理模块都是固定增益放大，则求得的G(s)为常数。

3、求模拟低通滤波器的传递函数。

二阶巴特沃斯模拟滤波器的传递函数为：

$H_1\left(s\right)=\frac{w_c^2}{s^2+\frac{w_c}{Q}s+w_c^2}---\left(5\right)$

其中w_c是滤波器的截止频率，Q是品质因数。

4、求混频模块之后到采样模块之间的传递函数H(s)，即放大调理环节和低通滤波环节的总传递函数，并求出对应的单位取样响应h(n)。

放大调理环节和低通滤波环节的总传递函数为：

H(s)＝H₁(s)G(s)(6)

在matlab中通过c2d()函数用双线性变换法将H(s)离散化，变换为数字滤波器的传递函数H(z)，离散频率可取m序列的跳变频率。再利用impz()函数求得数字滤波器的脉冲响应h(n)。

5、求观测矩阵Φ，正交字典矩阵Ψ和感知矩阵Θ

观测矩阵Φ是由伪随机序列p_c(n)和单位取样响应h(n)卷积而得，由于信号在频域稀疏，故正交字典矩阵Ψ为傅里叶变换矩阵。

求感知矩阵Θ的流程如下：

步骤1：输入单位取样响应h(n)；

步骤2：对h(n)翻转，这里是对h(n)进行倒序，求得h_r(n)，即：h(1)＝h_r(N),…,h(N)＝h_r(1)。

步骤3：用伪随机序列p_c(n)和h(n)卷积，即p_c(n)和h_r(n)依次移位、相乘、求和，求得M×N的Φ。每次移位的大小为T_s＝N/M，N为原信号的长度，也等于伪随机序列p_c(n)的长度，M为用downsample函数对滤波后信号低速采样时采到的数据个数。

forj＝1:M

ifj×T_s≤h(n)的长度

Φ的第j行，第1列至第j×T列为：

Φ(j,1:j×T_s)＝p_c(1:j×T_s).×h_r((n-j×T_s+1):n)(7)

else

Φ的第j行，第j×T_s-n+1列至第j×T_s列为：

Φ(j,(j×T_s-n+1):(j×T_s))＝p_c((j×T_s-n+1):j×T_s).×h_r(1:n)(8)

end

end

步骤4：求正交基矩阵Ψ。调用MATLAB自带的函数dftmtx可以求得N×N的傅里叶变换矩阵Ψ：Ψ＝dftmtx(N)；

步骤5：求感知矩阵Θ。根据Θ＝ΦΨ，将求得的观测矩阵Φ与傅里叶变换矩阵Ψ相乘，得到M×N的感知矩阵Θ。

6、用OMP信号重构算法重构原信号：压缩感知理论有多种信号重构算法，正交匹配追踪(OrthogonalMatchingPursuit，OMP)算法的步骤如下：

步骤1：设定算法输入，感知矩阵Θ，观测信号y，稀疏度K。

步骤2：初始化各参数，重构信号残差r₀＝y，信号支撑集

步骤3：迭代，在第l次循环(l≥1)，运行以下(a)～(d)步。

(a)利用相关计算寻找信号支撑索引：

λ_l＝argmax_j＝1,…,N|<r_l-1,θ_j>|(9)

(b)将寻找到的信号支撑索引加入信号支撑集：

Λ_l＝Λ_l-1∪{λ_l}(10)

(c)更新残差：

$r_l=y-{\mathrm{\&theta;}}_{{\mathrm{\&Lambda;}}_l}\left({\mathrm{\&theta;}}_{{\mathrm{\&Lambda;}}_l}^{+}y\right)---\left(11\right)$

(d)若l≥K，算法结束。

步骤4：输出重构信号

$\hat{x}={\mathrm{\&theta;}}_{{\mathrm{\&Lambda;}}_l}^{+}y,{\hat{x}}_{\{1,\mathrm{...},N\}-{\mathrm{\&Lambda;}}_l}=0---\left(12\right)$

以下采用具体的仿真实验验证本发明的效果：

FPGA产生跳变频率为10kHz的m序列，用信号发生器产生稀疏度为2的测试信号，频率分量为3kHz和4kHz。经混频、低通滤波、采样和信号重构后，原信号信号含有3kHz和4kHz两个频点，数据点连线形成的波形如图6所示，频谱如图7所示；所示。重构信号数据点连线波形如图8所示，重构信号频谱如图9所示。

由上述仿真实验可见：

1、本发明利用频域稀疏信号频谱稀疏的特性，在保证信号恢复效果的前提下，降低了频域稀疏信号的采样率，当信号所含频率分量的最大值越大时，这种效果越明显，能够使实际采样率远小于信号的乃奎斯特频率。

2、能够降低对AD器件的速率的要求，采样数据量小，节省存储空间，降低成本，同时有利于数据的传输。

Claims (1)

1.基于压缩感知的稀疏信号欠采样方法，它由以下步骤实现：

步骤一、采用嵌入在FPGA中m序列发生器产生m序列；并采用FPGA同步产生触发信号；

步骤二、将步骤一产生的m序列采用信号调理电路进行信号调理，获得调理后的m序列；

步骤三、将步骤二获得的调理后m序列与被测稀疏信号采用乘法器进行混频，获得混频信号；

步骤四、采用低通滤波器对步骤三获得的混频信号进行低通滤波，获得低通滤波后信号；

步骤五、采用步骤一中产生触发信号触发采样电路，并采用采样电路对步骤四获得的低通滤波后信号进行采样，获得采样结果；

步骤六、采用上位机对步骤五获得的采样结果进行信号重构，获得原始被测稀疏信号的重构结果；

其特征在于，步骤六中所述的采用上位机对步骤五获得的采样结果进行信号重构的具体方法是：

步骤A、根据公式：

H(s)＝H₁(s)G(s)

获得混频到采样之间的传递函数H(s)；

式中：H₁(s)为低通滤波器的传递函数，所述低通滤波器为二阶巴特沃斯模拟低通滤波器；

G(s)为混频到采样之间各级运算放大电路传递函数的连乘积；

H₁(s)是根据公式：

$H_1\left(s\right)=\frac{w_c^2}{s^2+\frac{w_c}{Q}s+w_c^2}$

获得的；

式中：w_c是低通滤波器的截止频率，Q是品质因数；s是符号变量；

G(s)是根据公式：

G(s)＝G₁(s)G₂(s)…G_n(s)

获得的；

式中：G₁(s)是第一个运算放大电路传递函数；G₂(s)是第二个运算放大电路的传递函数；G_n(s)是第n个运算放大电路的传递函数；n为正整数；

步骤B、采用双线性变换法对步骤A获得的混频到采样之间的传递函数H(s)进行离散化处理，获得数字滤波器传递函数H(z)，并求得数字滤波器传递函数H(z)的脉冲响应h(n)；

步骤C、对步骤B中的数字滤波器传递函数H(z)的脉冲响应h(n)进行翻转，获得翻转结果h_r(n)；

步骤D、将步骤一中的m序列p_c(n)和步骤C获得的翻转结果h_r(n)做卷积运算，获得观测矩阵Φ；

步骤E、根据公式：

Θ＝ΦΨ

获得感知矩阵Θ；

式中：Ψ为正交基矩阵，所述Ψ根据公式：

Ψ＝dftmtx(N)

获得；所述dftmtx为MATLAB软件的自带函数；

步骤F、根据步骤E获得的感知矩阵Θ通过正交匹配追踪算法实现信号重构。