基于FFT并二次修正的Rife頻率估計算法

謝 勝,陳 航,于 平

(1.西北工業(yè)大學航海學院,陜西西安 710072;2.91388部隊,廣東湛江 524022)

0 引言

在魚雷攻擊潛艇過程中,確定魚雷引信的炸點位置非常重要,這對于目標毀傷效果有直接的影響。工作區(qū)域主要在近場區(qū)的魚雷主動聲引信,回波的體目標效應明顯,引信接收到的信號為多個散射點回波信號的組合,回波信號的頻譜不再是單一的頻率,而是占據(jù)一定譜帶的、離散且隨時間變化的頻譜[1-2]。為高精度地提取潛艇目標回波各亮點多普勒頻率信息以便完成后續(xù)運動參數(shù)的估計,國內(nèi)外學者提出了一些快速高精度頻率估計的方法[3-4]。

基于DFT的譜分析法采用快速算法(FFT),運算速度快,特別適合實時信號處理,但由于利用FFT進行譜分析只能對有限多個樣本進行處理,不可避免地存在由于時域截斷產(chǎn)生的能量泄漏,使譜峰值變小、精度降低,即經(jīng)FFT得到的離散頻譜其幅值、相位和頻率都可能產(chǎn)生較大的誤差,于是人們又提出了基于FFT的各種插值算法。目前常用的基于FFT插值算法有:KAY法[5],雙線幅度 Rife法[6],結(jié)合幅度與相位的雙線Quinn法[7]等。Rife法和Quinn法由于算法簡單、計算量小,在工程實際中得到廣泛應用。文獻[8]分析了Rife法和Quinn法的頻率估計方差,指出:當信號頻率在接近DFT量化頻率中心區(qū)域時,Rife法和Quinn法標準差接近理論頻率估計方差下限;而當信號頻率接近DFT量化頻率且信噪比較低時,兩者估計誤差較大。對比Quinn法和Rife法的估計方差可知:在同等條件下Quinn法比 Rife法具有更小的估計方差。文獻[9]分析Rife法的性能,提出了一種修正 Rife(MRife)算法,通過對信號進行頻移,解決了當信號頻率位于量化頻率點附近時估計精度降低的問題,在適當?shù)男旁氡?0 d B以上)時,M-Rife算法在整個頻段上頻率估計均方誤差十分接近克拉美羅限CRLB;但當信噪比較低且當信號頻率位于DFT量化頻率附近時,M-Rife算法可能出現(xiàn)頻移方向出錯而導致估計誤差增大,需要進行二次修正。

針對M-Rife算法在低信噪比條件下估計性能下降的問題,本文提出一種基于FFT并二次修正的Rife頻率估計算法。

1 M-Rife算法

加性高斯白噪聲污染的信號表示為:

式中:a,f c,φ0分別為振幅、頻率和初相;Δt為采樣間隔;N為樣本數(shù);w(n)為實部和虛部相互獨立的、方差為2σ2的零均值復高斯白噪聲。

Rife算法公式可表示如下:

令Rife算法插值項為:

式中,|X(k0)|為X(n)的FFT最大譜線的幅值,|X(k0+r)|為FFT次大譜線的幅值,k0為FFT最大譜線位置。當|X(k0+1)|≤|(X(k0-1)|,r=-1;反之r=1,T=N?Δt。當信號頻率位于FFT量化頻率附近時,Rife算法在確定次大譜線時易受噪聲干擾發(fā)生插值方向錯誤,估計的誤差將可能大于DFT算法。

M-Rife算法建立在Rife算法基礎(chǔ)上,利用了Rife算法當信號頻率位于FFT兩相鄰量化頻率中心區(qū)域估計精度高的特點。其實現(xiàn)原理是:首先在對信號樣本作FFT基礎(chǔ)上用Rife算法進行頻率估計,如果估計出的頻率位于FFT兩相鄰量化頻率中心區(qū)域即|A|∈(ε,λ),其中0≤ε≤λ≤1/2,則以此估計值作為最終的頻率估計結(jié)果,否則對原始信號進行頻移p個量化單位以使信號頻率位于FFT兩相鄰量化頻率中心區(qū)域,再用Rife算法估計頻率。M-Rife算法在對信號進行頻移之前根據(jù)式(2)中r取值來確定頻移方向,若r=-1,向左頻移,反之向右頻移。也就是說M-Rife算法根據(jù)次大譜線位置來確定頻移方向。當信噪比較低且當信號頻率位于DFT量化頻率附近時,可能出現(xiàn)頻移方向出錯而導致估計誤差增大。M-Rife算法流程圖如圖1所示。

圖1 M-Rife算法流程圖Fig.1 flow chart of M-Rife algorithm

2 基于二次修正的Rife算法

針對M-Rife算法在低信噪比條件下估計性能下降的問題,提出一種基于FFT并二次修正的Rife頻率估計算法。由于利用FFT進行譜分析只能對有限多個樣本進行處理,不可避免由于時域截斷產(chǎn)生的能量泄漏導致FFT得到的離散頻譜其幅值、相位和頻率都可能產(chǎn)生較大的誤差。因此,改進算法首先通過計算機仿真方法搜索特定頻段的最佳修正因子m opt以對FFT最大譜線位置號進行初步修正?；贔FT估計的信號頻率可表示為 f=(k0+δ)?Δf,式中δ∈[-0.5,0.5],k0為FFT最大譜線位置號,Δf為FFT的頻率分辨率。假設(shè)被估計頻率在特定的頻段為均勻分布,在該頻段內(nèi)均勻選取N個頻率點,按照式(1)定義的輸入信號作FFT,根據(jù)公式f=(k0+δ)?Δf估計頻率,每個頻率點做1 000次蒙特卡羅仿真。以估計頻率的均方根誤差最小為準則求出各頻率點的最佳修正因子δi,然后將N個最佳修正因子δi進行平均處理得到該頻段內(nèi)最佳修正

量m opt:

在求得mopt的基礎(chǔ)上校正的Rife法公式可表示為:

前面分析已指出,M-Rife算法在Rife法基礎(chǔ)上根據(jù)次大譜線位置來確定頻移方向。當信噪比較低且當信號頻率位于DFT量化頻率附近時,可能出現(xiàn)頻移方向出錯而導致估計誤差增大。與M-Rife算法不同,改進算法在對原始樣本進行頻移之前,采用FFT最大譜線左右兩側(cè)譜線的相位信息即Quinn法來確定頻移方向,并且對原信號樣本作頻移后用校正的Rife法公式估計頻率時以Quinn法估計結(jié)果來確定插值方向。

Quinn算法利用FFT主瓣內(nèi)的次大譜線與最大譜線系數(shù)復數(shù)值之比的實部進行頻率估計。Quinn算法原理描述如下:設(shè)FFT最大譜線位置號k0,則k0-1,k0+1分別位于最大幅度值的兩側(cè),且其中一個為次大值。定義:

于是Quinn算法的頻率插值項可表示為:

改進算法的流程圖如圖2所示。

圖2 改進算法流程圖Fig.2 Flow chart of the improved algorithm

式中,p(p＞0)為頻移量,對新樣本作FFT,由式(3)求出 A′,然后將 A′減去p 得到A=A′-p,最后將k0+m opt和A代入校正的Rife法式(5)得到最終的頻率估計值。此時式(5)中r取值由Quinn法估計結(jié)果來確定,這與M-Rife算法是不同的。

其原理是:若估計的頻率位于FFT兩相鄰量化頻率中心區(qū)域即|A|∈(ε,λ),則用m opt對k0進行修正,得到k0+mopt,再與A一起代入校正的Rife法公式(5)估計頻率,此估計值作為本次實測樣本最終的頻率估計值(此時公式r取值與M-Rife算法同);否則若估計的頻率不位于FFT兩相鄰量化頻率中心區(qū)域,則根據(jù) Quinn 法公式(6)、(7)計算 ^δ1,^δ2 值并作判斷 :若＞0,＞0,令 r=1,向右頻移 ;否則令r=-1,向左頻移。頻移得到新樣本為:

3 仿真分析

在Matlab中對上面介紹的改進算法作了Monte Caro仿真。按照式(1)定義的輸入信號,相關(guān)參數(shù)設(shè)置如下:輸入信噪比SNR=a2/2σ2,信號初始相位φ0取[0,2π]區(qū)間均勻分布的隨機數(shù),噪聲均值為0、方差σ2=0.5。FFT量化頻率的中心區(qū)域為[0.3,0.5],頻移量 p取 0.3。對每個頻率點作1 000次Monte Caro仿真。

3.1 不同信號頻率的估計誤差

圖3為魚雷攻擊目標示意圖,魚雷與目標多普勒頻移:f d=2f c V/C,其中V為魚雷與目標連線方向相對速度,C為海水的聲速,f c為魚雷聲引信載頻。取脫靶量d=2 m,魚雷與目標在水平方向初始距離為X0=30 m,魚雷速度V TS=25 m/s,C=1 500 m/s,信噪比SNR=3 dB,聲引信載波頻率f 0=350 kHz,采樣頻率為 f s=1.4 MHz,信號樣本長度N=256,則FFT的頻率分辨率為Δf=1 376.2 Hz。假定點目標固定不動,則多普勒頻移范圍為[-15 000 Hz,15 000 Hz]。

圖3 魚雷攻擊目標示意圖Fig.3 Torpedo attacks taget

圖4 為用改進算法估計的魚雷與目標多普勒頻移圖,其中圖4(b)、(c)分別為對圖4(a)左右兩半部分曲線進行局部放大的效果圖。從圖4可看出,改進算法在整個頻段范圍內(nèi)頻率估計均值都非常接近真實值,估計多普勒頻移曲線與真實多普勒頻移曲線非常吻合。

圖4 改進算法估計多普勒頻移曲線圖Fig.4 The Doppler frequency shif t curveestimated by the improved algorithm

表1給出了魚雷與目標水平距離分別為6 m、5 m、…、0 m,信噪比為3 dB,樣本長度 N分別為256、512改進算法的估計均方根誤差。從表1可知,當N=256時改進算法估計均方根誤差不到FFT頻率分辨率Δf=1 376.2 Hz的1%,估計精度很高。當N=512時,估計均方根誤差最大值不大于36 Hz?？梢姼倪M算法在炸點附近的頻率估計精度能夠滿足反潛魚雷主動聲引信提取目標回波各亮點多普勒頻率的應用要求。

表1 SNR=3 d B,N=256、512,X=6 m、5 m、…、0 m條件下改進算法估計均方根誤差(RMSE)Tab.1 RMSE of the improved algorithm for SNR=3 dB,N=256,512,X=6 m,5 m,…,0 m

3.2 與M-Rife算法性能比較

圖5給出了與文獻[9]相同仿真參數(shù)條件下改進算法與M-Rife法的估計均方根誤差。圖5(a)—(c)細實線分別取自文獻[9]表2—4中當信噪比分別6 d B,0 d B,-3 dB時11個頻率點估計均方根誤差數(shù)據(jù)作圖而得到,粗實線表示改進算法的均方根誤差。從圖5(a)、(b)可看出,當信噪比SNR在0 dB以上時,改進算法和M-Rife法的均方根誤差均十分接近CRLB,兩者估計性能相當;而當SNR=-3 d B時,改進算法的估計穩(wěn)定性要明顯好于MRife算法。這是由于在低信噪比時且當信號頻率接近于FFT量化頻率時M-Rife法易因頻移方向錯誤導致估計誤差偏大,而改進算法利用了相位信息來判斷頻移方向,降低了頻移方向的錯誤概率,估計誤差得到有效控制。顯然,在低信噪比條件下,改進算法估計性能要優(yōu)于M-Rife法。

圖5 當SNR=6 d B,0 d B,-3 d B時,改進算法和M-Rife算法的均方根誤差Fig.5 RMSE of the improved algorithm and M-Rife for SNR=6 d B,0 d B,-3 d B

3.3 不同信噪比條件下估計誤差

為檢驗改進算法在不同信噪比條件下的估計性能,選取一個接近FFT量化頻率的頻率點 f 1=f 0+0.05Δf進行仿真,取 f 0=50 MHz,采樣頻率為fs=200 MHz,N=256,則FFT頻率分辨率Δf=781.2 k Hz。信噪比區(qū)間為[-8 dB,8 dB],間隔為1 d B,在每個SNR點作1 000次蒙特卡羅仿真,估計頻率方差下限為:

圖6給出了改進算法在不同信噪比條件下的均方根方差與CRLB比較。

從圖6可看出,盡管被估計信號頻率非常接近FFT量化頻率,但改進算法在很寬信噪比范圍內(nèi)均方根方差都十分接近CRLB。例如當SNR=-8 d B時,改進算法的均方根誤差約為CRLB的1.14倍。當SNR=-3 d B時,改進算法的均方根誤差約為CRLB的1.09倍;而在同等條件下M-Rife法其均方根誤差達到CRLB的1.4倍。隨著信噪比提高,改進算法的均方根誤差越來越接近CRLB。可見,改進算法比M-Rife法具有更寬信噪比門限。綜上所述,改進算法整體性能要優(yōu)于M-Rife法。

圖6 改進算法在不同信噪比條件均方根誤差Fig.6 RMSE of the improved algorithm in different SNR

4 結(jié)論

針對M-Rife算法在低信噪比時估計性能下降問題,本文提出了一種基于FFT并二次修正的Rife頻率估計算法。該算法首先利用修正因子對FFT最大譜線進行修正,然后用校正的Rife法公式估計頻率。若被估計信號頻率接近于FFT量化頻率,則對信號進行頻移以使信號頻率位于兩相鄰量化頻率的中心區(qū)域,再用校正的Rife法公式估計頻率。頻移方向由Quinn法來確定。計算機仿真結(jié)果表明:改進算法不僅在整個被估計的頻段內(nèi)具有較高估計精度,均方根誤差接近于CRLB,而且具有低信噪比門限,整體性能要優(yōu)于M-Rife法。改進算法計算量不大,便于硬件實現(xiàn),具有較高工程應用價值,能夠滿足反潛魚雷主動聲引信提取目標回波多普勒頻率的應用要求。

[1]Hao L,Paknys R.Comparison of near-field scattering for finite and infinitiearrays of parallel conducting strips TM incidence[J].IEEE Transactions on Antennas and Propagations,2005,53(11):3 735-3 740.

[2]Yu D,Fu D.A novel method of plannar near-field scattering measurements[C]//IEEE Internatinal Symposium on Microwave,Antenna,Propagation and EMC Technologies for Wireless Communications.USA:IEEE,2005:483-486.

[3]丁康,謝明.離散頻譜三點卷積幅值修正法的誤差分析[J].振動工程學報,1996,9(1):92-98.

[4]段虎明,秦樹人,李寧.離散頻譜的修正方法綜述[J].振動與沖擊,2007,26(11):138-145.DUAN Huming,QIN Shuren,LI Ning.Review of correction methods for discrete spectrum[J].Journal of Vibration and Shock,2007,26(11):138-145.

[5]Steven Kay.A Fast and Accurate Single Frequency Estimator[J].IEEE Transactions on Acoustics Speech and Signal Proeessing,1989,37(12):56-62.

[6]Rife D C,Vincent G A.Use of the Discrete Fourier Transform in the Measurement of Frequencies and Levels of Tones[J].Bell Sys Teeh J,1970,4(9):197-228.

[7]Quinn B G.Estimating frequency by interpolation using Fourier coefficients[J].IEEE Trans-SP,1994,42(5):1 264-1 268.

[8]齊國清.幾種基于FFT的頻率估計方法精度分析[J].振動工程學報,2006,19(1):86-92.QI Guoqing.Accuracy analysis and comparison of some FFT-based frequency estimators[J].Journal of Vibration Engineering,2006,19(1):86-92.

[9]鄧振淼,劉渝,王志忠.正弦波頻率估計的修正Rife算法[J].數(shù)據(jù)采集與處理,2006,21(4):473-477.DENG Zhenmiao,LIU Yu,WANG Zhizhong.Modified rife algorithm for frequency estimation of sinusoid wave[J].Journal of Data Acquisition&Processing,2006,21(4):473-477.