基于全相位濾波的密集成分的頻率估計

楊俊東， 余 江， 黃 銘， 張 果， 葛孚華

(1. 云南大學(xué) 信息學(xué)院，昆明 650504； 2. 昆明理工大學(xué) 信息工程與自動化學(xué)院，昆明 650500)

頻率測量是振動工程中最常見的問題之一。許多物理量的測量(如大氣動力學(xué)中的風速測量[1]、陣列信號處理中的波達方向估計[2]、振動分析中的轉(zhuǎn)速測量[3]等)都可轉(zhuǎn)化為頻率估計問題；頻率估計也是頻譜分析理論中最基本的問題，而其中基于FFT(Fast Fourier Transform)的經(jīng)典譜分析方法，相比于各種現(xiàn)代譜估計方法(如AR模型法[4]、多重信號分類法[5]、旋轉(zhuǎn)不變參數(shù)估計法[6]等)，具有其所不能比擬的運算效率高、耗費內(nèi)存小的優(yōu)勢，故在儀表設(shè)計、電力諧波分析、旋轉(zhuǎn)機械故障診斷等工程領(lǐng)域中得到非常廣泛的應(yīng)用，基于FFT的頻率估計器一直是學(xué)術(shù)研究的熱點問題。

影響FFT頻率估計精度主要有兩方面因素：頻譜泄漏和柵欄效應(yīng)。頻譜泄漏可通過加窗來改善，對于單復(fù)指數(shù)成分，理想的傅里葉譜是沖激函數(shù)(即對應(yīng)單根譜線)，但因存在譜泄漏，其直接FFT的可能結(jié)果是產(chǎn)生五根以上的大能量值譜線，加窗FFT可起到抑制泄漏的效果，表現(xiàn)在其主要能量通常集中在三根譜線上；而柵欄效應(yīng)則可以通過頻譜校正的途徑來解決，具體說來，就是借助內(nèi)插措施，對峰值譜附近的譜線做進一步修正和細化來提升估計精度，設(shè)計不同的內(nèi)插算法就產(chǎn)生多種頻率校正器(如Quinn校正器、能量重心校正器、Macleod校正器、Candan校正器、比值校正器、AM譜估計器、Tsui估計器等)。

文獻[7]指出，頻譜校正誤差來源于三個方面：系統(tǒng)固有誤差、譜間干擾誤差和噪聲干擾誤差。需指出，經(jīng)過文獻[8-15]在內(nèi)插算法設(shè)計方面的不斷完善，其系統(tǒng)誤差做得越來越小，精度越來越高，可以說幾乎不存在提升的空間(如AM估計器的頻率估計均方誤差僅為克拉美-羅限的1.014 7倍，再降低幾乎很困難)。然而，需指出的是，文獻[8-15]的頻譜校正算法都是針對單頻復(fù)指數(shù)信號而推導(dǎo)出的，即沒有考慮存在譜間干擾誤差的情況，故只適用于單頻信號或者間隔較遠的多頻信號情況，對于頻率間隔較近的密集譜信號，由于存在嚴重的譜間干擾，其頻率估計精度會急劇降低。因而如何提升密集譜校正精度是亟待解決的問題。

密集譜校正具有很高的工程意義，因為密集譜實際上是和短樣本情況緊密聯(lián)系在一起的：在采樣速率fs不變的情況下，其FFT頻率分辨單位Δf=fs/N(N為樣本長度)，因而當工程上無法采集到足夠樣本時(如地震波測量和爆破勘探中，N只能取小數(shù)值)， 則Δf變大， 這時會導(dǎo)致兩個頻率成分的間距|f1-f2|/Δf減小，從而短樣本情況趨于變?yōu)轭l率估計難度更高的密集譜情況。

本文提出基于全相位濾波的密集譜頻率估計算法，該算法以FFT為中心的三根譜線的相位譜值的線性偏離程度作為密集譜的判據(jù)，借助全相位點通濾波器依次去除所關(guān)注的成分的帶外干擾，進而結(jié)合比值頻譜校正器，可以大大減小系統(tǒng)誤差，提升頻率估計精度。具有較高的實用價值。

1 密集譜頻率估計問題

1.1 單頻信號的直接FFT譜和加窗FFT譜特性

若以采樣速率fs對頻率為f1的單頻復(fù)指數(shù)信號采集N個樣點， 其離散序列可表示為

(1)

令FFT的頻率分辨率Δω=2π/N， 相應(yīng)地，x1(n)對應(yīng)的數(shù)字角頻率ω1可表示為

(2)

令矩形窗傅里葉幅度譜函數(shù)為RN(ω)， 不難證明，x1(n)的直接FFT譜X1(k)表達式為

(3)

令某對稱窗序列w(n)的傅里葉幅度譜W(ω)， 則對x1(n)做加窗FFT得到的離散譜表達式為

(4)

令信號頻率f1=5.2 Hz，幅值a1=1， 初相θ1=60°，采樣速率fs=32 Hz，樣本長度N=32，則Δf=1 Hz，其直接FFT和加漢寧窗FFT的振幅譜，如圖1所示。

圖1 直接FFT和加窗FFT的振幅譜圖Fig.1 Amplitude spectra for direct FFT and windowed FFT

從圖1可知，加窗FFT相比于直接FFT， 其譜泄漏范圍大為減少(由于W(ω)的旁瓣衰減高于RN(ω)的緣故)。進而對于密集譜估計場合，為盡量減小譜間干擾誤差，只能用加窗FFT。

1.2 密集成分的直接FFT譜和加窗FFT譜特性

在單頻成分基礎(chǔ)上，假定信號包含兩個余弦成分

x1(n)=a1cos(ω1n+θ1), 0≤n≤N-1

(5)

x2(n)=a2cos(ω2n+θ2), 0≤n≤N-1

(6)

從而復(fù)合信號表示為

x(n)=x1(n)+x2(n)

(7)

令x1(n)的參數(shù)與圖1情況相同， 對x2(n)做如下參數(shù)設(shè)置：f2=6.9 Hz，幅值a2=1，初相θ2=45°， 則單頻信號x1(n)，x2(n)和復(fù)合信號x(n)的直接FFT和加漢寧窗FFT的振幅譜和相位譜，如圖2所示。

圖2 密集譜成分及復(fù)合信號的振幅譜和相位譜Fig.2 Amplitude spectra & phase spectra for individual dense components and their composite signal

從圖2(a)的振幅譜可看出， |X1(k)|主峰位于k=5處， 在左右相鄰的兩根譜線(即k=4，k=6)處產(chǎn)生大幅值旁譜泄漏；而|X2(k)|主峰位于k=7處，在相鄰的兩根譜線(即k=6，k=8)處產(chǎn)生大幅值旁譜泄漏。很明顯，由于兩頻率成分相隔太近，兩者在共同位置k=6處產(chǎn)生了譜重疊，這必然引起譜間干擾而降低譜校正精度。因而，降低譜間干擾是提升密集譜校正精度的關(guān)鍵。

2 基于全相位濾波的密集譜頻率估計算法

2.1 基于相位譜線性度偏離檢測的密集譜識別判據(jù)

顯然，密集譜校正的前提是能識別出密集成分。這單單依據(jù)振幅譜是不夠的：以圖2(a)為例， 從其復(fù)合振幅譜|X(k)|中仍能判斷出k=5和k=6的兩處譜峰，但卻不能推斷出兩者的相互影響程度。因而，除了振幅譜檢測外，引入相位譜檢測是必要的，故本文提出基于相位譜線性度偏離檢測的密集譜識別方法，原理如下。

從式(4)可知，加窗后的單頻復(fù)指數(shù)信號的FFT相位譜的表達式為

(8)

具體說來，圖2(b)給出了兩個單頻信號的相位譜φ1(k)，φ2(k)和復(fù)合信號的相位譜φ(k)(為盡量保持數(shù)值連續(xù)性，對這三個相位譜值都做了解纏繞處理)，可看出：單頻相位譜φ1(k)在k=4,k=5,k=6三根譜線上(即在虛框內(nèi))可連成直線；類似地，單頻相位譜φ2(k)則在k=6,k=7,k=8三根譜線上(即在虛框內(nèi))也可連成直線。

具體說來，先計算如下兩個差分相位

(9)

再求取其均值，即

(10)

進而定義如下相位線性偏離因子

(11)

2.2 基于全相位濾波的譜成分分離

為提升密集譜校正精度，不妨采用濾波方法分離各個密集成分，這對濾波器性能提出如下要求：①為分離出密集成分，濾波器帶寬必須足夠窄；②為精確定位所關(guān)注頻率的位置，濾波器的中心頻點必須可以精確控制；③為不改變各成分之間的相對時延，濾波器需具有線性相位特性；④濾波器設(shè)計應(yīng)足夠簡單。

為滿足以上需求，本文選取Huang等提出的全相位窄帶線性相位FIR濾波器來實現(xiàn)密集譜分離。該濾波器的優(yōu)勢在于：一個N階的全相位濾波器蘊含了N個子濾波器的相互補償?shù)倪^程，因而傳輸特性好；而且全相位濾波器系數(shù)可推出解析公式[16-19]，經(jīng)過Huang等的改進，對于其窄帶濾波情況，其濾波器設(shè)計非常簡單，僅依據(jù)下式即可確定所有濾波器系數(shù)

(12)

(13)

其中，式(13)的歸一化因子C為

(14)

仍以上述例子說明對復(fù)合密集譜信號的濾波過程：圖3(a)為前述的復(fù)合信號加窗FFT譜，為保留第1個譜成分，濾除k=6處的譜間干擾，將ωp=5Δω代入式(12)可得到全相位濾波器系數(shù)g(n)， 其對應(yīng)的傳輸曲線|G(jw)|如圖3(b)所示，從圖3(b)可知， |G(jw)|不僅將中心頻點精確定位在ωp=5Δω處，而且還可以抑制k=6處的譜間干擾。圖3(c)給出了濾波前復(fù)合信號的32個樣點波形, 圖3(d)給出了濾后的的94個樣點波形。

圖3 全相位濾波器傳輸曲線及濾波波形Fig.3 Transfer curve of the all-phase filter and the waveforms

需指出：由于輸入信號長度為N, 式(12)的全相位濾波器的系數(shù)長度為2N-1， 而FIR濾波輸入輸出之間為線性卷積關(guān)系，故圖3(d)的濾波輸出為N+(2N-1)-1=3N-2=94個樣點，文獻[20]指出，全相位濾波器存在(N-1)個樣點的群延時，故應(yīng)取從n=N=32～n=2N-1=63的N個樣本段做為有效的濾波輸出y(n)(如虛框所界定)。對比圖3(d)和圖3(c)可知，濾波輸出的有效樣本段y(n)相比于輸入x(n)波形更為平穩(wěn)，有利于提升后續(xù)的譜校正精度。進而對y(n)做漢寧窗FFT，即可得到頻譜|Y(k)|。

類似地，令ωp=7Δω，重復(fù)以上過程，也可分離出第2個頻率成分。

2.3 比值頻譜校正法

針對密集譜情況，顯然應(yīng)選取使用譜線數(shù)少的校正算法。而比值校正法因僅需兩根譜線(峰值譜和幅度次大的旁譜線)故成為首選。其譜校正過程如下

(1) 算出加窗FFT峰值譜和次高譜的幅度比值

(15)

u=(2-v)/(1+v)

(16)

(17)

2.4 算法總結(jié)

基于以上分析，將提出的算法流程總結(jié)如圖4所示。

圖4 密集譜頻率估計算法流程Fig.4 Dataflow of frequency estimator for dense components

3 仿真試驗

3.1 無噪試驗

分別用原有比值校正法、Tsui校正法(因其單頻情況做頻率估計噪聲估計方差接近克拉美羅限而選作對比ξ=0.5)和本文提出的基于全相位濾波的密集譜校正法(設(shè)置閾值)對圖2的兩個密集成分的信號進行頻率校正，其校正結(jié)果如表1所示。

表1 三種校正法的頻率估計結(jié)果

從表1可知，由于引入了相位線性偏離度判決和全相位濾波，本文提出的校正器比原比值校正器和Tsui估計器的精度有所提升，第1個成分的校正誤差分別從0.255 9 Hz，0.152 1 Hz降低到-0.105 7 Hz，第2個成分的校正誤差分別從-0.336 9 Hz，-0.176 2 Hz降低到-0.105 7 Hz。

這驗證了以上提及的結(jié)論：對于現(xiàn)有的校正器，當未非密集譜情況時，其單頻成分的頻率校正均具有很高的精度；當所關(guān)注的成分周圍存在干擾譜時，因沒有考慮去除干擾的濾波措施，原有各校正器的內(nèi)插公式準確度會降低。

3.2 加噪測量試驗

(18)

表2給出了三種信噪比情況下的RMSE(Root Mean Square Error)統(tǒng)計結(jié)果。

表2 三種校正法的均方根誤差

從表2的RMSE統(tǒng)計結(jié)果，可得出如下結(jié)論：

(1) 總體上，無論在哪種信噪比情況下，本文方法在估計密集成分時(即f1=5.1 Hz與f2=6.8 Hz為密集成分，f3=15.1 Hz與f4=17.2 Hz也為密集成分)，其均方根誤差都低于原比值法和Tsui方法，這驗證了本文提出的相位線性偏離度因子在辨識密集成分和全相位濾波在分離密集成分方面是有效的。

(2) 進一步觀察，可發(fā)現(xiàn)SNR越高，本文方法均方誤差相比于原比值法和Tsui方法的均方根誤差的相對比例就越小，即精度改善越明顯。這是因為，如前所述，校正器的誤差包括系統(tǒng)誤差、譜間干擾誤差和噪聲誤差三部分，由于都是基于比值法校正，兩者的系統(tǒng)誤差可近似認為相等，而由于本文引入了全相位濾波，故譜間干擾誤差相比于傳統(tǒng)比值法要減小。故隨著信噪比升高，在同樣的低噪聲干擾情況下，這就越加凸顯出全相位濾波在降低譜間干擾誤差方面對精度改善的貢獻。

(3) 另外可看出，對于離密集譜較遠的成分f5=22.3 Hz，本文方法的RMSE值反而高于原比值法和Tsui方法。這也很容易得到解釋：一方面，對于這種非密集譜，譜間干擾其實可以忽略，如前所述，在沒有譜間干擾情況下，對于現(xiàn)有的頻譜校正器來說，幾乎不存在精度提升的空間；另一方面，由于噪聲干擾是隨機的，總體說來，相位譜相比于振幅譜而言，對噪聲的敏感度會更高些，因而對于非密集譜情況，理論上根本不需要引入相位線性偏離度判據(jù)和全相位窄帶濾波的場合，反而出現(xiàn)小概率的判據(jù)出錯情況而增大總體的RMSE誤差。

4 結(jié) 論

本文提出基于全相位窄帶濾波的密集譜頻率估計法，該方法通過所提出的譜峰附近的相位線性偏離度因子，可以辨識密集譜是否存在，進而借助全相位窄帶濾波，可實現(xiàn)各密集頻率成分的有效分離，從而降低譜間干擾誤差，再結(jié)合現(xiàn)有的頻譜校正法(如比值校正法)，即可提升密集譜的頻率估計精度。

鑒于頻率估計在雷達、光學(xué)、聲納、電力、地震勘探等領(lǐng)域的應(yīng)用廣泛性，本文提出的基于全相位濾波的頻率估計器具有較廣闊的應(yīng)用前景。