黃 超，索繼東，于 亮

(1.大連海事大學(xué)信息科學(xué)技術(shù)學(xué)院，遼寧大連116026;2.大連理工大學(xué)城市學(xué)院電子與自動(dòng)化學(xué)院，遼寧 大連116600;3.大連理工大學(xué) 軟件學(xué)院，遼寧 大連116600)

1 引言

對(duì)淹沒(méi)于噪聲中的正弦信號(hào)進(jìn)行頻率估計(jì)，無(wú)論在理論中還是在實(shí)際應(yīng)用中都具有非常重要的研究?jī)r(jià)值，國(guó)內(nèi)外眾多學(xué)者對(duì)此做了大量的研究，主要從頻域和時(shí)域兩個(gè)角度進(jìn)行分析，如基于FFT的頻率估計(jì)算法[1－2]。本文從時(shí)域的角度對(duì)頻率估計(jì)算法進(jìn)行研究。傳統(tǒng)基于時(shí)域的自相關(guān)方法有:V.Pisarenko利用少量的自相關(guān)函數(shù)進(jìn)行頻率估計(jì)[3]，此類(lèi)方法計(jì)算簡(jiǎn)單，但其性能不高;K.Lui等人在在文獻(xiàn)[3]的基礎(chǔ)上加以改進(jìn)，并定義新的自相關(guān)函數(shù)[4]，但由于選取的自相關(guān)函數(shù)系數(shù)較低，估計(jì)性能受到限制;K.Lui等人在文獻(xiàn)[4]的基礎(chǔ)上加以修正，采用較高序號(hào)的自相關(guān)系數(shù)[5]，使得估計(jì)性能得到一定的提升，但是在信噪比偏低時(shí)，估計(jì)性能不夠理想;為了充分挖掘自相關(guān)函數(shù)包含的頻率信息，Rim，H.C.So，Yan Cao等人盡可能地利用多個(gè)自相關(guān)函數(shù)進(jìn)行頻率估計(jì)[6－8]，使得估計(jì)性能得到很大提升，但同時(shí)算法的計(jì)算量也增大;為了進(jìn)一步提升頻率估計(jì)性能，K.Lui提出了兩步自相關(guān)算法[9]，Yan Cao等人提出了基于擴(kuò)展自相關(guān)的頻率估計(jì)算法[10]，在低信噪比估計(jì)性能得到進(jìn)一步改善，在中高信噪比時(shí)，頻率估計(jì)方差接近克拉美羅下界(CRLB)[11]。然而，在現(xiàn)有很多算法中，估計(jì)性能的提升都是利用更多的自相關(guān)系數(shù)或者多步自相關(guān)函數(shù)，這樣必然會(huì)帶來(lái)算法計(jì)算量的增大，即存在頻率估計(jì)精度與算法復(fù)雜度的矛盾問(wèn)題。

為了解決上述問(wèn)題，本文同樣從自相關(guān)函數(shù)的角度，充分利用自相關(guān)函數(shù)包含的頻率信息，推導(dǎo)了一種新的自相關(guān)函數(shù)相位的頻率估計(jì)式，并且針對(duì)頻率估計(jì)范圍與頻率估計(jì)精度之間的矛盾問(wèn)題，提出了一種消除相位模糊的方法。計(jì)算機(jī)仿真表明，在信噪比較高時(shí)，估計(jì)方差接近克拉美羅下界(CRLB)，與TSA算法相比，在估計(jì)性能相同條件下，本文算法的計(jì)算量大大降低，具有很好的應(yīng)用價(jià)值。

2 正弦信號(hào)頻率估計(jì)算法原理

設(shè)混有高斯加性白噪聲的單頻正弦信號(hào)表示式為

其中，a、ω0、θ分別為信號(hào)的幅度、頻率、相位，η(t)為均值為零、方差為σ2的加性高斯白噪聲。對(duì)其在觀察時(shí)間0≤t≤T內(nèi)進(jìn)行采樣N個(gè)樣本值，于是得到離散序列為

對(duì)上式定義其自相關(guān)函數(shù)為

當(dāng)N足夠大時(shí)，由文獻(xiàn)[10]可知，式(3)可寫(xiě)為

Pisarenko[3]算法利用低階自相關(guān) r1、r2進(jìn)行頻率估計(jì)，如下式所示:

Pisarenko算法雖然計(jì)算量小，但是頻率估計(jì)方差和CRB相差較大。為了進(jìn)一步減小頻率估計(jì)偏差，K.Lui等人在PHD算法的基礎(chǔ)上加以改進(jìn)，提出了K.Lui算法[4]。其算法僅對(duì)Pisarenko算法的r1、r2進(jìn)行重新定義，改進(jìn)了r1、r2的估計(jì)性能，如下式所示:

于是，式(5)可改寫(xiě)為

為了進(jìn)一步提高載波頻率估計(jì)性能，TSA算法[9]利用更多的自相關(guān)函數(shù)進(jìn)行頻率估計(jì):

其中:

從上式可以看出，TSA算法經(jīng)過(guò)自相關(guān)函數(shù)的多次組合，充分挖掘自相關(guān)函數(shù)包含的頻率信息進(jìn)行頻率估值，理論分析及仿真結(jié)果表明，SNR≥0 dB時(shí)，頻率估計(jì)方差接近CRLB，且頻率估計(jì)范圍寬;但是其算法的計(jì)算量較大(與O(N3)成正比)，不利于實(shí)時(shí)信號(hào)的處理。

3 本文提出的新的頻率估計(jì)式

目前，基于自相關(guān)函數(shù)的頻率估計(jì)算法中，很多都是通過(guò)利用多個(gè)自相關(guān)系數(shù)，使得頻率估計(jì)的性能得到提升，如 Yan 算法[8，10]、TSA 算法[9]等，但同時(shí)會(huì)帶來(lái)計(jì)算量的增加，不利于實(shí)時(shí)通信的信號(hào)處理。針對(duì)頻率估計(jì)精度與算法復(fù)雜度的矛盾問(wèn)題，本文從自相關(guān)函數(shù)相位的角度，推導(dǎo)了一種新的頻率估計(jì)式，較好地解決了上述矛盾問(wèn)題。

由正弦信號(hào)的三角函數(shù)特性可知

于是，可得下式:

為了進(jìn)一步提高頻率估計(jì)性能，對(duì)上式進(jìn)行展開(kāi):

其中，p≥1，p≤q≤N－2m－1。于是，可得出本文新的頻率估計(jì)式:

式中，ξm為整數(shù)。

為了研究參數(shù)(p，q，m)對(duì)頻率估值的影響，對(duì)上式進(jìn)行了仿真。在仿真中，當(dāng)(m，p，q·m)在(0，π)范圍內(nèi)時(shí)，此時(shí) ξm=0，令信號(hào)幅度 a=1，θ=0，ω0=0.01π，SNR=10 dB，N=200，p=1(蒙特卡洛仿真100次)。圖1給出了m=1時(shí)不同q值時(shí)的頻率估計(jì)方差，從圖1可以看出，隨著q值增大，估計(jì)性能不斷提升，當(dāng)q為50左右時(shí)，估計(jì)性能逐漸接近最佳值。圖2給出了q=1時(shí)不同m取值的估計(jì)方差，從圖2可以看出，當(dāng)m為60左右時(shí)，頻率估計(jì)性能最佳。

由上式可得頻率估計(jì)值為

圖2 不同m值的估計(jì)性能Fig.2 Mean square error versus m

以上仿真中，假設(shè)頻率范圍為(0，π/m)。從圖1和圖2可以得出，為了提高頻率估計(jì)精度，(p，q，m)參數(shù)的取值很關(guān)鍵，不難得出(p=1，q=50，m=60)的性能會(huì)好于(p=1，q=50，m=1)及(p=1，q=1，m=60)的性能;而m＞1時(shí)會(huì)縮小頻率估計(jì)范圍，為了擴(kuò)大估計(jì)范圍，即m，p，q∈(0，π)，此時(shí)ξm有2m+1種可能的值，即存在相位模糊問(wèn)題，為了準(zhǔn)確地找出 ξm值，消除相位模糊，可采取以下措施(令 m=1，2，…，ε，):

(1)m=1，此時(shí)不存在相位模糊問(wèn)題，即ξ1=0，可求出頻率估計(jì)值1，p，q;(2)m=2，求出

時(shí)對(duì)應(yīng)的 ξ2，再根據(jù)式(13)求出頻率估計(jì)值(|·|min為取最小值);

(3)以此類(lèi)推;

(4)m=ε，求出

時(shí)對(duì)應(yīng)的ξε，再根據(jù)式(13)求出頻率估計(jì)值ε，p，q。

綜上可知，q值的增大會(huì)帶來(lái)頻率估計(jì)精度的提升，但同時(shí)會(huì)導(dǎo)致計(jì)算量過(guò)大，因此，在實(shí)際應(yīng)用中，應(yīng)綜合考慮估計(jì)性能和計(jì)算量之間矛盾的問(wèn)題。

4 計(jì)算量及性能分析

圖1 不同q值的估計(jì)性能Fig.1 Mean square error versus q

4.1 計(jì)算機(jī)仿真及性能分析

本節(jié)通過(guò)計(jì)算機(jī)MATLAB仿真實(shí)驗(yàn)(蒙特卡洛仿真100次)來(lái)驗(yàn)證本文提出算法的估計(jì)性能。圖3給出SNR=10 dB時(shí)不同頻率處的性能比較，本文頻率估計(jì)算法的性能(m=50，p=1，q=70)在 ω0∈[0.15π，0.85π]達(dá)到 CRLB，頻率估計(jì)范圍與 TSA算法[9]相當(dāng)，都寬于 K.Lui算法[5]的頻率估計(jì)法范圍。圖4給出了ω0=0.3π、N=200時(shí)不同信噪比下各種算法的性能比較，可以看出，在SNR＞0 dB時(shí)，隨著信噪比增大，本文算法(m=50，p=1，q=50，m=50，p=1，q=70 及m=50，p=1，q=80)與TSA 算法的估計(jì)性能相當(dāng)，都接近于CRLB，且性能都好于K.Lui算法[5];在 SNR＜0 dB 時(shí)，隨著 SNR 的減小，估計(jì)性能偏離CRLB的程度從大到小依次為K.Lui算法[5]、本文算法(m=50，p=1，q=50)、本文算法(m=50，p=1，q=70)、本文算法(m=50，p=1，q=80)和 TSA 算法[9]。

圖3 不同頻率處的性能比較Fig.3 Mean square error versus ω0

圖4 不同信噪比下的性能比較Fig.4 Mean square error versus SNR

綜上分析可得，在信噪比較高時(shí)，本文算法中q值可以選擇較小(q值越小，計(jì)算量越小)，如q=50，其估計(jì)性能也接近CRLB;在信噪比較低時(shí)，本文算法中q值可以選擇較大，如q=80，其估計(jì)性能偏離CRLB最小。

4.2 算法復(fù)雜度比較

從表1可以看出，本文算法的計(jì)算量(與 O(N2)成正比)高于 K.Lui算法[5]，低于 TSA算法[9]的計(jì)算量(與 O(N3)成正比)。由式(8)、式(9)、式(13)可以看出，K.Lui算法[5]只利用幾個(gè)自相關(guān)系數(shù)進(jìn)行頻率估計(jì)，本文算法利用多個(gè)自相關(guān)系數(shù)且通過(guò)遞推的方式求出頻率值，而TSA算法[9]是在一步自相關(guān)的基礎(chǔ)上再作自相關(guān)運(yùn)算，其計(jì)算量最大。從圖4的仿真結(jié)果來(lái)看，本文算法(m=50，p=1，q=80)與TSA算法的頻率估計(jì)性能非常接近，但本文算法的計(jì)算量約為T(mén)SA算法[9]的2/5，非常有利于信號(hào)的實(shí)時(shí)處理。

表1 算法計(jì)算量比較Table1 Comparison of arithmetic operations

5 結(jié)論

充分利用自相關(guān)函數(shù)包含的頻率信息進(jìn)行頻率估計(jì)，可以使估計(jì)性能得到提升。本文針對(duì)加性高斯白噪聲的正弦信號(hào)，提出了基于加窗自相關(guān)函數(shù)相位噪聲的頻率估計(jì)算法。文中推導(dǎo)了自相關(guān)函數(shù)相位噪聲的頻率估計(jì)式，并針對(duì)頻率估計(jì)范圍與頻率估計(jì)精度之間的矛盾問(wèn)題，提出了一種消除相位模糊的方法。仿真結(jié)果表明，在信噪比較高時(shí)，估計(jì)方差接近克拉美羅下界(CRLB)，且在估計(jì)性能相同條件下，與TSA算法相比，本文算法的計(jì)算量大大降低，在工程上具有很好的應(yīng)用價(jià)值。

然而，在現(xiàn)有的眾多算法中，低信噪比時(shí)的估計(jì)性能還是不夠理想，如何更好地挖掘自相關(guān)函數(shù)包含的頻率信息，降低頻率估計(jì)的信噪比閥值，將是我們進(jìn)一步研究的內(nèi)容。

[1]王宏偉，趙國(guó)慶，齊非林.一種實(shí)時(shí)精確的正弦波頻率估計(jì)算法[J].數(shù)據(jù)采集與處理，2009，24(2):208－211.WANG Hong－wei，ZHAO Guo－qing，QI Fei－lin.Real－Time and Accurate Single Frequency Estimation Approach[J].Journal of Data Acquisition ＆ Processing，2009，24(2):208－211.(in Chinese)

[2]路偉濤，楊文革，洪家財(cái)，等.一種新的正弦信號(hào)頻率和初相估計(jì)方法[J].電訊技術(shù)，2012，52(9):1459－1464.LU Wei－tao，YANG Wen－ge，HONG Jia－cai，et al.A Novel Method for Frequency and Initial Phase Estimation of Single－ tone Signals[J].Telecommunication Engineering，2012，52(9):1459－1464.(in Chinese)

[3]Pisarenko V.The retrieval of harmonicsby linear prediction[J].Geophysical Journal of the Royal Astronomical Society，1973，33(5):347－366.

[4]Lui K，So H C.Improved Variant of Pisarenko Harmonic Decomposition for Single Sinusoidal Frequency Estimation[J].IEICE Transactions on Fundamentals of Electronics，Communications and Computer Sciences，2007，90(11):2604－2607.

[5]Lui K，So H C.Modified Pisarenko Harmonic Decomposition for Single－tone frequency estimation[J].IEEE Transactions on Signal Processing，2008，56(7):3351－3356.

[6]Elasmi－Ksibi R，Besbes H，L pez－Valcarce R，et al.Frequency estimation of real－valued single－tone in colored noise using multiple autocorrelation lags[J].IEEE Transactions on Signal Processing，2010，90(7):2303－2307.

[7]So H C，Chan K W.Reformulation of Pisarenko harmonic decomposition method for single－tone frequency estimation[J].IEEE Transactions on Signal Processing，2004，52(4):1128－1135.

[8]Cao Yan，Wei Gang，Chen Fang－Jiong.An exact analysis of Modified Covariance frequency estimation algorithm based on correlation of single－tone[J].IEEE Transactions on Signal Processing，2012，92(4):2785－2790.

[9]Lui K，So H C.Two－stage autocorrelation approach for accurate single sinusoidal frequency estimation[J].IEEE Transactions on Signal Processing，2008，88(7):1852－1857.

[10]Cao Yan，Wei Gang，Chen Fang－Jiong.A ClosedformExpanded Autocorrelation Method for Frequency Estimation of A Sinusoid[J].IEEE Transactions on Signal Processing，2012，92(4):1852－1857.

[11]Funga H W，Alex C，Kotb K H Li，et al.ParameterEstimation of A Real Single Tone from Short Data Records[J].IEEE Transactions on Signal Processing，2004，84(3):601－617.