加性噪声中未知分量数的谐波恢复方法研究

Research on Harmonic Retrieval Method in Additive Noise without Components Number

摘要: 谐波恢复问题是信号处理中的经典问题，对于被加性噪声污染的谐波信号，现有的多重信号分类(MUSIC)方法可以实现良好的恢复效果，但在谐波分量数未知的情况下无法发挥其超分辨率性能。本文结合实值MUSIC算法的超分辨率和逐步LSE算法的逐步估计的特点，提出了一种逐步MUSIC算法，该算法在分量数已知的情况下具有良好的估计性能。在分量数未知的情况下，我们利用该算法估计出一个模型集，并采用BIC准则估计谐波分量数。最后通过仿真实验对所提方法的性能进行了验证。

Abstract: The retrieval of harmonic signal is a classic problem in signal processing. Consider harmonic signals in additive noise, among the proposed schemes in related literatures, the MUSIC (Multiple Signal Classification) algorithm can achieve good retrieval performance. However, it can’t play its super-resolution performance when the number of harmonic components is unknown. In this paper, we combine the super-resolution of real-value MUSIC algorithm and the step-by-step estimation of step-by-step LSE algorithm, a step-by-step MUSIC algorithm is proposed. Our algorithm has good estimation performance when the number of components is known. When the number of components is unknown, we use the algorithm to estimate a model set, and estimate the number of harmonic components with the Bayesian Information Criterion (BIC). Finally, the performance of the proposed method is verified by simulation experiments.

1. 引言

谐波信号是一种特殊的随机信号，它的分析与处理是统计信号处理中的一个典型问题。其参数估计在很多领域都有广泛的应用，包括：雷达目标波达角估计 [1]、传感器阵列信号处理 [2]、图像处理 [3] 等。根据背景噪声的复杂程度可以分为加性噪声中的谐波参数估计与乘性和加性噪声中的谐波参数估计两大类。当加性噪声为白噪声时，已有一些算法如超逼近方法 [4]、旋转不变子空间(Estimation of Signal Parameters via Rotational Invariance Techniques, ESPRIT)方法 [5]、预滤波ESPRIT方法 [6]、多重信号分类(Multiple Signal Classification, MUSIC)方法 [7] 和高阶统计量方法 [8] [9] [10] 等，这些算法都具有高分辨率，可以取得良好的参数估计效果。

谐波参数估计问题中，分量数是谐波信号的一个重要参数，它是建立谐波模型以及估计频率等参数的基础。考察具有代表性的MUSIC算法，准确的分量数是发挥其超分辨率的前提，而对分量数的欠估计或过估计，会导致谱峰的估计偏差。

本文考虑MUSIC算法和Prasad等 [11] 提出的逐步LSE算法的逐步估计相结合，提出了一种基于改进MUSIC算法的加性噪声中的谐波参数估计算法，该算法不需要预先判定谐波分量数，同时在低信噪比、样本量小的情况下依然保持较高的频率分辨能力。利用该算法可逐步的估计出一个模型集，并利用模型选择方法中的BIC准则估计谐波分量数。最后通过仿真实验对所提方法的性能进行了验证。

2. 实谐波模型

实际应用中的谐波信号通常为多个实值正弦信号与噪声信号的叠加，即

$x\left(t\right)={\displaystyle \underset{k=1}{\overset{p}{\sum }}{A}_k\cos \left({\omega }_{k}t+{\phi }_k\right)}+n\left(t\right)$ (1)

其中，p称为谐波的分量个数， $A_k,{\omega }_k,{\phi }_k$ 分别为第k个谐波分量的振幅、频率和相位； $t=1,\cdots ,N$。N为采样数； $n\left(t\right)$ 为加性噪声。

利用三角函数和差公式可将(1)式表示为

$x\left(t\right)={\displaystyle \underset{k=1}{\overset{p}{\sum }}\left[a_k\cos \left({\omega }_{k}t\right)+b_k\sin \left({\omega }_{k}t\right)\right]}+n\left(t\right)$ (2)

其中 $a_k=A_k\cos \left({\phi }_k\right)$， $b_k=-A_k\sin \left({\phi }_k\right)$。

关于上述谐波信号模型的假设如下：

AS1： ${\omega }_k\in \left(-\pi ,0\right)\cup \left(0,\pi \right),k=1,\cdots ,p$，并且 ${\omega }_k$ 互不相等；

AS2： ${\phi }_k$ 是 $\left(-\pi ,\pi \right]$ 上的确定性常数；

AS3：加性噪声 $\left\{n\left(t\right),t=1,\cdots ,N\right\}$ 为i.i.d。零均值且方差为 ${\sigma }^2$ 的高斯随机变量序列。

为分析方便，定义一个长度为 $m\left(m>2p\right)$ 的观测量和噪声量分别表示为观测向量 $y\left(t\right)$ 和噪声向量 $n\left(t\right)$ ：

$\{\begin{array}{l}y\left(t\right)={\left[y\left(t\right),y\left(t+1\right),\cdots ,y\left(t+m-1\right)\right]}^{\text{T}}\\ n\left(t\right)={\left[n\left(t\right),n\left(t+1\right),\cdots ,n\left(t+m-1\right)\right]}^{\text{T}}\end{array}$ (3)

则观测向量可表示为：

$y\left(t\right)=C\left(\omega \right)s\left(t\right)+n\left(t\right)$ (4)

式中：

$C\left(\omega \right)={\left[B\left({\omega }_1\right),B\left({\omega }_2\right),\cdots ,B\left({\omega }_p\right)\right]}_{m\times 2p}$ (5)

$B\left({\omega }_k\right)={\left[\begin{array}{cccc}1& \cos \left({\omega }_k\right)& \cdots & \cos \left[\left(m-1\right){\omega }_k\right]\\ 0& -\sin \left({\omega }_k\right)& \cdots & -\sin \left[\left(m-1\right){\omega }_k\right]\end{array}\right]}^{\text{T}}$ (6)

$s\left(t\right)={\left[\begin{array}{c}{A}_1\cos \left({\omega }_{1}t+{\phi }_1\right)\\ A_1\sin \left({\omega }_{1}t+{\phi }_1\right)\\ ⋮\\ A_p\cos \left({\omega }_{p}t+{\phi }_p\right)\\ A_p\sin \left({\omega }_{p}t+{\phi }_p\right)\end{array}\right]}_{2p\times 1}$ (7)

3. 基于逐步MUSIC算法的谐波参数估计

文献 [12] 提出了一种实数形式的MUSIC算法，相较于常规的MUSIC算法，无需将信号转换为复频率，降低了计算量和计算复杂度。

因为 $n\left(t\right)$ 与 $y\left(t\right)$ 相互独立，于是可对观测向量 $y\left(t\right)$ 的自相关矩阵 $R_{yy}=E\left\{y\left(t\right)y^{\text{T}}\left(t\right)\right\}$ 做特征值分解为：

$R_{yy}=C\left(\omega \right)R_{ss}{C}^{\text{T}}\left(\omega \right)+{\sigma }^2{I}_M={\displaystyle \underset{k=1}{\overset{M}{\sum }}{\lambda }_k{\mu }_k{\mu }_k^{\text{T}}}=U_s{\Sigma }_s{U}_s^{\text{T}}+{\sigma }^2{U}_n{U}_n^{\text{T}}$ (8)

其中 $R_{ss}=E\left\{s\left(t\right)s^{\text{T}}\left(t\right)\right\}$ 是 $s\left(t\right)$ 的自相关矩阵， $I_M$ 是 $m\times m$ 单位矩阵。 ${\lambda }_k$ 和 ${\lambda }_k$ 分别为 $R_{yy}$ 第k个特征值及其对应的特征向量，且满足 ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_{2p+1}=\cdots ={\lambda }_m={\sigma }^2$。 ${\Sigma }_s$ 为 $R_{yy}$ 的前2p个特征向量构成的对角矩阵。 $U_s$ 称为 $R_{yy}$ 的信号子空间，由 $R_{yy}$ 的前2p个特征向量构成； $U_n$ 称为 $R_{yy}$ 的噪声子空间，由 $R_{yy}$ 的后 $m-2p$ 个特征矢量构成。

由于特征向量相互正交，可得 $U_s$ 与 $U_n$ 正交，将式(8)右乘 $U_n$ 得：

$R_{yy}{U}_n=C\left(\omega \right)R_{ss}{C}^{\text{T}}\left(\omega \right)+{\sigma }^2{U}_n={\sigma }^2{U}_n$ (9)

于是有：

$C\left(\omega \right)R_{ss}{C}^{\text{T}}\left(\omega \right)U_n=O\Rightarrow U_n^{\text{T}}C\left(\omega \right)R_{ss}{C}^{\text{T}}\left(\omega \right)U_n=O$ (10)

进而可以得到：

$C^{\text{T}}\left(\omega \right)U_n=O$ (11)

即 $C^{\text{T}}\left(\omega \right)$ 与 $U_n$ 正交。

将式(5)代入式(11)即有：

$B^{\text{T}}\left(\omega \right)U_n=O_{2\times \left(m-2p\right)},\omega ={\omega }_1,\cdots ,{\omega }_p$ (12)

显然，当 $\omega \ne {\omega }_1,\cdots ,{\omega }_p$ 时， $B^{\text{T}}\left(\omega \right)U_n\ne O$。

将式(12)改写为标量形式，即定义如下实值MUSIC伪谱函数：

$P\left(\omega \right)=\frac{1}{\text{tr}\left[B^{\text{T}}\left(\omega \right)U_n{U}_n^{\text{T}}B\left(\omega \right)\right]}$ (13)

于是，若谐波模型具有p个分量，则可选取实值MUSIC伪谱函数 $P\left(\omega \right)$ 所对应的p个峰值点作为频率的估计。

在求解得各谐波频率的估计后，可将非线性的信号模型转为简单的线性模型，进而采用最小二乘方法估计各谐波振幅。

引入以下记号：

$Y={\left[y\left(1\right),y\left(2\right),\cdots ,y\left(N\right)\right]}^{\text{T}}$

$\alpha ={\left(a_1,b_1,\cdots ,a_p,b_p\right)}^{\text{T}}$

$E={\left[n\left(1\right),n\left(2\right),\cdots ,n\left(N\right)\right]}^{\text{T}}$

$D\left({\omega }_k\right)={\left[\begin{array}{cccc}\cos \left({\omega }_k\right)& \cos \left(2{\omega }_k\right)& \cdots & \cos \left(N{\omega }_k\right)\\ \sin \left({\omega }_k\right)& \sin \left(2{\omega }_k\right)& \cdots & \sin \left(N{\omega }_k\right)\end{array}\right]}^{\text{T}}$

$G\left(\omega \right)={\left[D\left({\omega }_1\right),D\left({\omega }_2\right),\cdots ,D\left({\omega }_p\right)\right]}_{m\times 2p}$

则谐波模型(2)可改写为如下的矩阵形式：

$Y=G\left(\omega \right)\alpha +E$ (14)

于是，在求解得频率 $\omega$ 的估计值 $\hat{\omega }$ 后，上式的最小二乘解为：

$\hat{\alpha }={{\left[G^{\text{T}}\left(\omega \right)G\left(\omega \right)\right]}^{-1}{G}^{\text{T}}\left(\omega \right)Y|}_{\omega =\hat{\omega }}$ (15)

在实际应用中，通常会出现低信噪比、样本量小或频率差小的情形，此时信号往往会被淹没掉，导致一些谱峰消失。此外，MUSIC算法需要给定准确的分量数才可以发挥其超分辨率，在分量数未知的情况下，若估计的分量数比实际的要少(即欠估计)，进行信号子空间和噪声子空间的划分时，由于部分信号子空间会被错分到噪声子空间中去，导致信号子空间的维数降低，从而也会使得一些谱峰消失。同样的，若估计的分量数比实际的要多(即过估计)，那么部分噪声子空间会被错分到信号子空间中去，导致噪声子空间的维数降低，信号子空间维数变多，此时谱峰会增加，即出现一些伪峰。

为在分量数未知的情况下也可以取得良好的凸峰效果，文献 [13] 考虑利用自相关矩阵的全部特征向量，利用权值控制不同特征向量对特征谱的作用，而适当的权系数可以减弱信号特征向量的作用，加弱权值的信号特征向量即相当于常规的MUSIC方法中只利用噪声特征向量的情形。

定义数据自相关矩阵的加权特征向量矩阵为：

$U=\left[\frac{1}{{\lambda }_1^c}{\mu }_1,\frac{1}{{\lambda }_2^c}{\mu }_2,\cdots ,\frac{1}{{\lambda }_m^c}{\mu }_M\right]$ (16)

其中，指数c用来调节加权的效果。

于是新的MUSIC空间谱表示为：

$H\left(\omega \right)=\frac{1}{\text{tr}\left[B^{\text{T}}\left(\omega \right)U{U}^{\text{T}}B\left(\omega \right)\right]}$ (17)

利用上式，我们结合Prasad等 [11] 提出的逐步LSE算法，提出了一个简单的逐步估计算法，在谐波分量数未知情况下也可以取得良好的估计效果。

对原始数据 $Y$ 的自相关矩阵做特征值分解，并构造加权左奇异矩阵 $U^{\left(1\right)}$。接着搜索对应的MUSIC空间谱的最大峰值，并利用式(15)求解得到振幅的估计，即：

${\hat{\omega }}_1=\underset{\omega }{\text{argmax}}\frac{1}{\text{tr}\left[B^{\text{T}}\left(\omega \right)U^{\left(1\right)}{U}^{\left(1\right)}{}^{{}^{\text{T}}}B\left(\omega \right)\right]}$ (18)

${\hat{\alpha }}^{\left(1\right)}={{\left[D^{\text{T}}\left(\omega \right)D\left(\omega \right)\right]}^{-1}{D}^{\text{T}}\left(\omega \right)Y|}_{\omega ={\hat{\omega }}_1}$ (19)

则有 ${\hat{\alpha }}^{\left(1\right)}={\left[{\hat{a}}_1^{\left(1\right)},{\hat{b}}_1^{\left(1\right)}\right]}^{\text{T}}$， $\left({\hat{a}}_1^{\left(1\right)},{\hat{b}}_1^{\left(1\right)},{\hat{\omega }}_1\right)$ 是 $\left(a_1^0,b_1^0,{\omega }_1^0\right)$ 的估计量。接着对原始数据做调整，即如下残差序列：

$Y^{\left(1\right)}=Y-B\left({\hat{\omega }}_1\right){\hat{\alpha }}^{\left(1\right)}$ (20)

对调整后的数据 $Y^{\left(1\right)}$ 的自相关矩阵做特征值分解，并构造加权特征向量矩阵 $U^{\left(2\right)}$。接着搜索对应的MUSIC空间谱的最大峰值，即：

${\hat{\omega }}_2=\underset{\omega }{\text{argmax}}\frac{1}{\text{tr}\left[B^{\text{T}}\left(\omega \right)U^{\left(2\right)}{U}^{\left(2\right)}{}^{{}^{\text{T}}}B\left(\omega \right)\right]}$ (21)

根据已有的频率估计量 ${\hat{\omega }}^{\left(2\right)}=\left({\hat{\omega }}_1,{\hat{\omega }}_2\right)$，利用式(22)求解前两个分量的振幅的估计量：

${\hat{\alpha }}^{\left(2\right)}={{\left[G^{\text{T}}\left(\omega \right)G\left(\omega \right)\right]}^{-1}{G}^{\text{T}}\left(\omega \right)Y^{\left(1\right)}|}_{\omega ={\hat{\omega }}^{\left(2\right)}}$ (22)

则有 ${\hat{\alpha }}^{\left(2\right)}={\left[{\hat{a}}_1^{\left(2\right)},{\hat{b}}_1^{\left(2\right)},{\hat{a}}_2^{\left(2\right)},{\hat{b}}_2^{\left(2\right)}\right]}^{\text{T}}$。接着对原始数据做调整，即如下残差序列：

$Y^{\left(2\right)}=Y-C\left({\hat{\omega }}^{\left(2\right)}\right){\hat{\alpha }}^{\left(2\right)}$ (23)

若谐波分量数p已知，我们将继续这一过程到第p步。加性噪声方差 ${\sigma }^2$ 的估计给出为：

${\hat{\sigma }}_p^2=\frac{1}{N}{\left[Y^{\left(p\right)}\right]}^{\text{T}}{Y}^{\left(p\right)}$ (24)

若p未知，我们将拟合到一个q阶模型，并给出对应的加性噪声方差的估计。然而q可能不等于p，且通常是一个大于p的整数。于是，我们可根据BIC准则选择分量数，针对加性高斯白噪声中的谐波模型，可给出BIC的形式如下：

$\text{BIC}\left(k\right)=N\ln {\hat{\sigma }}_k^2+3k\ln N$ (25)

对 $k=1,\cdots ,q$，选取使得上式最小的k作为谐波分量数的估计。

4. 仿真实验

在上节中，提出了加性噪声中谐波恢复的一种新方法——逐步MUSIC算法。本节将采用Mont Carol方法进行仿真实验，以展示我们算法的有效性。所有仿真实验均有MATLAB编写运行。

实验一 考虑以下谐波信号：

$x_1\left(t\right)=\underset{k=1}{\overset{3}{{\displaystyle \sum }}}\left[a_k\cos \left({\omega }_{k}t\right)+b_k\sin \left({\omega }_{k}t\right)\right]+n(\; t\; )$

其中， $a_1=5$， $b_1=3$， ${\omega }_1=0.4$ ； $a_2=4$， $b_2=2$， ${\omega }_2=0.43$ ； $a_3=3$， $b_3=1.8$， ${\omega }_3=0.5$。 $n\left(t\right)$ 是零均值的高斯白噪声。

在谐波分量数已知的情况下，为考察本文算法在不同噪声情形下的性能，取不同信噪比分别为0和5，其中信噪比定义为：

$\text{SNR}=10{\log }_{10}\frac{\min \left\{a_k^2+b_k^2\right\}}{2{\sigma }^2}$

为考察本文算法在不同样本下的性能，取不同样本量分别为。

实验参数设置上，取观测向量长度为0.4倍样本量，频率估计采用峰值搜索法，取步长为0.0001。逐步MUSIC算法的权重参数c取0.5。进行200次Mont Carol实验。

Table 1. Frequency estimation results with a SNR of 0 dB

表1. 信噪比为0 dB时的频率估计结果

Table 2. Frequency estimation results with a SNR of 5 dB

表2. 信噪比为5 dB时的频率估计结果

在谐波分量数已知的情况下，从实验结果(表1和表2)得知，逐步MUSIC算法在高信噪比、样本量大时具有和实值MUSIC算法相当的估计性能；而在低信噪比、样本量小时，逐步MUSIC算法的估计性能超过了实值MUSIC算法。

实验二 考虑以下谐波信号：

其中，，，；，，；，，。是零均值的高斯白噪声。

在谐波分量数未知的情况下，我们将对拟合谐波模型，取不同信噪比分别为−5 dB和0 dB，取样本量从逐次递增10个到，最后采用BIC准则估计谐波分量数。

实验参数设置同实验一，进行500次Mont Carol实验。

Figure 1. Probability of correct estimation of the number of components under different SNR and sample sizes

图1. 不同信噪比和样本量情况下的分量数估计正确的概率

从图1的实验结果得知，在信噪比为0 dB时，逐步MUSIC算法相较于实值MUSIC算法具有更好的估计性能，且随着样本量增大分量数估计正确的概率接近于1。而在信噪比为−5 dB时，逐步MUSIC算法仍具有不错的估计准确率，但此时实值MUSIC算法已经完全失效。

综合实验一与实验二，可以得出本文所提出的逐步MUSIC算法相较于实值MUSIC算法具有更好的恢复效果。

5. 结论

本文将加性噪声中的实值MUSIC算法与逐步LSE算法相结合，提出了一种新的谐波恢复方法——逐步MUSIC算法，在谐波分量数未知的情况下，可利用该算法逐步估计出一个模型集，进而运用BIC准则估计出谐波分量数。仿真实验说明，与原有的实值MUSIC算法相比，逐步MUSIC算法具有更好的谐波恢复效果。