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

陳毅濱

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

0 引言

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

1 Kay 算法分析

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

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

如果把噪聲添加進(jìn)去觀測(cè)信號(hào)可以表示為：

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

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

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

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

通過(guò)抽樣定律可以得到ω的值比π 要小，也就是說(shuō)式（4）的計(jì)算并沒(méi)有相位模糊，所以能估計(jì)信號(hào)的頻率。

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

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

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

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

2 正弦波信號(hào)頻率和初相估計(jì)

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

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

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

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

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

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

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

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

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

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

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

同時(shí)，還能得到初相估計(jì)的最大似然方程：

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

3 算法改進(jìn)

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

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

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

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

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

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

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

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

4 測(cè)試結(jié)果

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

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

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

圖1 對(duì)比測(cè)試結(jié)果

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

5 結(jié)語(yǔ)

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