復白噪聲中復正弦波頻率估計方法研究

王 樂 蘇小敏 杜 林 李春化

(西安電子工程研究所 西安 710100)

1 引言

對于隨機信號，由于其無始無終、能量無限，對其作傅里葉變換不收斂，因而不能像確定性信號那樣確定此類信號的頻譜。雖然此類信號能量無限，但是其功率未必無限，因此對于隨機信號而言，常用功率譜描述其頻率特性。功率譜估計分為兩大類:參數(shù)化方法和非參數(shù)化方法。非參數(shù)化方法依賴于傳統(tǒng)的傅里葉變換法，其包括 BT法、周期圖法等［1］，參數(shù)化方法包括MUSIC算法、AR模型算法、最大熵等估計算法。利用MUSIC算法、AR模型、最大熵等現(xiàn)代譜估計算法可以對正弦信號頻率進行精確的估計［2］，但是這些算法復雜、計算量大，難以實時處理從而影響其應用。

對于任意形式的確定性信號x(n)，對其進行傅里葉分解，都可以分解為直流與許多余弦(或正弦)分量線性組合［1］。然而在基于離散傅里葉變換(DFT)的經(jīng)典譜估計中，由時域截斷引起的頻譜泄露是不可避免的，能量泄露會引起較大的估計誤差［3］，這時我們需要對頻譜進行校正以期獲得精確的譜峰位置。本文提出一種功率加權求平均的頻譜校正算法，并通過仿真實驗在基于FFT分析基礎上與單譜線法、Rife法進行對比，然后討論了觀察時間對FFT分析的影響以及信號采樣頻率的選擇，最后針對短持續(xù)時間信號，給出頻譜細化方法(Zoom-FFT)，并通過實驗與FFT方法進行對比。

2 頻譜校正算法

設正弦信號為:

式中，A、f0、φ0代表信號的振幅、頻率和初相。對其采樣，得到離散序列:

其中，Δt=1/fs代表采樣間隔，fs為采樣頻率，M為采樣點數(shù)。對(2)式進行離散傅里葉變換(DFT)得:

直接計算DFT，其計算量是與變換區(qū)間長度N的平方成正比的。當N較大時，計算量非常的大，因此直接用DFT算法進行譜分析和信號的實時處理是不切實際的。工程實際中常用快速傅里葉變換(FFT)來計算DFT，其離散頻譜譜峰值對應的頻率作為信號頻率的初步估計，即:

其中，k表示頻譜峰值對應的位置。由于其僅僅利用了峰值譜線，又被稱為單譜線幅度法。

我們知道，對一個信號作DFT，前提要求其必須是有限長度的，因此必須對信號作加窗處理，窗函數(shù)的引入會導致頻譜泄漏。又因為作離散傅里葉變換相當于對整個頻率范圍(0～fs)內(nèi)作N點均勻采樣，在均勻采樣過程中我們并不一定能采得幅度峰值，這就導致了k取值的不準確性。該頻率估計方法分辨率較低，最大頻率誤差可以達到 ± fs/2N［2］。

對于基于基帶FFT法利用單根譜線估計頻率造成誤差過大問題，人們提出了諸多頻譜校正方法。其中較著名、效果較好的是雙線幅度法(Rife法)。Rife法利用最大譜線和次大譜線對單譜線方法進行修正，當信噪比較高時估計性能很好。但當信噪比較低或者待估計頻率靠近最大譜線時，會導致次大譜線判斷的失誤，造成估計誤差比僅僅用單譜線粗略估計誤差還要大，此時則要利用插值細化技術對Rife法作進一步改進，關于 Rife法詳見參考文獻［4～5］。

本文提出了一種利用最大譜線及其鄰近譜線，對其進行功率加權并取平均的方法來彌補單譜線估計的不足。方法示意圖見圖1。

圖1 功率加權求平均校正頻譜

圖1 中虛線部分表示頻譜包絡，其中頻譜峰值對應頻率為f0，假設在其左右各取兩根相鄰譜線，分別對應頻率值為 f-2、f-1、f1、f2，所取五根譜線如圖所示，對應幅值分別為 l-2、l-1、l0、l1、l2，則信號頻率估值為:

上述主要介紹了在整個頻域(0～fs)范圍內(nèi)，對序列信號進行頻域上一維搜索，找出頻譜峰值，用其所對應的頻率值作為原信號的頻率估計值。這種估計方法原理清晰、易于理解，但計算量較大。若要提高估計精度，延長觀測時間是一種行之有效的方法。采用MATLAB仿真工具進行驗證，假設采樣頻率fs=2048Hz，采樣間隔 Ts=1/fs，信號頻率 f0=125.65Hz，信號幅度A=2.0，噪聲設為正態(tài)分布隨機噪聲，方差為0.4，信噪比SNR=10dB，當觀測時間分別為 t1=0.1s、t2=0.5s、t3=1s、t4=1.5s、t5=3s，采用單譜線法、Rife法和功率加權求平均方法所估計頻率為，結果表 1 所示。

表1 頻譜校正對比實驗估計頻率結果

估計頻率值減去頻率真值，所得差值再取絕對值，得到各方法的估計誤差，如圖2所示。

圖2 校正頻譜后的估計誤差

由圖2可以看出:a.雙線幅度法(Rife法)并不一定優(yōu)于單譜線幅度法;

b.當增加觀察時間時，柵欄效應引起的頻譜泄漏影響減小，頻率估計誤差減小，估計頻率精度增高。

當觀察時間固定，采樣頻率提高時，采樣點數(shù)也相應增多，但這種情況下，頻率分辨率未增加，頻率估計效果是不會提高的。這是因為，雖然采樣頻率提高會帶來采樣點數(shù)N的增加，但同樣導致頻率軸上的折疊頻率(fs/2)增加，頻率分辨率fs/N未變。實驗驗證，假設采樣頻率 fs1=400Hz、fs2=800Hz、fs3=1200Hz、fs4=1600Hz、fs5=2000Hz，采樣間隔 Ti=1/fsi，信號頻率 f0=125.65Hz，信號幅度 A=2.0，噪聲設為正態(tài)分布隨機噪聲，方差為0.4，信噪比SNR=10dB，當采樣時間t=1s，采用上述單譜線法、Rife法和功率加權求平均方法所估計頻率為、，結果如表2所示。

表2 頻率遞增對比實驗估計頻率結果

估計頻率值減去頻率真值，所得差值再取絕對值得到各方法的估計誤差，如圖3所示。

由圖3可以看出，采樣頻率變化時頻率估計誤差變化不超過0.5Hz，這種些微變化可以認為是由隨機噪聲引起的，而與采樣頻率遞增無關。

然而通過上面仿真實驗，可以看出采樣頻率提高，無論哪種估計方法誤差變化趨勢是一樣的，這就提示我們頻率估計誤差跟采樣頻率之間并不是毫無關聯(lián)的。

在頻率估計中，采樣率的選擇實際上是有一定限制的，不能隨意更改采樣頻率的大小，否則會造成估計精度變差。換句話說，就是固定的采樣頻率只適合對一定頻率范圍內(nèi)的信號進行采樣，否則會給頻率估計帶來問題，以下進行仿真驗證。

圖3 頻率遞增后的估計誤差

假設采樣頻率fs=100MHz，采樣間隔Ts=1/fs，信號幅度A=1.0、噪聲Ψ(n)為正態(tài)分布隨機噪聲，方差為0.01，信噪比SNR=20dB，信號觀測時間t=1μs，信號頻率為f，則信號可以表示為:

通過仿真，ω取100個不同值所對應的頻率估計實驗結果如圖4所示，圖中橫坐標表示ω在0～π之間取值，縱坐標表示頻率估計誤差均方根，圖4(右)為去掉圖4(左)中始末兩點所得。

圖4 頻率估計誤差均方根隨信號頻率變化

由圖4可以看出:當信號角頻率ω位于0.2π～0.8π之時，即信號頻率時，頻率估計誤差均方根較小，估計效果較好，當信號頻率時，頻率估計均方根非常大，估計不穩(wěn)定，這驗證了采樣頻率的選取是有一定限制的結論，因此我們在作頻率估計時，涉及到采樣頻率時應該使2.5f＜fs＜10f，其中fs的上下限并不嚴格限定為10f和2.5f，僅代表fs在取值上注意脫離危險區(qū)域。

以上論述了幾種較簡單的頻率校正方法，通過仿真實驗證明了延長觀測時間有利于頻率估計結果的結論，并且討論采樣頻率與信號頻率之間的關系。

然而用增加觀察時間來提高頻率估計精度的方法并不是放之四海而皆準的，首先觀察時間的增加勢必帶來計算量的增大;其次在短持續(xù)時間信號中，觀測時間是有限的，不能通過增加觀測時間來提高頻率估計精度，這種情況下就需要改進算法，才能獲得理想的頻率估計精度。

3 頻譜細化算法

在許多實際應用中，我們并不是對整個頻譜感興趣，往往只關注其中的一個窄頻帶，并且需要對這一窄頻帶作局部放大以進行細致的觀察，為此我們應該在所關心的頻帶內(nèi)增加譜線，以增大譜線密度［6］。近年來，頻譜細化技術得到迅猛發(fā)展，常見的方法有基于復調制的Zoom-FFT法、線性調頻Z變換法(Chirp-Z變換，簡稱CZT)、Yip-Zoom法、相位補償細化法等。但從分析精度、計算效率、分辨率、靈活性等方面來看，Zoom-FFT法是一種非常有效的方法［7］，在工程中也得到了非常廣泛的應用。

Zoom-FFT法(以下簡稱ZFFT法)是提供高分辨率FFT分析的一種技術，它允許在給定的一組FFT譜線中選擇起始和終止頻率，在所得的窄頻帶上能以指定的、足夠高的采樣頻率去分析，可以成倍的提高頻率分辨率，這樣就可以對我們感興趣的頻率段進行細化。

在實際應用中，我們通常對信號做基帶FFT估計，得到信號的粗測頻率，然后以此粗測頻率作為窄頻率帶中心頻率，向左右延展確定起始頻率和終止頻率，起始頻率和終止頻率之間的窄頻帶作為細化分析頻帶，確定細化倍數(shù)，移頻再進行復FFT處理。Zoom-FFT處理過程流程見圖5所示。

圖5 Zoom-FFT處理框圖

算法具體步驟可以歸納為:

a.復調制移頻，即將需要細化的頻帶中心(FFT粗測所得)移至頻率軸原點，通常我們用數(shù)字下變頻來實現(xiàn);

b.低通濾波，我們需要濾出感興趣的頻段信號，加以分析。假設頻率細化倍數(shù)為d，則低通濾波器的截止頻率fc=fs/2d;

c.抽樣。信號被調制移頻和低通濾波之后，所得到的信號頻帶變窄，因而可以用較低的采樣頻率f's=fs/d進行重采樣，即對原采樣點每隔固定點抽樣一次;

d.復FFT處理和頻率回調。抽樣之后得到復序列信號，對其進行復FFT處理，此時分辨率提高了d倍，然后進行頻率調整，將譜線移至實際頻率處即可得到細化的頻譜。

仿真實驗假設采樣頻率fs=100MHz，采樣間隔Ts=1/fs，信號幅度A=2.0、噪聲為正態(tài)分布隨機噪聲，方差為0.4，信噪比SNR=10dB，信號采樣時間t=1μs，信號頻率 f=30.123MHz，頻帶中心 fa=30MHz，細化倍數(shù)為 d=10，Zoom-FFT分析帶寬為25MHz～35MHz，用基帶FFT和ZFFT法進行頻率估計，仿真結果見圖5和圖6所示。圖中橫坐標表示頻率軸，縱坐標表示幅值。用功率加權求平均的方法進行頻譜校正，做100次Monte-carlo實驗得到仿真結果:FFT方法估計結果=30.215MHz，同真值比較，誤差為Δ=0.092MHz;ZFFT方法頻率估計為，誤差為

由FFT與ZFFT頻率估計對比實驗可以看出，ZFFT方法頻率估計在信號主瓣內(nèi)采樣點數(shù)明顯多于FFT方法，這些意味著頻譜分辨率更高，頻率估計精度更大。

4 總結

頻譜校正算法可以很好的修正由獲取最大譜線位置不準確所帶來的頻率估計誤差。本文的主要工作是:首先，通過實驗對基于FFT分析的單譜線、雙譜線和功率加權等頻譜校正方法進行對比，得出功率加權算法優(yōu)于單譜線法和雙譜線法，這是因為功率加權法利用了更多的譜線信息;其次，通過實驗驗證了延長觀測時間利于頻率估計以及采樣頻率與信號頻率之間的關系;最后，針對短持續(xù)時間信號觀測時間有限的特點，介紹了一種行之有效的頻譜細化方法，并通過仿真實驗與FFT方法作對比，證明了ZFFT算法估計效果優(yōu)于FFT算法。

［1］陸光華，彭學愚，張林讓，毛用才.隨機信號處理［M］.西安:西安電子科技大學出版社，2003.

［2］王洪先，王玉軍，王宏偉.寬帶自相關接收機中的一種頻率估計方法［J］.火控雷達技術，2011，40(1):52-54.

［3］李天昀，葛臨東.基于ZFFT和Chirp-Z的頻譜細化分析中能量泄露的研究［J］.電子對抗技術，2003，18(5):11-15.

［4］STEVE KAY.A Fast and Accurate Single Frequency Estimator［J］.IEEE TRANS.ON ACOUSTICS，1989，37(12):1987-1989.

［5］王紀強，歐攀，張春熹，林志立.基于頻偏校正的頻率估計算法誤差分析［J］.北京航空航天大學學報，2010，36(7):849-852.

［6］王力，張冰，徐偉.基于MATLAB復調制ZOOMFFT算法的分析和實現(xiàn)［J］.艦船電子工程，2006，4:119-121.

［7］丁康，潘成灝，李巍華.ZFFT與Chirp-Z變換細化選帶的頻譜分析對比［J］.振動與沖擊，2006，25(6):9-12.