CN104052494B - 面向频域稀疏信号的信号重构方法 - Google Patents

Abstract

面向频域稀疏信号的信号重构方法，涉及一种频域稀疏信号重构方法，属于信号处理技术领域。本发明为了解决现有的方法计算量大、耗费时间长，所需数据存储空间大、功耗大，实时性不高、效率不高的问题。提出了面向频域稀疏信号的信号重构方法，该方法利用离散傅里叶变换矩阵的对称性和矩阵运算的一些特性，对压缩感知模型中计算感知矩阵的方式进行了改进，设计出一种新的感知矩阵计算方法，然后引入FFT变换，从而实现了计算量的减小和速度的提升。本发明可应用于信号处理领域。

Description

面向频域稀疏信号的信号重构方法

技术领域

本发明涉及一种基于FFT的压缩矩阵构造方法，尤其涉及面向频域稀疏信号的信号重构方法，属于信号处理技术领域。

背景技术

传统的信息获取和处理过程必须遵循奈奎斯特采样定理，即采样速率至少要大于原信号最高频率的2倍，这样才能从采样得到的离散信号中不失真地恢复出原始信号。然而随着信息技术的发展，这种以奈奎斯特采样定理为基础的信息获取和处理模式的弊端逐渐暴露出来。如对前端ADC的采样速率和处理速度要求高，采样数据量大，冗余度高等。

近年来出现的压缩感知(Compressed Sensing,CS)理论指出，在信号满足稀疏性的条件下，能够对信号进行全局观测，通过将压缩和采样合并进行的方式降低对信号的采样速率的要求，以远低于奈奎斯特采样率的速率进行采集，然后通过适当的重构算法恢复出原信号。CS理论是一个充分利用信号稀疏性或可压缩性的全新信号采集、编解码理论，为突破奈奎斯特采样定理的限制提供了一个新思路。CS理论的特殊思想对传统采样理论提出了巨大挑战，在信号处理领域带来了一场大的变革，引起广大信号处理专家以及数学家的重视，迅速成为信号处理领域的一个研究热点，在最近几年里压缩感知系统的设计吸引了大量国内外相关科研人员的广泛关注。

下面是对压缩感知基本原理的介绍。

对一个模拟信号进行奈奎斯特采样可以得到一列采样值，奈奎斯特采样定理指出，从这一列采样值可以无失真地恢复原模拟信号，所以这一列采样值组成的离散时间信号向量就是原模拟信号的等价表示。对于一个实值的有限长一维离散时间信号x，可以看作为一个空间N×1维(即N行1列)的列向量，它可以在一个基矩阵下进行展开，即：

x＝Ψα (1)

其中α是x被展开后的稀疏系数向量，Ψ是基矩阵。

x和α是同一个信号的等价表示，x是信号在时域的表示，α是信号在Ψ域的表示。如果α的非零个数比N小很多，则表明该信号是可压缩的。这里我们研究的是频域稀疏信号，即信号经过傅里叶变换后在频域是稀疏的，稀疏系数向量α就是信号的频谱，而且α中只有少量的若干值非零。对于频域稀疏信号而言，基矩阵Ψ就是离散傅里叶逆变换(InverseDiscrete Fourier Transform,IDFT)矩阵。

在信号具有频域稀疏性的前提下，可以用一个与IDFT矩阵Ψ不相关的观测矩阵Φ对信号x执行压缩观测，如公式(2)所示。观测矩阵Φ的维数为M×N(M＜＜N)，那么压缩观测后就可以得到M个观测值这些少量观测值中包含了重构信号x的足够信息。

y＝Φx (2)

将公式(1)代入公式(2)得到公式(3)，其中Θ＝ΦΨ，称为压缩矩阵。这样，对一个信号x进行压缩观测的过程就可以表示为图1的矩阵形式。

y＝ΦΨα＝Θα (3)

从信号的压缩观测中实现信号的恢复可以通过求解如下问题，获得一个唯一确定的解，即稀疏系数向量α，然后以将它与基矩阵Ψ相乘，就可以得到信号x＝Ψα。

argmin||α||₀subject to y＝Θα (4)

上述只是从理论方面进行的阐述，在实际应用中，随机解调是一种将压缩感知理论拓展到模拟域的技术和方法。其的结构如图2所示，将模拟信号x(t)与一个MLS序列p_c(t)相乘，即进行混频，然后经过一个模拟低通滤波器h(t)，最后，通过一个传统的低速ADC进行采样得到观测数据y(n)。

首先，在混频阶段，原信号x(t)与一个高速MLS序列p(t)在模拟乘法器(AD633)中相乘。我们实验采用的被测原信号是常见的多频点信号，其频谱由若干离散的频率分量组成，如图3、图4。我们所用的MLS序列是一种伪随机二值序列，其幅值近似随机地跳变为1或-1，跳变速率至少为信号最高频率的二倍。MLS序列的随机特性和高速的变化频率使其频谱看起来近似于带宽很宽的噪声，几乎覆盖整个频率轴，如图5、图6所示。这种频谱特点正好与原信号相反。众所周知，两个信号在时域相乘，在频域表现为频谱的卷积，频谱卷积后的带宽等于这两个信号的带宽之和。这样的宽带、类噪声频谱与原信号的稀疏频谱进行卷积就使得原信号的频谱被展开。形象地说，就好像将信号的频谱信息在整个频率轴上进行涂抹，这样频率轴上每一点都含有原信号的全局信息，并且由于MLS序列的参与，这一过程还对原信号进行了编码，这样，频率轴上每一点处的信息都具有了独特的标识。显然低频部分也包含信号的频谱信息，因此我们只需要滤除高频部分，采集低频信号这部分就能获得原信号的所有频谱信息并恢复原信号。

然后，在滤波阶段对混频之后的信号进行模拟低通滤波，滤掉高频部分，留下低频部分，使得信号的带宽变窄，这样，就可以用较低的采样率采集信号；接着，对滤波后的信号进行低速均匀采样，采样率需大于低通滤波器的截止频率的2倍；采样获得的数据值就构成观测值向量；最后利用MLS序列、系统的单位脉冲响应h(t)，IDFT矩阵构造压缩观测过程的压缩矩阵Θ，获得压缩矩阵Θ和采样数据以后，就可以信号重构算法恢复出原信号的波形和频谱了。

综上所述，混频和低通滤波的过程就是对信号进行压缩观测的过程，均匀采样环节是获取观测值的过程，信号重构环节则是利用压缩感知的一些重构算法实现对原信号的解调，恢复出原信号。通过随机解调可以实现对模拟信号的全局观测，大大降低信号的采样率。而压缩矩阵Θ是压缩感知模型中非常关键的一个部分，它决定了对原信号进行何种压缩观测，以及能否高概率重构出原信号。只有求得压缩观测过程的压缩矩阵Θ，才能利用一些优化算法从观测值y中恢复原信号。

频域稀疏信号，即信号经过离散傅里叶变换(Discrete Fourier Transform,DFT)后，幅值较大的频谱分量个数较少。这样的信号在生活中较常见，例如AM信号，雷达信号，音乐信号等等。当信号在傅里叶变换域中具有稀疏性时，上面压缩感知理论模型中的正交基矩阵Ψ就是离散傅里叶逆变换(IDFT)矩阵。压缩矩阵Θ的获取则是通过上面提到的公式Θ＝ΦΨ进行计算。IDFT矩阵Ψ的维数是N×N，观测矩阵Φ的维数是M×N，那么计算压缩矩阵Θ就需要M×N²次复数乘法运算，M×N×(N-1)次加法运算。其计算量为O(M×N²)，尤其是当N增大时，计算量，存储空间，运算时间将按几何级数急剧增加。例如，当信号x长度为1000，观测值y的点数为200时，则观测矩阵Φ的维数是200×1000，离散傅里叶变换矩阵Ψ的维数1000×1000,最后求得的压缩矩阵Θ维数为200×1000。取数据类型为double型，则需要的存储空间为5.6×10⁶Byte。在实际问题中，信号的长度N通常都较大，所以有必要开发快速算法，减小计算量。

压缩矩阵Θ是压缩感知模型中非常关键的一个部分，它是高概率重构原信号的条件之一。按照旧有的压缩矩阵计算公式(5)，需要计算观测矩阵Φ和IDFT矩阵Ψ的相乘，这种直接矩阵相乘的运算是将矩阵Φ的每一行元素与矩阵Ψ的每一列元素对应相乘再相加，当信号长度增加时，各个矩阵的维数也会随之增加，此时需要的计算量和内存会以几何级数的速度增加。在实际问题中，信号的长度通常都较大，这就意味着观测矩阵Φ和IDFT矩阵Ψ的都较大，从而暴露出原有的压缩矩阵存在计算量大，耗时，占用内存多等等一系列弊端。

Θ＝ΦΨ (5)

其中，Ψ为IDFT矩阵，如(6)所示。观测矩阵Φ，如(7)所示。

为了解决这些弊端，现利用DFT矩阵和IDFT矩阵的共轭关系、DFT矩阵的对称性和矩阵转置运算的规则，设计新的压缩矩阵Θ计算方法，然后引入FFT变换，从而实现了计算量的减小和速度的提升。

发明内容

本发明的目的是提出一种基于快速感知矩阵的频域稀疏信号重构方法，以解决现有的方法计算量大、耗费时间长，所需数据存储空间大、功耗大，实时性不高、效率不高的问题。

本发明为解决上述技术问题所采用的技术方案是：

本发明所述的一种基于快速感知矩阵的频域稀疏信号重构方法，包括以下步骤：步骤一、用模数转换器采集滤波器的输出信号，获得一系列采样值，记为y(m),(m＝1,2,...,M)，根据MLS序列本身的产生模式和输出采样率fs，计算出采样时间t内输入给乘法器的MLS序列的一系列值，记为p(n),(n＝1,2,...,N)，

步骤二、给乘法器的一个输入端输入1V直流信号，给另一个输入端输入矩形脉冲信号，信号的低电平为0V，高电平为1V，且高电平持续时间为0.1ms；与此同时，用模数转换器采集低通滤波器的输出信号，采样率和采样时间与步骤一中的fs和t相同，采集的结果为乘法器和低通滤波器的脉冲响应，记为h(n),(n＝1,2,...N),N＝fs×t；

步骤三、用脉冲响应h(n)构造矩阵H，用MLS序列p(n)构造矩阵P；将矩阵H和矩阵P相乘，计算观测矩阵Φ，方法如下：

步骤三一、定义一个N行N列的全零矩阵H0，用h(n)的第1个元素替换掉矩阵H0的第一行的第1个零元素；然后将h(n)前2个元素倒序排列，替换掉矩阵H0的第二行的前2个零元素；以此类推，将h(n)前i个元素倒序排列，替换掉矩阵H0的第i行的前i个零元素，如下面的公式所示：

步骤三二、定义一个N行N列的全零矩阵P0，用p(1)到p(n)代替P中的对角线元素，结果如下所示：

步骤三三、将矩阵H和P相乘得到观测矩阵Φ，即Φ＝HP；

步骤四、对N×N的观测矩阵Φ进行如下处理：每隔C行抽取一行，共抽取M行，组成新的M×N的观测矩阵Φ，C为整数；

步骤五、根据公式(1)得到傅里叶逆变换矩阵即IDFT矩阵Ψ：

步骤六、利用IDFT矩阵Ψ和DFT矩阵的共轭性质将IDFT矩阵转换为DFT矩阵，在式(2)的左右两边同时进行共轭操作，得到式(3)

Θ＝ΦΨ (2)

其中Θ是感知矩阵，

其中为离散傅里叶变换(DFT)矩阵，如式(4)所示，

因为本方法针对的是频域稀疏信号，是一种经过离散傅里叶变换后具有稀疏性的信号，所以所用的基矩阵Ψ是IDFT矩阵。又因为IDFT矩阵是DFT矩阵的共轭矩阵，FFT算法中使用的变换矩阵是DFT矩阵而不是IDFT矩阵，所以可以利用IDFT矩阵和DFT矩阵的共轭性质将IDFT矩阵转换为DFT矩阵，此步骤是后续其他步骤的基础。另外，针对频域稀疏信号，这里使用的观测矩阵Φ具有随机特性，共轭操作不会影响其随机特性。

步骤七、在式(3)两边同时进行转置操作，得到式(5)，

再根据矩阵的转置运算规则，将式(5)变为式(6)，

其中是的转置矩阵，是的转置矩阵，该转置操作的重要作用在于使得和的相对位置不同于和的相对位置，将调整到了之前。FFT算法要求DFT矩阵必须在另一个矩阵的前面，这样才满足FFT的使用条件。虽然该转置操作并没有将DFT矩阵提前到的前面，只是将DFT矩阵的转置矩阵提前到的前面，但这样的操作却为FFT的使用创造了一个必要条件。此步骤和步骤八具有紧密联系，步骤八中将DFT矩阵转置矩阵转变为DFT矩阵为FFT的使用创造了另外一个必要条件。这两个步骤结合起来才满足了FFT使用的全部要求，故该步骤具有重要作用。

步骤八、由式(4)可知，DFT矩阵是对称阵，故有式(7)所示关系，

将式(7)带入式(6)后可得式(8)，

此步骤的作用是利用DFT矩阵的对称性将DFT矩阵转置矩阵转变为DFT矩阵构成FFT使用的另一个必要条件，该条件和步骤二中创造的条件结合起来满足了FFT使用的全部要求。

步骤九、对的每一列都进行FFT操作，结果就是

步骤十、对步骤九得到的进行共轭转置得到压缩矩阵Θ，如式(9)所示，其中(·)^*表示共轭转置，此步骤的作用是将计算结果转换为最终需要的结果，对共轭操作和转置操作没有顺序要求，即孰先孰后没有影响。

步骤十一、将步骤一中获得的采样值y(m)和步骤十中获得的感知矩阵Θ作为参数，利用正交匹配追踪算法重构出原信号系数向量正交匹配追踪算法的步骤如下：

(1)初始化各参数，残差r₀＝y，信号支撑集支撑矩阵计数变量l＝1；

(2)解决如下优化问题，

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

其中θ_j(j＝1,…,N)是感知矩阵Θ的第j列元素组成的列向量，问题解决过程为：计算残差r_l-1和感知矩阵Θ每一列(θ_j)的内积，记录最大内积对应的列向量θ_j和θ_j在的Θ中的位置j，位置j即为本次找到的支撑集元素λ_l；

(3)将找到的支撑集元素λ_l添加到信号支撑集中，如(11)式所示，将θ_j添加到支撑矩阵中，如(12)式所示，

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

(4)更新残差，如(13)式所示，

其中是Θ_Λl的伪逆，

(5)计数变量l加1，若l≤K，则跳回到(2)，否则，执行(6)；

(6)输出重构信号的系数向量

步骤十二、利用步骤十一中重构出的系数向量和基矩阵Ψ得到重构信号

本发明的有益效果是：

一、IDFT矩阵Ψ的维数是N×N，观测矩阵Φ的维数是M×N，那么计算压缩矩阵Θ时如果采用原来的矩阵相乘法，则需要M×N²次复数乘法运算，M×N×(N-1)次加法运算。而采用本发明所述方法，则所需的计算量是次复数乘法，M×Nlog₂N次加法。计算量比较如表格1所示。

二、由于一次加法和一次复数乘法的运算量相比很小，所以图8，图9中只比较复数乘法次数，从中可以看出，在M一定时，矩阵相乘法的复数乘法次数随着N的增大以指数方式急剧增加，相比之下，本发明方法复数乘法次数则受N的影响很小，增大不明显。可见，和现有的方法相比本发明的计算量小、耗费时间短，所需数据存储空间小、功耗小，实时性更好、效率更高。

附图说明

图1压缩观测过程的矩阵表示；图2随机解调结构图；图3原信号部分波形；图4原信号频谱；图5 MLS序列部分波形；图6 MLS序列频谱；图7矩阵相乘法与本发明方法的复数乘法次数比较；图8矩阵相乘法与本发明方法的复数乘法次数比较；图9矩阵相乘法与本发明方法构造的压缩矩阵之间的复数误差大小比较。

具体实施方式

具体实施方式一：本实施方式所述的步骤一、用模数转换器采集滤波器的输出信号，获得一系列采样值，记为y(m),(m＝1,2,...,M)，根据MLS序列本身的产生模式和输出采样率fs，计算出采样时间t内输入给乘法器的MLS序列的一系列值，记为p(n),(n＝1,2,...,N)，

步骤二、给乘法器的一个输入端输入1V直流信号，给另一个输入端输入矩形脉冲信号，信号的低电平为0V，高电平为1V，且高电平持续时间为0.1ms；与此同时，用模数转换器采集低通滤波器的输出信号，采样率和采样时间与步骤一中的fs和t相同，采集的结果为乘法器和低通滤波器的脉冲响应，记为h(n),(n＝1,2,...N),N＝fs×t；

步骤三、用脉冲响应h(n)构造矩阵H，用MLS序列p(n)构造矩阵P；将矩阵H和矩阵P相乘，计算观测矩阵Φ，方法如下：

步骤三一、定义一个N行N列的全零矩阵H0，用h(n)的第1个元素替换掉矩阵H0的第一行的第1个零元素；然后将h(n)前2个元素倒序排列，替换掉矩阵H0的第二行的前2个零元素；以此类推，将h(n)前i个元素倒序排列，替换掉矩阵H0的第i行的前i个零元素，如下面的公式所示：

步骤三二、定义一个N行N列的全零矩阵P0，用p(1)到p(n)代替P中的对角线元素，结果如下所示：

步骤三三、将矩阵H和P相乘得到观测矩阵Φ，即Φ＝HP；

步骤四、对N×N的观测矩阵Φ进行如下处理：每隔C行抽取一行，共抽取M行，组成新的M×N的观测矩阵Φ，C为整数；

步骤五、根据公式(1)得到傅里叶逆变换矩阵即IDFT矩阵Ψ：

步骤六、利用IDFT矩阵Ψ和DFT矩阵的共轭性质将IDFT矩阵转换为DFT矩阵，在式(2)的左右两边同时进行共轭操作，得到式(3)

Θ＝ΦΨ (2)

其中Θ是感知矩阵，

其中为离散傅里叶变换(DFT)矩阵，如式(4)所示，

因为本方法针对的是频域稀疏信号，是一种经过离散傅里叶变换后具有稀疏性的信号，所以所用的基矩阵Ψ是IDFT矩阵。又因为IDFT矩阵是DFT矩阵的共轭矩阵，FFT算法中使用的变换矩阵是DFT矩阵而不是IDFT矩阵，所以可以利用IDFT矩阵和DFT矩阵的共轭性质将IDFT矩阵转换为DFT矩阵，此步骤是后续其他步骤的基础。另外，针对频域稀疏信号，这里使用的观测矩阵Φ具有随机特性，共轭操作不会影响其随机特性。

步骤七、在式(3)两边同时进行转置操作，得到式(5)，

再根据矩阵的转置运算规则，将式(5)变为式(6)，

其中是的转置矩阵，是的转置矩阵，该转置操作的重要作用在于使得和的相对位置不同于和的相对位置，将调整到了之前。FFT算法要求DFT矩阵必须在另一个矩阵的前面，这样才满足FFT的使用条件。虽然该转置操作并没有将DFT矩阵提前到的前面，只是将DFT矩阵的转置矩阵提前到的前面，但这样的操作却为FFT的使用创造了一个必要条件。此步骤和步骤三具有紧密联系，步骤三中将DFT矩阵转置矩阵转变为DFT矩阵为FFT的使用创造了另外一个必要条件。这两个步骤结合起来才满足了FFT使用的全部要求，故该步骤具有重要作用。

步骤八、由式(4)可知，DFT矩阵是对称阵，故有式(7)所示关系，

将式(7)带入式(6)后可得式(8)，

此步骤的作用是利用DFT矩阵的对称性将DFT矩阵转置矩阵转变为DFT矩阵构成FFT使用的另一个必要条件，该条件和步骤二中创造的条件结合起来满足了FFT使用的全部要求。

步骤九、对的每一列都进行FFT操作，结果就是

步骤十、对步骤九得到的进行共轭转置得到压缩矩阵Θ，如式(9)所示，其中(·)^*表示共轭转置，此步骤的作用是将计算结果转换为最终需要的结果，对共轭操作和转置操作没有顺序要求，即孰先孰后没有影响。

步骤十一、将步骤一中获得的采样值y(m)和步骤十中获得的感知矩阵Θ作为参数，利用正交匹配追踪算法重构出原信号系数向量正交匹配追踪算法的步骤如下：

(1)初始化各参数，残差r₀＝y，信号支撑集支撑矩阵计数变量l＝1；

(2)解决如下优化问题，

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

其中θ_j(j＝1,…,N)是感知矩阵Θ的第j列元素组成的列向量，问题解决过程为：计算残差r_l-1和感知矩阵Θ每一列(θ_j)的内积，记录最大内积对应的列向量θ_j和θ_j在的Θ中的位置j，位置j即为本次找到的支撑集元素λ_l；

(3)将找到的支撑集元素λ_l添加到信号支撑集中，如(11)式所示，将θ_j添加到支撑矩阵中，如(12)式所示，

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

(4)更新残差，如(13)式所示，

其中是的伪逆，

(5)计数变量l加1，若l≤K，则跳回到(2)，否则，执行(6)；

(6)输出重构信号的系数向量

步骤十二、利用步骤十一中重构出的系数向量和基矩阵Ψ得到重构信号

结合图1解释本实施方式。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一中的y(m)的获取过程为：

被测信号输入乘法器的一个输入端，乘法器的另一个输入端用于输入MLS序列；乘法器的输出端与低通滤波器的输入端相连，低通滤波器的输出端与模数转换器的输入端相连，模数转换器的输出端与上位机的输入端相连；

将一路被测模拟信号x(t)与一路MLS序列p(t)输入到乘法器中进行相乘，相乘后的信号经过低通滤波器滤波，滤波后的信号被模数转换器采样，采样时所用的采样率为fc，采样时间为t，得到一系列采样值，记为y(m),m＝1,2,...,M，M＝fc×t。结合图2解释本实施方式。其它步骤与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤一中的p(n)的获取过程为：

MLS序列是根据设定好的模式计算出来的一系列数值，然后借助任意波形发生器将这些数值按设定的采样率fs逐个输出，形成一路MLS序列，由于其产生模式和输出采样率都已知，故采样时间t内MLS序列的所有值都能计算出来，记为p(n),(n＝1,2,...,N)，N＝fs×t；

MLS序列的输出采样率fs和采样滤波器输出信号所用的采样率fc有如下关系：fs＝C×fc，同时自然有如下关系N＝C×M，C是一个整数。结合图3至图6解释本实施方式。其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤九的操作过程为：

步骤九一、定义i为循环次数变量，给i赋初值1，即i＝1；

步骤九二、判断i是否大于M，如果小于等于M，顺序执行步骤九三和步骤九四，否则执行步骤九五，其中M为观测矩阵Φ的行数；

步骤九三、其中和的下标表示该元素对应和中的第几列第几行；

步骤九四、将i的值加1，即i＝i+1，然后跳到步骤九二；

步骤九五、退出此循环。其它步骤及参数与具体实施方式一至三之一相同。

本发明的实验验证如下：

设本发明方法的计算量和矩阵相乘法的计算量的比值为R，则当N＝1000时，本发明方法的计算量仅为矩阵相乘法的计算量的0.5％。

在计算时间、占用内存方面的比较如表格2、表格3所示，从中可以看出本发明的方法在这两方面的表现都明显优于原来的矩阵相乘法。

另外在计算精度方面也进行了相应的实验，将矩阵相乘法构造的压缩矩阵减去本发明方法构造的压缩矩阵以计算它们之间的偏差矩阵。由于压缩矩阵的元素都是复数，所以求得的偏差矩阵的元素也都是复数，如图9所示，横坐标表示偏差的实部大小，纵坐标表示偏差的虚部大小。进行多次实验，计算偏差矩阵每个元素的均值。从图中可以看出偏差矩阵的所有元素其均值都极小，在2×10^-17以内，说明本发明方法构造的压缩矩阵与原来的矩阵相乘法精度相当。

表格1比较本发明的方法和矩阵相乘法的计算量

表格2比较本发明的方法和矩阵相乘法的计算时间

表格3比较本发明的方法和矩阵相乘法的占用内存

Claims (4)

1.面向频域稀疏信号的信号重构方法，其特征在于所述方法包括以下步骤：

步骤一、用模数转换器采集滤波器的输出信号，获得一系列采样值，记为y(m)，m＝1,2,...,M，根据MLS序列本身的产生模式和输出采样率fs，计算出采样时间t内输入给乘法器的MLS序列的一系列值，记为p(n)，n＝1,2,...,N，

步骤二、给乘法器的一个输入端输入1V直流信号，给另一个输入端输入矩形脉冲信号，且高电平持续时间为0.1ms；与此同时，用模数转换器采集低通滤波器的输出信号，采样率和采样时间与步骤一中的fs和t相同，采集的结果为乘法器和低通滤波器的脉冲响应，记为h(n)，n＝1,2,...N，N＝fs×t；

步骤三、用脉冲响应h(n)构造矩阵H，用MLS序列p(n)构造矩阵P；将矩阵H和矩阵P相乘，计算观测矩阵Φ，方法如下：

步骤三一、定义一个N行N列的全零矩阵H0，用h(n)的第1个元素替换掉矩阵H0的第一行的第1个零元素；然后将h(n)前2个元素倒序排列，替换掉矩阵H0的第二行的前2个零元素；以此类推，将h(n)前i个元素倒序排列，替换掉矩阵H0的第i行的前i个零元素，如下面的公式所示：

步骤三二、定义一个N行N列的全零矩阵P0，用p(1)到p(n)代替P中的对角线元素，结果如下所示：

步骤三三、将矩阵H和P相乘得到观测矩阵Φ，即Φ＝HP；

步骤四、对N×N的观测矩阵Φ进行如下处理：每隔C行抽取一行，共抽取M行，组成新的M×N的观测矩阵Φ，C为整数；

步骤五、根据公式(1)得到离散傅里叶逆变换矩阵即IDFT矩阵Ψ：

$\mathrm{\&Psi;}=\left[\begin{array}{cccccc}1& 1& 1& 1& ...& 1\\ 1& \stackrel{\&OverBar;}{w}& {\stackrel{\&OverBar;}{w}}^2& {\stackrel{\&OverBar;}{w}}^3& ...& {\stackrel{\&OverBar;}{w}}^{N-1}\\ 1& {\stackrel{\&OverBar;}{w}}^2& {\stackrel{\&OverBar;}{w}}^4& {\stackrel{\&OverBar;}{w}}^6& ...& {\stackrel{\&OverBar;}{w}}^{2\left(N-1\right)}\\ 1& {\stackrel{\&OverBar;}{w}}^3& {\stackrel{\&OverBar;}{w}}^6& {\stackrel{\&OverBar;}{w}}^9& ...& {\stackrel{\&OverBar;}{w}}^{3\left(N-1\right)}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 1& {\stackrel{\&OverBar;}{w}}^{N-1}& {\stackrel{\&OverBar;}{w}}^{2\left(N-1\right)}& {\stackrel{\&OverBar;}{w}}^{3\left(N-1\right)}& ...& {\stackrel{\&OverBar;}{w}}^{\left(N-1\right)\left(N-1\right)}\end{array}\right],\stackrel{\&OverBar;}{w}=e^{j\frac{2\mathrm{\&pi;}}{N}}---\left(1\right)$

步骤六、利用IDFT矩阵Ψ和DFT矩阵的共轭性质将IDFT矩阵转换为DFT矩阵，在式(2)的左右两边同时进行共轭操作，得到式(3)

Θ＝ΦΨ (2)

其中Θ是感知矩阵，

$\stackrel{\&OverBar;}{\mathrm{\&Theta;}}=\stackrel{\&OverBar;}{\mathrm{\&Phi;}}\stackrel{\&OverBar;}{\mathrm{\&Psi;}}---\left(3\right)$

其中为离散傅里叶变换矩阵，如式(4)所示，

$\stackrel{\&OverBar;}{\mathrm{\&Psi;}}=\left[\begin{array}{cccccc}1& 1& 1& 1& ...& 1\\ 1& w& w^2& w^3& ...& w^{N-1}\\ 1& w^2& w^4& w^6& ...& w^{2\left(N-1\right)}\\ 1& w^3& w^6& w^9& ...& w^{3\left(N-1\right)}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 1& w^{N-1}& w^{2\left(N-1\right)}& w^{3\left(N-1\right)}& ...& w^{\left(N-1\right)\left(N-1\right)}\end{array}\right],w=e^{-j\frac{2\mathrm{\&pi;}}{N}}---\left(4\right)$

步骤七、在式(3)两边同时进行转置操作，得到式(5)，

${\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}^T={\left(\stackrel{\&OverBar;}{\mathrm{\&Phi;}}\stackrel{\&OverBar;}{\mathrm{\&Psi;}}\right)}^T---\left(5\right)$

再根据矩阵的转置运算规则，将式(5)变为式(6)，

${\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}^T={\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}^T{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}^T---\left(6\right)$

其中是的转置矩阵，是的转置矩阵；

步骤八、由式(4)可知，DFT矩阵是对称阵，故有式(7)所示关系，

${\stackrel{\&OverBar;}{\mathrm{\&Psi;}}}^T=\stackrel{\&OverBar;}{\mathrm{\&Psi;}}---\left(7\right)$

将式(7)带入式(6)后可得式(8)，

${\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}^T=\stackrel{\&OverBar;}{\mathrm{\&Psi;}}{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}^T---\left(8\right);$

步骤九、对的每一列都进行FFT操作，结果就是

步骤十、对步骤九得到的进行共轭转置得到压缩矩阵Θ，如式(9)所示，其中(·)^*表示共轭转置，

$\mathrm{\&Theta;}={\left({\stackrel{\&OverBar;}{\mathrm{\&Theta;}}}^T\right)}^{*}=\left[\begin{array}{cccccc}{\mathrm{\&theta;}}_{1,1}& {\mathrm{\&theta;}}_{1,2}& {\mathrm{\&theta;}}_{1,3}& {\mathrm{\&theta;}}_{1,4}& ...& {\mathrm{\&theta;}}_{1,N}\\ {\mathrm{\&theta;}}_{2,1}& {\mathrm{\&theta;}}_{2,2}& {\mathrm{\&theta;}}_{2,3}& {\mathrm{\&theta;}}_{2,4}& ...& {\mathrm{\&theta;}}_{2,N}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\&theta;}}_{M,1}& {\mathrm{\&theta;}}_{M,2}& {\mathrm{\&theta;}}_{M,3}& {\mathrm{\&theta;}}_{M,4}& ...& {\mathrm{\&theta;}}_{M,N}\end{array}\right]={\left(\left[\begin{array}{cccc}{\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{1,1}& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{2,1}& ...& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{M,1}\\ {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{1,2}& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{2,2}& ...& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{M,2}\\ {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{1,3}& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{2,3}& ...& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{M,3}\\ {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{1,4}& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{2,4}& ...& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{M,4}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{1,N}& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{2,N}& ...& {\stackrel{\&OverBar;}{\mathrm{\&theta;}}}_{M,N}\end{array}\right]\right)}^{*}---\left(9\right);$

步骤十一、将步骤一中获得的采样值y(m)和步骤十中获得的感知矩阵Θ作为参数，利用正交匹配追踪算法重构出原信号系数向量正交匹配追踪算法的步骤如下：

(1)初始化各参数，残差r₀＝y，信号支撑集支撑矩阵计数变量l＝1；

(2)解决如下优化问题，

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

其中θ_j是感知矩阵Θ的第j列元素组成的列向量，其中j＝1,…,N，问题解决过程为：计算残差r_l-1和感知矩阵Θ每一列(θ_j)的内积，记录最大内积对应的列向量θ_j和θ_j在的Θ中的位置j，位置j即为本次找到的支撑集元素λ_l；

(3)将找到的支撑集元素λ_l添加到信号支撑集中，如(11)式所示，将θ_j添加到支撑矩阵中，如(12)式所示，

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

${\mathrm{\&Theta;}}_{{\mathrm{\&Lambda;}}_l}={\mathrm{\&Theta;}}_{{\mathrm{\&Lambda;}}_{l-1}}\&cup;\left\{{\mathrm{\&theta;}}_l\right\}---\left(12\right)$

(4)更新残差，如(13)式所示，

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

其中是的伪逆，

(5)计数变量l加1，若l≤K，则跳回到(2)，否则，执行(6)；

(6)输出重构信号的系数向量

$\widehat{\mathrm{\&alpha;}}={\mathrm{\&Theta;}}_{{\mathrm{\&Lambda;}}_l}^{+}y,and{\widehat{\mathrm{\&alpha;}}}_{\{1,...,N\}-{\mathrm{\&Lambda;}}_l}=0---\left(14\right)$

步骤十二、利用步骤十一中重构出的系数向量和IDFT矩阵Ψ得到重构信号

$\hat{x}=\mathrm{\&Psi;}\widehat{\mathrm{\&alpha;}}---\left(15\right).$

2.如权利要求1所述的面向频域稀疏信号的信号重构方法，其特征在于步骤一中的y(m)的获取过程为：

将一路被测模拟信号x(t)与一路MLS序列p(t)输入到乘法器中进行相乘，相乘后的信号经过低通滤波器滤波，滤波后的信号被模数转换器采样，采样时所用的采样率为fc，采样时间为t，得到一系列采样值，记为y(m),m＝1,2,...,M，M＝fc×t。

3.如权利要求2所述的面向频域稀疏信号的信号重构方法，其特征在于步骤一中的p(n)的获取过程为：

MLS序列是根据设定好的模式计算出来的一系列数值，然后借助任意波形发生器将这些数值按设定的采样率fs逐个输出，形成一路MLS序列，故采样时间t内MLS序列的所有值都能计算出来，记为p(n)，n＝1,2,...,N，N＝fs×t；

MLS序列的输出采样率fs和采样滤波器输出信号所用的采样率fc有如下关系：fs＝C×fc，同时自然有如下关系N＝C×M，C是一个整数。

4.如权利要求3所述的面向频域稀疏信号的信号重构方法，其特征在于步骤九的操作过程为：

步骤九一、定义i为循环次数变量，给i赋初值1，即i＝1；

步骤九二、判断i是否大于M，如果小于等于M，顺序执行步骤九三和步骤九四，否则执行步骤九五，其中M为观测矩阵Φ的行数；

步骤九三、其中和的下标表示该元素对应和中的第几列第几行；

步骤九四、将i的值加1，即i＝i+1，然后跳到步骤九二；

步骤九五、退出此循环。