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

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

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

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

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

現(xiàn)結(jié)合時(shí)域和頻域自相關(guān)的特點(diǎn)，提出一種基于窄帶自相關(guān)的實(shí)信號(hào)頻率估計(jì)算法。該算法在頻域進(jìn)行譜峰搜索后，利用信號(hào)的窄帶功率譜來計(jì)算自相關(guān)，以消除窄帶外的噪聲對(duì)頻率估計(jì)的影響，然后再用簡單的時(shí)域自相關(guān)的MC算法來得到頻率估計(jì)的閉式解。最后給出與時(shí)域自相關(guān)MC算法，以及頻域全頻帶的自相關(guān)MC算法和頻域三點(diǎn)內(nèi)插算法[15]相比較的仿真結(jié)果。

1 基于自相關(guān)的實(shí)信號(hào)頻率估計(jì)

1.1 實(shí)信號(hào)的自相關(guān)函數(shù)

假設(shè)原實(shí)正弦信號(hào)為

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

(1)

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

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

(2)

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

(3)

式中：τ為自相關(guān)函數(shù)的下標(biāo)，τ=0，1,…,N-1，有噪聲情況下，實(shí)信號(hào)的自相關(guān)函數(shù)為

(4)

當(dāng)N足夠大時(shí)，自相關(guān)函數(shù)的期望為

(5)

式(5)中：δτ,0為克羅內(nèi)克(Kronecker)函數(shù)：

(6)

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

1.2 基于原信號(hào)的MC算法

由于正弦信號(hào)s(n)具有線性預(yù)測(cè)特性：s(n)=2cosω0s(n-1)-s(n-2)，也可以時(shí)延后寫成：

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

(7)

基于這種線性預(yù)測(cè)，對(duì)含噪的實(shí)正弦信號(hào)模型，可定義預(yù)測(cè)誤差函數(shù)e(n)為

(8)

(9)

1.3 基于信號(hào)自相關(guān)的MC算法

1.3.1 基于時(shí)域自相關(guān)的MC算法

由于正弦信號(hào)的自相關(guān)函數(shù)仍為與原信號(hào)同頻率的正弦信號(hào)，仍然具有線性預(yù)測(cè)特性。對(duì)照式(7)，自相關(guān)函數(shù)顯然也滿足：

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

(10)

定義誤差函數(shù)：

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

(11)

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

(12)

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

1.3.2 基于頻域全頻帶自相關(guān)的MC算法

根據(jù)Weiner-Khintchine定理，信號(hào)的功率譜和其自相關(guān)函數(shù)服從一對(duì)傅里葉變換關(guān)系，自相關(guān)函數(shù)可以看成功率譜的傅里葉反變換，即：

(13)

2 基于頻域窄帶自相關(guān)的MC算法

2.1 算法思想

(14)

(15)

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

(16)

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

2.2 算法流程及計(jì)算量分析

該算法的實(shí)現(xiàn)流程如下：

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

3 算法仿真及性能比較

3.1 不同補(bǔ)零倍數(shù)的仿真

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

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

3.2 自相關(guān)系數(shù)不同起始的仿真

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

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

圖3 不同n1時(shí)算法均方誤差的對(duì)比，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時(shí)，這幾種情況下的均方誤差都能達(dá)到CRB，而頻域(全頻帶)自相關(guān)的MC算法(p=4)只有在SNR≥18 dB后才有比較好的性能。然而隨著信噪比的增大，自相關(guān)去噪的效果在減弱，反而是本身的自相關(guān)估計(jì)占據(jù)更為重要的作用，而序號(hào)比較小的比序號(hào)大的估計(jì)更為精確，同時(shí)自相關(guān)系數(shù)越多算法估計(jì)性能越好，如(n1,n2)=(2,256)；而自相關(guān)系數(shù)用的比較少的時(shí)候，算法估計(jì)性能變差，特別是用高序號(hào)自相關(guān)系數(shù)比較多，低序號(hào)自相關(guān)系數(shù)比較少時(shí)，性能變差，如(n1,n2)={(40,120)，(40,256)，(120,256)}比(n1,n2)=(2,120)差。

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

3.3 不同算法的性能比較

這里給出本文頻域窄帶自相關(guān)的MC算法與基于時(shí)域的原信號(hào)的MC算法、時(shí)域自相關(guān)的MC算法、和基于頻域的頻域全頻帶自相關(guān)的MC算法、文獻(xiàn)[15]提出的頻域三點(diǎn)內(nèi)插算法進(jìn)行比較。這里自相關(guān)的系數(shù)下標(biāo)選用n1=2，n2=N-2，為了粗估計(jì)譜峰的統(tǒng)一，選用M=4N,Δ=3。

圖5給出了不同信噪比情況下的比較，其中ω0=0.3π?？梢钥闯霰疚乃惴ㄔ赟NR≥-7 dB下均方誤差達(dá)到CRB，估計(jì)器的性能遠(yuǎn)遠(yuǎn)超過基于時(shí)域的原信號(hào)的MC算法、時(shí)域自相關(guān)的MC算法；在SNR<10 dB時(shí)，性能優(yōu)于基于頻域的頻域全頻帶自相關(guān)的MC算法，體現(xiàn)出低信噪比下去除帶外噪聲帶來的好處；由于頻域三點(diǎn)內(nèi)插算法是基于原信號(hào)頻譜的，而基于頻域的頻域全頻帶自相關(guān)的MC算法和本文算法都補(bǔ)取了3倍信號(hào)長度的零，因此可以看到本文算法一直都優(yōu)于它，基于頻域的頻域全頻帶自相關(guān)的MC算法在SNR>5 dB后性能也超過它。

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

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

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

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

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

4 結(jié)論

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