一種基于窄帶自相關的實信號頻率估計算法

曹 燕， 王一歌， 李欣雯， 趙明劍， 丁泉龍

(華南理工大學電子與信息學院，廣州 510641)

正弦信號分為復正弦信號和實正弦信號。實正弦信號在頻域中除了含有正頻率外，還有“負頻率”[1]，會在一定程度上造成頻譜的混疊，因此其頻率估計比復正弦信號復雜。在實際工程應用中，涉及的信號多數為被高斯白噪聲污染的實正弦信號，由于不能直接應用復正弦信號的相關方法，因此如何從混有高斯白噪聲的實信號樣點中提取信號的頻率一直以來都備受關注。

正弦信號的頻率估計有基于時域特性和基于變換域特征的兩種方法。基于時域的方法相對簡單，特別是利用少量的點在實時計算上有優勢[2]。基于變換域的頻率估計方法通常是考慮信號的頻域特性，在估計性能上有一定的優勢，但實現較為復雜。由于正弦信號的自相關仍然是同頻率的正弦信號，且去除部分噪聲影響，所以，在時域中基于自相關的頻率估計算法頗多，例如Pisarenko諧波分解(pisarenko harmonic decomposer，PHD)[3]算法，改進協方差(modified covariance，MC)[4]算法、以及基于自相關的一些改進頻率算法[5-10]，如利用更多的自相關系數[8]或是多步自相關[9-10]來提升估計性能。文獻[4]把MC運用到自相關序列，提出了基于自相關線性預測誤差最小的估計器。

基于頻域的頻率估計方法一般先進行頻率粗估計，再進行頻率細化。粗估計是尋找幅度譜的譜峰，可利用(快速傅里葉變換,fast fourier transform,FFT)實現，但是會出現 “柵欄”效應[11]，導致頻率估計精度降低。為了提高精度，可對原始序列補零后再做FFT。頻率細化的一種方法是基于插值的方法[12-16]。Quinn[12]提出利用最大譜峰和相鄰次大譜峰對應的傅里葉系數進行插值來實現復正弦信號的頻率估計。文獻[15]在不對原始實信號補零的情況下，利用最大譜峰和相鄰次大譜峰三點來內插，得到頻率的精估計。

現結合時域和頻域自相關的特點，提出一種基于窄帶自相關的實信號頻率估計算法。該算法在頻域進行譜峰搜索后，利用信號的窄帶功率譜來計算自相關，以消除窄帶外的噪聲對頻率估計的影響，然后再用簡單的時域自相關的MC算法來得到頻率估計的閉式解。最后給出與時域自相關MC算法，以及頻域全頻帶的自相關MC算法和頻域三點內插算法[15]相比較的仿真結果。

1 基于自相關的實信號頻率估計

1.1 實信號的自相關函數

假設原實正弦信號為

s(n)=acos(ω0n+φ)，n=0,1,…，N-1

(1)

式(1)中：a、ω0、φ分別表示信號的振幅、頻率和相位。混入噪聲的信號為

x(n)=s(n)+z(n)

(2)

式(2)中：z(n)為高斯白噪聲，其均值為0，方差為σ2。它的自相關函數為

(3)

式中：τ為自相關函數的下標，τ=0，1,…,N-1，有噪聲情況下，實信號的自相關函數為

(4)

當N足夠大時，自相關函數的期望為

(5)

式(5)中：δτ,0為克羅內克(Kronecker)函數：

(6)

可以看出，在有噪聲或者無噪聲的情況下，正弦信號的自相關函數仍為與原信號同頻率的正弦信號。而高斯白噪聲的自相關函數的期望只有在零點位置存在，其余位置均為0。所以，利用非零序號的自相關可以估計原信號的頻率，同時還可以在一定程度上降低噪聲的影響。

1.2 基于原信號的MC算法

由于正弦信號s(n)具有線性預測特性：s(n)=2cosω0s(n-1)-s(n-2)，也可以時延后寫成：

s(n+2)=2cosω0s(n+1)-s(n)

(7)

基于這種線性預測，對含噪的實正弦信號模型，可定義預測誤差函數e(n)為

(8)

(9)

1.3 基于信號自相關的MC算法

1.3.1 基于時域自相關的MC算法

由于正弦信號的自相關函數仍為與原信號同頻率的正弦信號，仍然具有線性預測特性。對照式(7)，自相關函數顯然也滿足：

Rτ+2≈2Rτ+1cosω0-Rτ

(10)

定義誤差函數：

eτ=Rτ+Rτ+2-2cosω0Rτ+1

(11)

的和取最小值，即可以求得:

(12)

式(12)用自相關序列{Rn1,Rn1+1,…,Rn2,Rn2+1,Rn2+2}來估計信號的頻率，其中n1、n2+2分別表示所需自相關系數的起始和截止下標。由式(6)可知，高斯白噪聲的自相關函數的期望只有在零點位置存在，其余位置均為0，也就是說序號小的自相關系數受噪聲的影響會大于序號大的自相關系數，但是從自相關的定義式(3)來看，當信號的長度N一定時，隨著自相關系數序號的增大，所用的信號變少，求和數變少，又會影響自相關本身定義的準確性，因此n1取值也不能太大。當自相關序列直接從時域，即從式(3)來求解時，該算法稱為時域自相關的MC算法。

1.3.2 基于頻域全頻帶自相關的MC算法

根據Weiner-Khintchine定理，信號的功率譜和其自相關函數服從一對傅里葉變換關系，自相關函數可以看成功率譜的傅里葉反變換，即：

(13)

2 基于頻域窄帶自相關的MC算法

2.1 算法思想

(14)

(15)

式(15)中：Δ為窄帶的前后寬度，表示選取譜線的范圍。式(15)和式(14)相比，限制了窄帶外的噪聲對頻率估計的影響，同時帶來計算量的減少。將式(15)代入式(12)，即可得到本文算法的頻率估計表達式：

(16)

此算法只用了窄帶里面的譜線，因此稱為頻域窄帶自相關的MC算法。

2.2 算法流程及計算量分析

該算法的實現流程如下：

該算法只用一次FFT，由于用的譜線少，式(15)的求和計算量少于IFFT的計算量，式(12)的計算量是一樣的，因此可以看出計算量比頻域全頻帶自相關的MC算法少。而式(3)計算自相關的計算量遠大于FFT，故該算法與時域自相關的MC算法相比，有明顯的計算優勢。

3 算法仿真及性能比較

3.1 不同補零倍數的仿真

本文提出的頻域窄帶自相關的MC算法是基于頻域進行分析，得到一個幅度譜譜峰，再通過自相關在時域的特點來進行頻率估計的，因此可以補零增加信號的長度來得到更為準確的譜峰，使最后的頻率估計更為精確。圖1是本文算法補p-1倍零的仿真結果，以及頻域全頻帶自相關的MC算法在p=1時候的仿真。注意p=1，表示補p-1=0倍零，即表示沒有補零。可以看出，在沒有補零(p=1)的時候，在低信噪比下，本文算法由于用窄帶功率譜表示自相關，去除了帶外的噪聲，性能優于頻域全頻帶自相關的MC算法，但是當噪聲影響不大，信噪比提高之后，頻域全頻帶自相關的MC算法由于用了全頻帶的信息，使得自相關系數更為準確，因而均方誤差低于本文算法。一旦補零后，本文算法的性能就得以提高，SNR≥-7 dB下均方誤差就達到CRB，且在隨著補零倍數的增多，性能越來越好。

圖1 本文算法不同p在不同信噪比下均方誤差的對比Fig.1 Mse vs SNR for different p of the proposed algorithm

3.2 自相關系數不同起始的仿真

從式(16)可知，不同的n1、n2對頻率估計的結果有影響，下面通過仿真來說明。圖2描述本文算法在n1=2、不同n2時的均方誤差比較。另外信號ω0=0.3π，SNR=20 dB，M=4N,Δ=3。為了保證自相關估計的估計精度，這里的n1=2，而n2的范圍是n1+1≤n2≤N-2。從圖2中可以看出本文算法在n2比較小的時候，即和n1相差不大時，所用的自相關系數比較少，性能比頻域(全頻帶)自相關的MC算法差，但是隨著n2的增大，本文算法的均方誤差一直下降，當n2接近200時，其性能超過頻域(全頻帶)自相關的MC算法，達到CRB，再增大n2，下降趨勢變得比較平坦。

圖2 本文算法在不同n2時均方誤差的對比，SNR=20 dB，n1=2Fig.2 Mse vs n2 with SNR=20 dB，n1=2 for the proposed algorithm

圖3 不同n1時算法均方誤差的對比，SNR=20 dB，n2=254Fig.3 Mse vs n1 with SNR=20 dB，n2=254 for the proposed algorithm

圖4列舉了(n1,n2)={(2,256)，(2,120)，(40,120)，(40,256)，(120,256)}這幾種情況下本文算法(p=4)的性能。 可以看出在SNR接近-7 dB時，這幾種情況下的均方誤差都能達到CRB，而頻域(全頻帶)自相關的MC算法(p=4)只有在SNR≥18 dB后才有比較好的性能。然而隨著信噪比的增大，自相關去噪的效果在減弱，反而是本身的自相關估計占據更為重要的作用，而序號比較小的比序號大的估計更為精確，同時自相關系數越多算法估計性能越好，如(n1,n2)=(2,256)；而自相關系數用的比較少的時候，算法估計性能變差，特別是用高序號自相關系數比較多，低序號自相關系數比較少時，性能變差，如(n1,n2)={(40,120)，(40,256)，(120,256)}比(n1,n2)=(2,120)差。

圖4 不同n1,n2，本文算法在不同信噪比下的均方誤差對比Fig.4 Mse vs SNR for different n1,n2 for the proposed algorithm

3.3 不同算法的性能比較

這里給出本文頻域窄帶自相關的MC算法與基于時域的原信號的MC算法、時域自相關的MC算法、和基于頻域的頻域全頻帶自相關的MC算法、文獻[15]提出的頻域三點內插算法進行比較。這里自相關的系數下標選用n1=2，n2=N-2，為了粗估計譜峰的統一，選用M=4N,Δ=3。

圖5給出了不同信噪比情況下的比較，其中ω0=0.3π。可以看出本文算法在SNR≥-7 dB下均方誤差達到CRB，估計器的性能遠遠超過基于時域的原信號的MC算法、時域自相關的MC算法；在SNR<10 dB時，性能優于基于頻域的頻域全頻帶自相關的MC算法，體現出低信噪比下去除帶外噪聲帶來的好處；由于頻域三點內插算法是基于原信號頻譜的，而基于頻域的頻域全頻帶自相關的MC算法和本文算法都補取了3倍信號長度的零，因此可以看到本文算法一直都優于它，基于頻域的頻域全頻帶自相關的MC算法在SNR>5 dB后性能也超過它。

圖6給出了SNR=10 dB時不同頻率處的性能比較，也可以明顯看到本文估計器的均方誤差遠遠低于基于時域的原信號的MC算法、時域自相關的MC算法和基于頻域的三點內插算法，在ω0∈[0.15π,0.85π]都達到CRB；在ω0∈[0.4π,0.6π]，本文算法和頻域全頻帶自相關的MC算法有相比擬的性能，而在其他頻率處都超過其性能。

圖5 不同信噪比下不同算法的均方誤差對比Fig.5 Mse vs SNR for different estimator

圖6 不同頻率ω0下不同算法的均方誤差對比，SNR=10 dBFig.6 Mse vs ω0 for different estimator with SNR=10 dB

最后，圖7給出了信號長度N對估計器性能的影響，同樣ω0=0.3π，SNR=10 dB。可以看到隨著N的增加，本文算法、時域自相關的MC算法、和基于頻域的三點內插算法和頻域全頻帶MC算法的性能得到了很大的提升，而MC估計器的性能趨于平緩。在N<150時，本文估計器的均方誤差比頻域全頻帶MC估計器低一點，且隨著N的增大，這種優勢也越來越明顯，前者可以達到CRB，而后者逐漸遠離。在不同長度下，本文算法都優于與之比較的其他算法。

圖7 不同信號長度N下不同算法的均方誤差對比，SNR=10 dBFig.7 Mse vs N for different estimator with SNR=10 dB

4 結論

結合時域和頻域自相關的特點，給出了一種基于窄帶自相關的實信號頻率估計算法，該算法通過在時域上避開選擇受噪聲影響比較大的低序號自相關系數，在頻域上利用信號的窄帶功率譜來計算自相關，去除帶外的噪聲，這樣來提高估計器的估計性能。仿真結果也表明該算法性能優于傳統的時域算法：MC算法和基于時域自相關的MC算法，以及傳統的頻域算法：頻域全頻帶自相關的MC算法和頻域三點內插算法，在-7 dB信噪比時就能逼近CRB界。特別是本文算法和基于時域自相關的MC算法、頻域全頻帶自相關的MC算法相比，都是基于自相關的算法，但是在性能上和計算量上卻有比較明顯的優勢。