CN110896308B - 一种单音信号重构方法 - Google Patents

Abstract

本发明公开了一种单音信号重构方法，对于原始单音信号x的测量向量y和测量矩阵，按照频率搜索步长，在频率f的取值范围内选取频率，根据选取的频率计算均方误差，将计算得到的均方误差与迭代过程中前一步的误差比较，直到迭代匹配出均方误差最小值，再将均方误差最小值所对应的选取的频率作为单音信号频率的估计值，最后，根据所述单音信号频率的估计值构造出原始单音信号x的估计值；每一次迭代过程中的频率搜索步长小于前一次迭代过程中的频率搜索步长，同时每一次迭代过程中的频率大于前一次迭代过程中的频率。本发明从单音信号特性出发，通过迭代逐步缩小稀疏表示矩阵的频率范围，提高与信号的匹配程度，从而解决了稀疏表示矩阵的失配问题。

Description

一种单音信号重构方法

技术领域

本发明涉及从压缩感知测量数据中重构原始数据的方法，具体涉及从单音信号的压缩感知测量数据中重构原始数据的方法。

背景技术

压缩感知技术利用信号的稀疏特性进行处理得到测量数据，在后续处理时，再利用重构算法从测量数据中恢复出原始数据。测量数据与原始数据相比数据量大幅度降低，可以有效减轻数据存储和传输压力。可见，压缩感知中的数据测量和重构实际上是数据的压缩和解压缩过程，属于信源编解码范畴。

压缩感知的测量过程可表示为，

y＝Φx (1)

其中x为待测量的原始数据，它是一个N维列向量，即x＝[x₁,x₂,...,x_N]^T；Φ为测量矩阵，其维数为M×N，并且M＜N；y是得到的测量数据，它是一个M维列列向量。可见，测量过程将一个N维列向量转化为一个M维列列向量，而M＜N，从而实现数据压缩。

鉴于应用压缩感知的前提条件是信号具有稀疏特性，因此假设x用(2)式表示，

x＝Ψα (2)

其中ψ为N×N表示基矩阵，α＝[α₁,α₂,…,α_N]^T为信号x在Ψ上的线性投影系数。如果系数向量α中非零元素的数量为K，且K＜＜N，则称x在基Ψ上是稀疏的，K为信号的稀疏度，并称Ψ为x的稀疏表示矩阵。

将(2)式代入(1)式得到，

y＝ΦΨα＝Θα (3)

其中，Θ称为恢复矩阵，是一个M×N矩阵。在Θ满足一定条件下，α可求解下述最优化问题精确重构，

||α||₀表示α中非零元素的个数。s.t.是“subject to(满足)”的缩写，即(4)表示：在满足y＝Θα的所有α中，非零元素个数最少的即为原始系数矩阵的估计。

求解上述最优化问题是非常困难的，所以一般将上述最优化问题转换为下述可求解的最优化问题

其中

为了能够从有限的测量值中恢复原始信号，恢复矩阵Θ必须满足约束等距性条件(Restricted Isometry Property，RIP)，即

且M满足

M≤cKlog(N/K) (7)

其中δ_K∈[0,1)为K阶约束等距常数；c是一个与重构精度有关的常数；||·||₂表示向量元素平和的开方。

重构出系数向量之后，即可以根据(2)式恢复出原始数据，即

在压缩感知的应用对象中，单音信号(正弦信号或余弦信号)是一种常见的频域稀疏信号，其稀疏度为2。频域稀疏信号的稀疏表示矩阵一般采用付利叶矩阵，在归一化采样频率下可表示为

Ψ＝[ψ₁ ψ₂ … ψ_N] (9)

其中，

可知，付利叶矩阵是由N个离散频率原子向量组成的，但是单音信号的频率是未知的而且是任意的，有可能不与付利叶矩阵中任何原子向量对应的频率相等，此时如果还利用此矩阵作为稀疏表示矩阵，则会导致稀疏度扩散，即稀疏度将远大于2，这就是所谓的稀疏表示矩阵失配现象。

稀疏表示矩阵失配导致单音信号的稀疏度增大，如果此时测量时仍按照稀疏度为2设计测量参数M，则无法从测量数据中精确恢复出原始数据，而如果增大M以适应稀疏度的增大，则会大幅度降低压缩感知的压缩效率，甚至会失去应用压缩感知技术的必要性。目前，处理测量矩阵失配问题，主要是采用冗余稀疏表示矩阵的方法，即增大稀疏表示矩阵原子向量的个数，以减小频率间隔，增大频率分辨率，但是由于频率是任意的，即使采用冗余稀疏表示矩阵也可能出现失配问题，因此无法完全解决这个问题。

发明内容

本发明旨在提出一种单音信号重构方法，该方法从单音信号特性出发，通过逐步缩小稀疏表示矩阵的频率范围，提高分辨率，动态构造与单音信号匹配的稀疏表示矩阵，从而实现与信号的匹配，解决单音信号重构中的稀疏表示矩阵失配问题，而且前端可用同一测量参数进行测量，而无需考虑信号频率问题。

本发明的技术方案如下：

一种单音信号重构方法，其特征在于：对于原始单音信号x的测量向量y和测量矩阵Φ，按照频率搜索步长Δf，在频率f的取值范围内选取频率，根据选取的频率计算均方误差，将计算得到的均方误差与迭代过程中前一步的误差比较，直到迭代匹配出均方误差最小值，再将均方误差最小值所对应的选取的频率作为单音信号频率的估计值，最后，根据所述单音信号频率的估计值构造出原始单音信号x的估计值每一次迭代过程中的频率搜索步长小于前一次迭代过程中的频率搜索步长，同时每一次迭代过程中选取的频率大于前一次迭代过程中选取的频率。因为，步长过小会使得所需运算次数太大，运算量增加，而步长过大又会造成精度很难满足要求；为此，本发明采取步长由大到小，逐步缩小搜索范围的迭代匹配重构方法，在计算量增加有限的情况下，提高估计精度。

所述频率f满足：f∈[f_a,f_b]，f_a表示频率搜索范围的低端频率，f_b表示频率搜索范围的高端频率。

所述频率搜索步长Δf满足：Δf＝(f_b-f_a)/L，其中，L为搜索因数，表示将频率搜索范围f_w，划分为L份，L为不小于2的整数。

所述频率f，进一步具体的取值范围可以是f∈[0,0.5]，即f_a＝0，f_b＝0.5。

所述单音信号重构方法的迭代匹配步骤如下，

(1)初始化：Δf＝0.5/L，f_w＝0.5，f_a＝0，f_b＝0.5；

(2)k为记数器，E_old表示迭代过程中前一步得到的误差，令k的初始值为0，E_old的初始值为inf，inf代表计算机所能表示的最大数；

(3)基于Δf、f_a、f_b，得到频率：f＝f_a+kΔf；

(4)构造本次迭代过程中的稀疏表示矩阵：Ψ_f＝[ψ_fψ_-f]；

(5)计算本次迭代过程中原始系数矩阵的估计值本次迭代过程中实单音信号的估计值/>本次迭代过程中测量向量的估计值/>

(6)根据步骤(5)中的三个估计值计算本次迭代过程中的均方误差E(f)，并将计算得到的均方误差E(f)和E_old比较；

(7)如果E(f)＜E_old，则令E_old＝E(f)，f_e＝f，否则转向(8)；

(8)根据步骤(7)，如果f＝f_b，则计算f_w＝f_w/2，f_a＝max(0,f_e-f_w/2)，f_b＝min(0.5,f_e+f_w/2)，Δf＝(f_b-f_a)/L，并再进一步判断Δf是否为0：

如果Δf＝0，则令最后构造得到/>其中Re{·}表示取实部；

如果Δf≠0，返回步骤(2)；

(9)根据步骤(7)，如果f≠f_b，则令k＝k+11，又转到步骤(3)重新取频率再次迭代计算，直到最终迭代计算出最小的E(f)，才计算f_w＝f_w/2，f_a＝max(0,f_e-f_w/2)，f_b＝min(0.5,f_e+f_w/2)，Δf＝(f_b-f_a)/L；再判断Δf是否为0，当Δf＝0，才令构造出/>

上述步骤中，计算原始系数矩阵的估计值的公式为：/>

上述步骤中，计算实单音信号的估计值的公式为：/>

上述步骤中，计算测量向量的估计值的公式为：/>

上述步骤中，计算均方误差E(f)的公式为：

所述原始单音信号x为实信号且不存在噪声，频率为f_c，采样率为归一化采样率，即采样率为1Hz。

本发明方法实现的原理为：

假设原始单音信号x为实信号且不存在噪声，频率为f_c，采样率为归一化采样率，即采样率为1Hz。对于任意频率f(f∈[0,0.5])可以构造实单音信号x_f，

x_f＝Ψ_fα_f (11)

Ψ_f是一个N×2维矩阵，可表示为

Ψ_f＝[ψ_f ψ_-f] (12)

其中，

ψ_f＝[1 e^j2πf/N ... e^j2πf(N-1)/N]^T (13)

ψ_-f＝[1 e^-j2πf/N ... e^{-j2πf(N-1)/N}]^T (14)

α_f为构成x_f的系数向量，可表示为

α_f＝[α₁(f) α₂(f)]^T (15)

显然，

x_f＝Ψ_fα_f＝α₁(f)ψ_f+α₁(f)ψ_-f (16)

利用同一测量矩阵对x_f进行测量，得到测量向量y_f，

y_f＝Φx_f＝ΦΨ_fα_f (17)

将y_f替换为y，并利用最小二乘法可求得α_f的估计值，

其中[·]^H表示对复矩阵进行转置，[·]^-1表示矩阵求逆。进而可得到x_f的估计值y_f的估计值/>

显然，当f＝f_c时，与y之间误差最小。如果利用均方误差来衡量/>与y之间的差异，定义为，

则可得到，

其中表示求取使E(f)最小的f值。

进而得到时的α_f的估计值/>为

为/>时的Ψ_f。最后得到x的估计值/>

其中Re{·}表示取实部。

式(22)是一个一维最优化问题，本发明采用频率搜索的方法，即在频率取值范围内，按照一定步长由小到大选取频率，按照(21)式计算误差，并将误差最小值所对应的频率作为单音信号频率的估计值。但是，在归一化采样频率下，f在0～0.5之间取值，如果步长过小的所需运算的次数太大，运算量增加，如果步长过大，精度又很难满足要求；为此，本发明采取步长由大到小，逐步缩小搜索范围的迭代匹配重构方法，在计算量增加有限的情况下，提高估计精度。

本发明的有益效果如下：

本发明从单音信号特性出发，逐步缩小步长逐步放大搜索频率，通过逐步缩小稀疏表示矩阵的频率范围，提高分辨率，以此动态构造与单音信号匹配的稀疏表示矩阵，实现与信号的匹配，解决单音信号重构中的稀疏表示矩阵失配问题，而且前端可用同一测量参数进行测量，而无需考虑信号频率问题。

附图说明

图1为实施例中重构误差随54次迭代次数的变化示意图。

图2为实施例中根据图1中的迭代过程，取前20个原始数据与重构数据对比示意图。

具体实施方式

本发明提供的一种单音信号重构方法，对于原始单音信号x的测量向量y和测量矩阵Φ，按照频率搜索步长Δf，在频率f的取值范围内选取频率，根据选取的频率计算均方误差，将计算得到的均方误差与迭代过程中前一步的误差比较，直到迭代匹配出均方误差最小值，再将均方误差最小值所对应的选取的频率作为单音信号频率的估计值，最后，根据所述单音信号频率的估计值构造出原始单音信号x的估计值每一次迭代过程中的频率搜索步长小于前一次迭代过程中的频率搜索步长，同时每一次迭代过程中选取的频率大于前一次迭代过程中选取的频率。因为，步长过小会使得所需运算次数太大，运算量增加，而步长过大又会造成精度很难满足要求；为此，本发明采取步长由大到小，逐步缩小搜索范围的迭代匹配重构方法，在计算量增加有限的情况下，提高估计精度。

具体的迭代匹配重构步骤如下，

输入：测量向量y，测量矩阵Φ，数据长度N，搜索因数L；

输出：重构数据

(1)初始化,Δf＝0.5/L，f_w＝0.5，f_a＝0，f_b＝0.5，Δf表示频率搜索步长，f_w表示频率搜索范围的大小，f_a表示频率搜索范围的低端频率，f_b表示频率搜索范围的高端频率；

(2)令k＝0,E_old＝inf，k为记数器，E_old表示迭代过程中前一步得到的误差，inf代表计算机所能表示的最大数；

(3)计算f＝f_a+kΔf；

(4)构造Ψ_f＝[ψ_f ψ_-f]；

(5)计算

(6)计算及/>

(7)计算均方误差

(8)如果E(f)＜E_old，则E_old＝E(f)，f_e＝f，否则转向(9)；

(9)如果f＝f_b，则转到(11)，否则转到(10)；

(10)k＝k+1，转到(3)；

(11)计算f_w＝f_w/2，f_a＝max(0,f_e-f_w/2)，f_b＝min(0.5,f_e+f_w/2),Δf＝(f_b-f_a)/L；

(12)如果Δf＝0转向(13)，否则转向(2)；

(13)

搜索因数L表式将频率搜索范围划分为L份，显然L为不小于2的整数。

基于上述基本方法，举出以下应用实例：

利用频率为1/7的单音信号来测试上述算法的性能，测试过程中L＝N＝1024，M＝16，测量矩阵采用高斯测量矩阵。图1给出了重构误差随迭代次数的变化，重构误差定义为

从图1中可以看出经过54次迭代后搜索步长达到0，迭代结束；此时的重构误差为-263dB，估计频率为0.142857142857143。图2给出了原始信号和重构信号的对比，为了便于比较，图2中只取了前20个数据。

从上述应用实例可以看出，本发明从单音信号特性出发，通过迭代逐步缩小稀疏表示矩阵的频率范围，提高与信号的匹配程度，从而解决了稀疏表示矩阵的失配问题。

Claims (6)

1.一种单音信号重构方法，其特征在于：对于频率为f_c的原始单音信号x的测量向量y和测量矩阵Φ，按照频率搜索步长Δf，在取值范围内选取频率f，构造频率为f的单音信号x_f，并利用测量矩阵Φ对x_f进行测量得到测量向量y_f，利用y替换y_f，再根据最小二乘法，得到x_f系数的估计值进而得到x_f的估计值/>在利用测量矩阵Φ得到y_f的估计值/>计算/>与y之间的均方误差，将均方误差最小值所对应频率作为单音信号频率的估计值f_e；然后，根据f_e确定新的频率范围及搜索步长，重复迭代匹配出均方误差最小值的过程，直到满足搜索条件为止；最后，根据所述单音信号频率的估计值f_e构造出原始单音信号x的估计值/>

所述频率f满足：f∈[f_a,f_b]，f_a表示频率搜索范围的低端频率，f_b表示频率搜索范围的高端频率；所述频率f的取值范围是f∈[0,0.5]，即f_a＝0，f_b＝0.5；

所述频率搜索步长Δf满足：Δf＝(f_b-f_a)/L，其中，L为搜索因数，表示将频率搜索范围f_w，划分为L份，L为大于或等于2的整数；

所述单音信号重构方法的迭代匹配步骤如下，

(1)初始化：Δf＝0.5/L，f_w＝0.5，f_a＝0，f_b＝0.5，记数器k＝0，误差的最小值E_old＝inf，inf代表计算机所能表示的最大数；

(2)基于Δf、f_a、f_b，得到频率：f＝f_a+kΔf；

(3)构造矩阵：Ψ_f＝[ψ_fψ_-f]；

(4)计算系数向量的估计值实单音信号的估计值/>本次迭代过程中测量向量的估计值/>

(5)计算与y之间的均方误差E(f)，并将计算得到的均方误差E(f)和E_old比较；

(6)根据步骤(5)的比较，如果E(f)＜E_old，则令E_old＝E(f)，f_e＝f，否则转到(7)；

(7)如果f＝f_b，转到(9)，否则转到(8)；

(8)令k＝k+1，转到(2)；

(9)计算f_w＝f_w/2，f_a＝max(0,f_e-f_w/2)，f_b＝min(0.5,f_e+f_w/2),

Δf＝(f_b-f_a)/L；

(10)进一步判断Δf是否为0：如果Δf＝0，则取f_c的估计值对应的系数矩阵向量估计值/>最后构造得到x的估计值/>其中Re{·}表示取实部，如果Δf≠0，返回步骤(2)。

2.根据权利要求1所述的一种单音信号重构方法，其特征在于：上述步骤中，计算原始系数矩阵的估计值的公式为：/>

3.根据权利要求1所述的一种单音信号重构方法，其特征在于：上述步骤中，计算实单音信号的估计值的公式为：/>

4.根据权利要求1所述的一种单音信号重构方法，其特征在于：上述步骤中，计算测量向量的估计值的公式为：/>

5.根据权利要求1所述的一种单音信号重构方法，其特征在于：上述步骤中，计算均方误差E(f)的公式为：

6.根据权利要求1所述的一种单音信号重构方法，其特征在于：所述原始单音信号x为实信号且不存在噪声，采样率为归一化采样率，即采样率为1Hz。