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

謝 勝,陳 航,于 平

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

0 引言

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

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

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

1 M-Rife算法

加性高斯白噪聲污染的信號(hào)表示為:

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

Rife算法公式可表示如下:

令Rife算法插值項(xiàng)為:

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

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

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

2 基于二次修正的Rife算法

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

量m opt:

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

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

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

于是Quinn算法的頻率插值項(xiàng)可表示為:

改進(jìn)算法的流程圖如圖2所示。

圖2 改進(jìn)算法流程圖Fig.2 Flow chart of the improved algorithm

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

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

3 仿真分析

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

3.1 不同信號(hào)頻率的估計(jì)誤差

圖3為魚雷攻擊目標(biāo)示意圖,魚雷與目標(biāo)多普勒頻移:f d=2f c V/C,其中V為魚雷與目標(biāo)連線方向相對(duì)速度,C為海水的聲速,f c為魚雷聲引信載頻。取脫靶量d=2 m,魚雷與目標(biāo)在水平方向初始距離為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,信號(hào)樣本長(zhǎng)度N=256,則FFT的頻率分辨率為Δf=1 376.2 Hz。假定點(diǎn)目標(biāo)固定不動(dòng),則多普勒頻移范圍為[-15 000 Hz,15 000 Hz]。

圖3 魚雷攻擊目標(biāo)示意圖Fig.3 Torpedo attacks taget

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

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

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

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

圖5 當(dāng)SNR=6 d B,0 d B,-3 d B時(shí),改進(jìn)算法和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 不同信噪比條件下估計(jì)誤差

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

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

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

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

4 結(jié)論

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

[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]丁康,謝明.離散頻譜三點(diǎn)卷積幅值修正法的誤差分析[J].振動(dòng)工程學(xué)報(bào),1996,9(1):92-98.

[4]段虎明,秦樹人,李寧.離散頻譜的修正方法綜述[J].振動(dòng)與沖擊,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]齊國(guó)清.幾種基于FFT的頻率估計(jì)方法精度分析[J].振動(dòng)工程學(xué)報(bào),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]鄧振淼,劉渝,王志忠.正弦波頻率估計(jì)的修正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.