一種低信噪比下正弦波頻率估計方法*

陳毅濱

（上海郵電設(shè)計咨詢研究院有限公司，上海 200092）

0 引言

正弦波信號頻率估計在數(shù)字信號處理領(lǐng)域中的應(yīng)用非常廣泛，也是國內(nèi)外學(xué)者研究的熱點[1-4]。目前，研究正弦波信號頻率估計的方法比較多，如子空間算法。此種方法優(yōu)點是分辨率高，但是運算量大。離散傅里葉變換分析法的優(yōu)點是運算量小，但是在小樣本下的情況估計精度不高[5-6]。Rife 算法雖然估計精度有所提高，但是要求在信噪比高的時候估計誤差才能和克拉美-羅限（Cramer-Rao Lower Bound，CRLB）差不多[7]。Kay 算法對信號長度要求較高，如果采樣點數(shù)不多，需要較高的信噪比門限。此類算法在低信噪比下的估計精度都比較低[8]。為了能夠?qū)崿F(xiàn)在低信噪比下對正弦波信號頻率估計，文章提出了一種新的方法。此方法利用Kay算法的特點，經(jīng)過理論分析和計算得到正弦波兩個參數(shù)的最大最大似然估計。這兩個參數(shù)是頻率和初相，同時得到了修正權(quán)值矩陣。這樣在低信噪比下和采樣點數(shù)不多的情況下可完成正弦波頻率的估計。

1 Kay 算法分析

設(shè)離散正弦波信號為：

式中：A代表幅值；φ代表初始相位；ω代表數(shù)字角頻率。

如果把噪聲添加進去觀測信號可以表示為：

式中：y(t)代表高斯白噪聲，其均值是0、方差是

w(t)的無模糊相位為：

式中：λ(t)為相位噪聲。

w(t)的一階相位差為：

通過抽樣定律可以得到ω的值比π 要小，也就是說式（4）的計算并沒有相位模糊，所以能估計信號的頻率。

式中：ψ是一階相位差分矢量；R代表的意義是單位向量，表達式為R=[1 1 … 1]N；ε的意義是相位差噪聲矢量[9]，表達式如式（6）所示。

根據(jù)Kay 算法的特點，λ(t)可以認(rèn)為是高斯噪聲，其均值可以認(rèn)為接近0、其方差是δy2/2A。因此，就可以計算頻率估計值：

式中：RNJopt=[J(1)J(2)…J(T-1)]，其最小二乘估計的加權(quán)行向量；Jopt表示最優(yōu)權(quán)值矩陣，表達式為Jopt=D-1。

D為ε的協(xié)方差矩陣，表達式為D=G(εεN)，還可以表示為：

2 正弦波信號頻率和初相估計

文章使用的是最大似然估計法[10]，是一種統(tǒng)計方法。假設(shè)w(t)的似然函數(shù)為：

式中：S[·]表示取信號的實部；X[·]表示去信號的虛部。如果把直角坐標(biāo){S[w(t)],X[w(t)]}變成極坐標(biāo)[|w(t)|,η(t)]，得到：

w(t)的模值可以表示為：

實際應(yīng)用中，頻率和初相都是固定不變的，所以y'(t)和y(t)的統(tǒng)計特性是一樣的，也就是說的統(tǒng)計特征與頻率及初相這兩參數(shù)沒有任何關(guān)系。因此，能把似然函數(shù)進一步簡化為：

式中：F|w(t)|表示|w(t)|的相關(guān)函數(shù)。通過上述分析計算可以看出，w(t)的似然函數(shù)主要是由概率密度函數(shù)決定的[11]。

|w(t)|和相位聯(lián)合概率密度函數(shù)可以表示為：

|w(t)|的邊緣概率密度函數(shù)可以表示成：

式中：B[·]代表的意義是第一類零階修正貝塞爾函數(shù)[12]。在分析計算時，如果把ω 和φ0的值確定，并且λ(t)是η(t)中僅有的隨機量，也就是說l[λ(t)||w(t)|]=l[η(t)||w(t)|]。因此就可以把條件概率密度函數(shù)表示：

因此，似然函數(shù)為：

通過式（14）、式（15）和式（16），對數(shù)似然函數(shù)可表示為：

式中：F′|w(t)|的值和參數(shù)ω以及φ0沒有任何關(guān)系。依據(jù)最大似然估計的基本原理，頻率估計的最大似然方程可表示為：

同時，還能得到初相估計的最大似然方程：

因為正弦函數(shù)自身的特點，其是一個非線性函數(shù)，所以沒有辦法實現(xiàn)頻率和初相估計的閉合解。通過式（18）和式（19）可以看出，頻率和初相的估計同時和兩個參數(shù)有關(guān)，這兩個參數(shù)是η(t)和|w(t)|。通過前面的論述可知，Kay 算法只能使用η(t)信息，并不能使用|w(t)|信息，所以低信噪比下的頻率估計誤差較大。

3 算法改進

把余弦函數(shù)的泰勒級數(shù)代入再次進行計算，同時還把4 次方相去除，就可以得到：

通過式（21）可以看出，如果觀測信號的模值是確定的，則信號相位噪聲可以看成是服從高斯分布規(guī)律的[13]，因此就可以把其均值和方差表示出來。

因為ψ(t)與采樣點w(t)和w(t+1)是相關(guān)的，所以ψ(t)的權(quán)值是由相位噪聲λ(t)和λ(t+1)決定的。|w(t)|和|w(t+1)|越小，λ(t)和λ(t+1)越大。在加權(quán)的時候，需要賦予ψ(t)更大的權(quán)重，反過來也是成立的。

相位差噪聲矢量協(xié)方差矩陣可以表示成：

實際上，E[λ2(t)]和C[λ(t)]的值是相等的，所以可以把式（24）進一步簡化成：

通過式（25）可以看出，Dmn的值由參數(shù)δ2、A 以及|w(m)|決定。實際上，參數(shù)δ2、A 都是不確定的，在設(shè)定權(quán)值矩陣的時候，只需確定相位差噪聲相對大小就可以。所以，式（25）可以進一步簡化成：

由式（26）就能得到修正后的權(quán)值矩陣[14]。

本文算法估計正弦波評率主要分4 步進行：第一步是計算觀測信號的模值；第二步是通過式（24）得到Dmn的值；第三步是把Dmn代入式（7），計算最優(yōu)權(quán)值矩陣；第四步是使用最小二乘估計法實現(xiàn)估計正弦波頻率。

4 測試結(jié)果

如果不知道正弦波的相位、幅度和頻率信息，則CRLB可以表示成：

式中：ψn表示意義是采樣周期；SNR的意義是信噪比，值為

w(t)在進行理論計算的時是沒有模糊的。然而，如果λ(t+1)-λ(t)比較大，就會出現(xiàn)相位模糊的情況[15]。運用MATLAB 做對比仿真實驗，參數(shù)ω=0.125π，采樣點數(shù)N 的值設(shè)定為11，蒙特卡洛為500 次，初始相位是不定值，圖1 就是測試結(jié)果。

圖1 對比測試結(jié)果

需要說明的是，縱坐標(biāo)為101gMSE，MSE 表示的意義是頻率估計均方誤差。由圖1 可知，當(dāng)信噪比大于6 dB 情況下，兩種算法都比較接近CRLB，但是在信噪比較低時，本文算法的優(yōu)勢非常明顯，在信噪比為1 dB 時接近CRLB。兩種算法都存在信噪比門限，但是本文算法的誤差惡化程度相對來說非常低，因為本文算法根據(jù)觀測信號的模值可以實現(xiàn)自動調(diào)整權(quán)值矩陣。

5 結(jié)語

文章論述了低信噪比下正弦波頻率估計方法，在Kay 算法基礎(chǔ)上進一步改進，利用修正后的權(quán)值矩陣的優(yōu)點進行頻率估計。通過測試結(jié)果可以看出，本文算法的正弦波頻率估計誤差更小，特別是在低信噪比下優(yōu)勢明顯。雖然也有信噪比門限，但是和傳統(tǒng)的Kay 算法相比，大大降低了誤差惡化程度。