一种基于多观测向量稀疏表示的复正弦波频率估计方法

摘要

本发明属于信号处理领域，尤其涉及复正弦波频率估计。本发明提供一种高精度复正弦波频率估计方法。由于复正弦波信号在频域具有稀疏性，基于观测数据稀疏表示的模型，首先针对网格失配的情况，通过一阶泰勒展开来逼近真实的频率向量，得到多观测向量下的稀疏表示逼近模型。然后在该模型的基础上用正交匹配追踪算法计算得到网格上离真实频率值最近的点，并且利用最小二乘方法计算频率的修正值。最后得到修正后复正弦波频率的估计值。利用本发明的估计方法，修正后的估计性能有明显的提高，即使在粗糙的网格，也有很高的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高。

说明书

技术领域

本发明属于信号处理领域，尤其涉及复正弦波频率估计。

背景技术

正弦波频率估计算法广泛应用于雷达、通信以及电子对抗等信号处理领域。其中经典的算法有：多重信号分类(Multiple Signal Classification，MUSIC)算法、旋转不变子空间(Estimation of Signal Parameters via Rotational Invariance Technique，ESPRIT)算法等子空间类算法和最大似然估计类算法(Maximum Likelihood，ML)等。然而，基于子空间理论的频率估计方法虽然实现了超分辨侧向，但是一旦接收数据不足或信噪比比较低时，这类方法不能有效地区分信号子空间和噪声子空间，其性能会急剧下降。而最大似然估计类算法由于要进行复杂的多维搜索而不具有实用性。

近年来，随着压缩感知和稀疏表示技术的发展，许多基于稀疏表示的频率估计方法被提出。最具代表性的为l₁-SVD算法，它利用l₁范数来重构稀疏信号，并且在多观测向量的条件下通过奇异值分解(Singular Value Decomposition，SVD)来减小数据矩阵的规模以及降低噪声的影响。然而，这些基于稀疏表示的方法通常假设频率位于离散化网格上，一旦真实的频率值不在预先设定的网格上，即在网格失配的情况下，其估计性能会急剧恶化。另一方面，虽然更密集的网格理论上可以减小重构误差，但是太过密集的离散化网格会使得过完备化字典原子间高度相关。在其他领域，通用一阶泰勒展开来逼近真实字典来解决字典失配的问题，例如，在DOA估计中，利用稀疏贝叶斯推论(Off-Grid Sparse Bayesian Inference，OGSBI)来重构稀疏信号，可以估计不在网格上的来波方向角。类似地，该方法也可以用于复正弦波频率估计，但是这种方法的缺点是估计结果容易受初值的影响，且计算量较大。

发明内容

本发明的目的在于，提出一种基于多观测向量稀疏表示的复正弦波频率估计方法，在频率网格失配的情况下，能够对频率值进行修正，提高估计精度。

本发明基于观测数据稀疏表示的模型，首先针对网格失配的情况，通过一阶泰勒展开来逼近真实的频率向量，得到多观测向量下的稀疏表示逼近模型，然后用正交匹配追踪算法和最小二乘法计算复正弦波频率的估计值。

本发明的目的通过如下步骤实现：

S1、由长度为P的复正弦信号x(i)加白噪声得到接收数据其中，i＝1，…，P，K为复正弦波信号的个数，为第k个信号的复幅度，f_k∈[0，1)为第k个信号的归一化频率，k＝1，2，…，K，初始相位是在[0，2π]均匀分布的随机变量，v(i)是零均值、方差为的高斯白噪声，当i≠k时，和相互独立，

S2、S1所述接收数据y(i)的向量形式为y＝A(f)s₀+v，其中，y＝[y(1) y(2) … y(P)]^T，频率向量$>a\left(f_k\right)={\left(\begin{array}{cccc}1& e^{j2\pi f_k}& ···& e^{j2\pi (P-1)f_k}\end{array}\right)}^T,$>幅度向量s₀＝[α₁ α₂ … α_K]^T，v＝[v(1) v(2) … v(P)]^T；

S3、将S2所述接收数据y(i)的向量形式y写成稀疏表示的模型Y，具体如下：

S31、由S2所述接收数据y(i)，得到信号向量y(i)＝[y(i) y(i+1) … y(i+M-1)]^T，从而得到数据矩阵其中，L+M-1＝P，M为信号向量的长度，L为信号向量的个数；

S32、把S1所述归一化频率f_k在[0，1)的范围上过完备化为一个离散的网格N代表网格数且满足N＞＞K；

S33、将S31所述信号向量y(i)写成多观测向量的稀疏表示形式其中，过完备矩阵频率向量$>a\left({\tilde{f}}_n\right)={\left(\begin{array}{cccc}1& e^{j2\pi {\tilde{f}}_{\mathrm{rk}}}& ...& e^{j2\pi (M-1){\tilde{f}}_n}\end{array}\right)}^T,$>是一个K-行稀疏矩阵，即S中只有K行元素为非零，其他行的元素都为零，记所述K行元素的索引值为支撑集Λ，v(i)＝[v(i) v(i+1) … v(i+M-1)]^T；

S4、在S3所述稀疏表示的模型Y的基础上，用一阶泰勒展开来逼近真实的频率向量$>a\left(f_k\right)\approx a\left({\tilde{f}}_{n_k}\right)+b\left({\tilde{f}}_{n_k}\right)(f_k-{\tilde{f}}_{n_k}),$>其中，存在第k个真实频率$>f_k\notin \{{\tilde{f}}_1,...,{\tilde{f}}_N\},$>且是网格上距离f_k最近的点，是的一阶导数，n_k∈{1，…，N}；

S5、根据S4所述a(f_k)构造过完备字典$>B=\left(\begin{array}{cccc}b\left({\tilde{f}}_1\right)& b\left({\tilde{f}}_2\right)& ...& b\left({\tilde{f}}_N\right)\end{array}\right),$>记Δ＝diag(δ)，δ＝[δ₁ δ₂ … δ_N]^T，

当$>f_k\ne {\tilde{f}}_{n_k}$>时，则$>{\delta }_n=f_k-{\tilde{f}}_{n_k},$>

当$>f_k={\tilde{f}}_{n_k},$>则δ_n＝0，

因此，观测模型可以写成$>Y=(A+\mathrm{B\Delta })S+V=\left[A\right|B]\left(\begin{array}{c}S\\ \mathrm{\Delta S}\end{array}\right)+V=[A\left|B\right]\left(\begin{array}{c}S\\ Q\end{array}\right)+V;$>

S6、利用S5所述观测模型Y通过正交匹配追踪(orthogonal matching pursuit，OMP)算法和最小二乘方法求解出真实频率值的估计值，具体为：

S61、用正交匹配追踪算法求解出网格上离真实频率值最近的点，具体如下：

S611、初始化记迭代次数t＝1，残差矩阵R＝Y，支撑集其中0_N×L表示维数为N×L的全零矩阵；

S612、计算C＝A^HR，λ＝arg max||C(i，：)||₂，Λ＝Λ∪λ，其中(·)^H表示矩阵的共轭转置，||C(i，：)||₂表示矩阵C第i行的l₂范数，arg max(·)表示最大值的索引值；

S62、用最小二乘方法计算S和Q在支撑集上元素的估计值和并更新残差矩阵$>R=Y-A_{\Lambda }{\hat{S}}_{\Lambda }-B_{\Lambda }{\hat{Q}}_{\Lambda },$>更新迭代次数t＝t+1；

S63、当迭代次数t≤K时，返回步骤S612，依次循环迭代；一旦t＞K，停止迭代，输出和把和的对应元素相除，得到维数为K×L的矩阵求矩阵Γ每一行的平均值，得到向量δ在支撑集上元素的估计值

S64、设S63所述在网格上的索引值为n₁，...，n_K，则修正后的频率估计值为$>{\hat{f}}_k={\tilde{f}}_{n_k}+{\hat{\delta }}_{n_k}.$>

本发明的有益效果是：

本发明可在粗糙的网格上进行复正弦波频率值的精确估计，避免了密集的网格带来的高计算量，提高了估计精度。利用本发明的估计方法，修正后的估计性能有明显的提高，即使在粗糙的网格，也有很高的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高。

附图说明

图1是本发明方法的流程图。

图2是本发明方法与其他频率估计方法的均方根误差随信噪比变化曲线图。

图3是本发明方法与其他频率估计方法的均方根误差随信号向量个数L变化曲线图。

具体实施方式

下面结合实施例和附图，详细说明本发明的技术方案。

如图1所示，一种基于多观测向量稀疏表示的高精度频率估计方法包括如下步骤：

S1、由长度为P的复正弦信号x(i)加白噪声得到接收数据$>y\left(i\right)=x\left(i\right)+v\left(i\right)=\underset{k=1}{\overset{K}{\Sigma }}{\alpha }_k{e}^{j2\pi f_k(i-1)}+v\left(i\right),$>其中，i＝1，…，P，K为复正弦波信号的个数，为第k个信号的复幅度，f_k∈[0，1)为第k个信号的归一化频率，k＝1，2，…，K，初始相位是在[0，2π]均匀分布的随机变量，v(i)是零均值、方差为的高斯白噪声，所述v(i)与信号相互独立，当i≠k时，和相互独立，

S2、S1所述接收数据y(i)的向量形式为y＝A(f)s₀+v，其中，y＝[y(1) y(2) … y(P)]^T，频率向量$>a\left(f_k\right)={\left(\begin{array}{cccc}1& e^{j2\pi f_k}& ...& e^{j2\pi (P-1)f_k}\end{array}\right)}^T,$>幅度向量s₀＝[α₁ α₂ … α_K]^T，v＝[v(1) v(2) … v(P)]^T，进行频率估计的目的是根据接收的P×1维数据，估计K个信号的归一化频率{f_k，k＝1，2，…，K}；

S3、将S2所述接收数据y(i)的向量形式y写成稀疏表示的模型Y，具体如下：

S31、由S2所述接收数据y(i)，得到信号向量y(i)＝[y(i) y(i+1) … y(i+M-1)]^T，从而得到数据矩阵其中，L+M-1＝P，M为信号向量的长度，L为信号向量的个数；

S32、把S1所述归一化频率f_k在[0，1)的范围上过完备化为一个离散的网格N代表网格数且满足N＞＞K；

S33、将S31所述信号向量y(i)写成多观测向量的稀疏表示形式其中，过完备矩阵频率向量$>a\left({\tilde{f}}_n\right)={\left(\begin{array}{cccc}1& e^{j2\pi {\tilde{f}}_{\mathrm{rk}}}& ...& e^{j2\pi (M-1){\tilde{f}}_n}\end{array}\right)}^T,$>是一个K-行稀疏矩阵，即S中只有K行元素为非零，其他行的元素都为零，记所述K行元素的索引值为支撑集Λ，v(i)＝[v(i) v(i+1) … v(i+M-1)]^T；

S4、在S3所述稀疏表示的模型Y的基础上，用一阶泰勒展开来逼近真实的频率向量$>a\left(f_k\right)\approx a\left({\tilde{f}}_{n_k}\right)+b\left({\tilde{f}}_{n_k}\right)(f_k-{\tilde{f}}_{n_k}),$>其中，$>f_k\notin \{{\tilde{f}}_1,...,{\tilde{f}}_N\}$>是网格上距离f_k最近的点，k∈{1，…，K}且$>{\tilde{f}}_{n_k},n_k\in \{1,...,N\},b\left({\tilde{f}}_{n_k}\right)=a^{\prime }\left({\tilde{f}}_{n_k}\right)$>是的一阶导数；

S5、根据S4所述a(f_k)构造过完备字典$>B=\left(\begin{array}{cccc}b\left({\tilde{f}}_1\right)& b\left({\tilde{f}}_2\right)& ...& b\left({\tilde{f}}_N\right)\end{array}\right),$>记Δ＝diag(δ)，δ＝[δ₁ δ₂ … δ_N]^T，

当$>f_k\ne {\tilde{f}}_{n_k}$>时，则$>{\delta }_n=f_k-{\tilde{f}}_{n_k},$>

当则δ_n＝0，

因此，观测模型可以写成$>Y=(A+\mathrm{B\Delta })S+V=\left[A\right|B]\left(\begin{array}{c}S\\ \mathrm{\Delta S}\end{array}\right)+V=[A\left|B\right]\left(\begin{array}{c}S\\ Q\end{array}\right)+V;$>

S6、利用S5所述观测模型Y通过正交匹配追踪(orthogonal matching pursuit，OMP)算法和最小二乘方法求解出真实频率值的估计值，具体为：

S61、用正交匹配追踪算法求解出网格上离真实频率值最近的点，具体如下：

S611、初始化记迭代次数t＝1，残差矩阵R＝Y，支撑集其中0_N×L表示维数为N×L的全零矩阵；

S612、计算C＝A^HR，λ＝arg max||C(i，：)||₂，Λ＝Λ∪λ，其中(·)^H表示矩阵的共轭转置，||C(i，：)||₂表示矩阵C第i行的l₂范数，arg max(·)表示最大值的索引值；

S62、用最小二乘方法计算S和Q在支撑集上元素的估计值和并更新残差矩阵$>R=Y-A_{\Lambda }{\hat{S}}_{\Lambda }-B_{\Lambda }{\hat{Q}}_{\Lambda },$>更新迭代次数t＝t+1；

S63、当迭代次数t≤K时，返回步骤S612，依次循环迭代；一旦t＞K，停止迭代，输出和把和的对应元素相除，得到维数为K×L的矩阵求矩阵Γ每一行的平均值，得到向量δ在支撑集上元素的估计值

S64、设S63所述在网格上的索引值为n₁，...，n_K，则修正后的频率估计值为$>{\hat{f}}_k={\hat{f}}_{n_k}+{\hat{\delta }}_{n_k}.$>

下面结合具体实施例来进行说明。

实施例1、估计值的均方根误差随信噪比变化

采用接收数据的长度为P＝207，信号向量长度为M＝8。两个相同功率的复正弦波的归一化频率为[0.253，0.418]，即K＝2。为了使真实的复正弦波频率不落在网格上，取离散化网格为{0，0.01，…，0.99}，间隔0.01。参考信噪比SNR从-5dB到20dB变化，间隔为5dB，每个信噪比进行1000次蒙特卡洛实验。

根据不同信噪比下的接收数据y(i)得到稀疏表示模型Y＝AS+V；

利用泰勒一阶展开对真实的频率向量进行逼近，进而得到多观测向量下的稀疏表示逼近模型$>Y=\left[A\right|B]\left(\begin{array}{c}S\\ Q\end{array}\right)+V;$>

通过正交匹配追踪算法计算得到网格上离真实频率值最近的点：

通过最小二乘方法求解上一步骤波达方向角的修正值

得到修正后的频率估计值

按照本发明的方法估计得到的频率值的均方根误差随信噪比变化的曲线如图2所示。图2可以看到，利用本发明的估计方法，修正后的估计性能有明显的提高，即使在粗糙的网格，且信噪比为0dB时都能达到0.005以内的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。而且本发明方法所需的计算时间大约为0.005秒左右，而其他估计方法的计算时间均在1秒左右，说明了本发明估计方法的高效性。用均方根误差(RMSE)来评估各算法的性能，其定义为：其中，Mon为蒙特卡洛实验次数，和f_k分别代表第m次蒙特卡洛实验估计得到的第k个频率值和第k个真实频率。

实施例2、估计值均方根误差随信号向量个数L变化

信号向量长度M＝8，信噪比固定为10dB。两个相同功率的复正弦波的归一化频率为[0.253，0.418]，即K＝2。为了使真实的复正弦波频率不落在网格上，取离散化网格为{0，0.01，…，0.99}，间隔0.01。信号向量个数L从20到200变化，间隔为20，向量个数L进行1000次蒙特卡洛实验。用均方根误差(RMSE)来评估各算法的性能，其定义为：其中，Mon为蒙特卡洛实验次数，和f_k分别代表第m次蒙特卡洛实验估计得到的第k个频率值和第k个真实频率。

根据不同信号向量个数L下的接收数据y(i)得到稀疏表示模型Y＝AS+V；

利用泰勒一阶展开对真实的频率向量进行逼近，进而得到多观测向量下的稀疏表示逼近模型$>Y=\left[A\right|B]\left(\begin{array}{c}S\\ Q\end{array}\right)+V;$>

通过正交匹配追踪算法计算得到网格上离真实频率值最近的点：

通过最小二乘方法求解上一步骤波达方向角的修正值

得到修正后的频率估计值

按照本发明的方法估计得到的频率值的均方根误差随信号向量个数L变化的曲线如图3所示。图3可以看到，利用本发明的估计方法，修正后的估计性能有明显的提高，即使在粗糙的网格，且信号向量个数L＝40时都能达到0.001以内的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。而且本发明方法所需的计算时间大约为0.005秒左右，而其他估计方法的计算时间均在1秒左右，说明了本发明估计方法的高效性。