黃翔東， 黎鳴詩， 羅 蓬， 馬 欣

(1. 天津大學(xué) 電氣自動(dòng)化與信息工程學(xué)院，天津 300072； 2. 國網(wǎng)河北省電力公司電力科學(xué)研究院， 石家莊 050021)

頻率估計(jì)的應(yīng)用非常廣泛(如多普勒效應(yīng)檢測[1]，振動(dòng)學(xué)中的轉(zhuǎn)速測量[2]都可轉(zhuǎn)化為頻率估計(jì)問題)。相比于基于子空間分解(如多重信號(hào)分類法[3]、旋轉(zhuǎn)不變參數(shù)估計(jì)法[4]等)的頻率估計(jì)算法，F(xiàn)FT(Fast Fourier Transform)因其具有運(yùn)算效率高的優(yōu)勢，因而基于FFT的頻率估計(jì)器一直是學(xué)術(shù)研究的熱點(diǎn)問題。?

Stoica等[5]指出：給定有限長樣本，最大似然意義上的頻率估計(jì)值位于該序列的DTFT(Discrete Time Fourier Transform)的峰值位置，但是DTFT的理想譜峰需要對頻率做高分辨率掃描才可搜索得到, 顯然頻率掃描會(huì)增加估計(jì)器的計(jì)算復(fù)雜度。為降低復(fù)雜度, 通常的做法是對序列做FFT(而不是DTFT)，進(jìn)一步對FFT峰值譜附近的譜線采取內(nèi)插、迭代等措施，對FFT譜峰附近的譜值做進(jìn)一步修正和細(xì)化來提升頻率估計(jì)精度。從而不同的內(nèi)插算法就可以產(chǎn)生不同的頻率估計(jì)器(如Quinn估計(jì)器[6]、能量重心估計(jì)器[7]、Macleod估計(jì)器[8]、Jacobsen估計(jì)器[9]、Candan估計(jì)器[10]、相位差估計(jì)器[11-14]、比值估計(jì)器[15]、AM譜估計(jì)器[16]、細(xì)化譜估計(jì)器[17]、Tsui估計(jì)器[18]等)。

需指出，對于單頻信號(hào)的FFT頻率估計(jì)器，經(jīng)過不斷完善，其精度幾乎不存在提升的空間(如文獻(xiàn)[16]提出的AM估計(jì)器的頻率估計(jì)方差僅為克拉美-羅限的1.014 7倍，而文獻(xiàn)[18]提出的Tsui估計(jì)器的估計(jì)方差幾乎逼近克拉美-羅限[19](Cramer-Rao lower Bound，CRB)。但是當(dāng)對多頻信號(hào)做FFT譜分析時(shí)，由于FFT固有的譜泄漏，各頻率成分之間會(huì)嚴(yán)重的譜間干擾，導(dǎo)致Quinn等的估計(jì)器精度都會(huì)大幅度降低。

因而，提升多頻信號(hào)的頻率估計(jì)精度的關(guān)鍵在于：選擇比FFT更優(yōu)良的抑制譜泄漏性能的譜分析方法，基于此設(shè)計(jì)出具體的內(nèi)插措施。王兆華等[20-21]指出，全相位FFT(all-phase FFT，apFFT)相比于傳統(tǒng)FFT，具有更優(yōu)良的抑制譜泄漏性能，黃翔東等[22]還把離散的全相位FFT譜分析延伸到連續(xù)的全相位DTFT譜分析(all-phase DTFT, apDTFT)領(lǐng)域，從而拓寬了其理論內(nèi)涵。故本文用全相位FFT替代傳統(tǒng)FFT譜分析，通過分析信號(hào)頻偏值與譜峰附近apDTFT譜的內(nèi)在聯(lián)系，設(shè)計(jì)出一種迭代多頻內(nèi)插估計(jì)器。仿真實(shí)驗(yàn)表明：給定同樣數(shù)目的觀測樣本，本文估計(jì)器比Tsui估計(jì)器具有更高的估計(jì)精度。由于工程應(yīng)用中的多頻信號(hào)比單頻信號(hào)更為廣泛，故本文估計(jì)器具有較高的應(yīng)用價(jià)值。

1 全相位FFT和全相位DTFT基本原理

1.1 全相位FFT譜分析

王兆華等[23]提出的apFFT譜分析簡化流程，如圖1所示。

圖1 apFFT譜分析簡化流程(N=4)

圖1的apFFT譜分析分為兩個(gè)簡單步驟：

步驟1全相位數(shù)據(jù)處理，用長為2N-1的卷積窗wc(n)對輸入數(shù)據(jù)x(n)加權(quán)，然后將間隔為N的數(shù)據(jù)兩兩疊加而得到y(tǒng)(0),y(1),…,y(N-1)；

步驟2對y(0),y(1),…,y(N-1)做FFT即得全相位離散譜Y(k)。

圖1中的卷積窗wc(n)是由長為N的窗f(n)和翻轉(zhuǎn)的窗b(n)通過卷積得到，即

wc(n)=f(n)*b(-n),n∈[-N+1,N-1]

(1)

本文取f(n)=b(n)=RN(n)，其中RN是長為N的矩形窗。

令Δω=2π/N，ω0=(k*+δ)Δω,其中k*∈Z+,δ∈[-0.5,0.5)。黃翔東等[22]已證明，對于單頻復(fù)指數(shù)信號(hào)

x(n)=A·ej(ω0n+θ0),n∈[-N+1,N-1]

(2)

其歸一化的apFFT譜值為

k=0…,N-1

(3)

式(3)的平方項(xiàng)使得旁譜線Y(k),k≠k*,相對于峰值譜Y(k*)出現(xiàn)大幅度衰減，使得峰值譜更突出，因而apFFT具有比傳統(tǒng)FFT更優(yōu)良的抑制譜泄漏性能。

1.2 全相位DTFT譜分析

黃翔東等提出的apDTFT的簡化流程，如圖2所示。

圖2 apDTFT譜分析的簡化流程

圖2的apDTFT輸出為

(4)

式中:ω為頻率掃描變量。因而不同于離散的apFFT譜，apDTFT譜Y(jω)為連續(xù)譜。

黃翔東等已證明，對于式(2)中的單頻復(fù)指數(shù)信號(hào)，其apDTFT譜的理論值Y(jω)為

(5)

聯(lián)立式(3)，式(5)，有

Y(k)=Y(jω)|ω=kΔω,k=0,…,N-1

(6)

因而apFFT譜Y(k)可看作是在ω∈[0,2π)內(nèi)對apDTFT譜Y(jω)的等間隔離散采樣結(jié)果。兩者都是全相位譜分析理論的重要組成部分。

特別地，易證明在距離譜峰的兩個(gè)對稱頻點(diǎn)ωL=ω0-0.5Δω、ωR=ω0+0.5Δω處的apDTFT值同為

(7)

2 全相位譜分析的迭代多頻內(nèi)插估計(jì)器設(shè)計(jì)

2.1 apDTFT峰值譜附近的頻偏估計(jì)機(jī)理

ωL=kLΔω=(kc-0.5)Δω

(8)

ωR=kRΔω=(kc+0.5)Δω

(9)

把式(8)和式(9)代入式(5)，可推導(dǎo)出這兩個(gè)頻點(diǎn)處的apDTFT譜幅值分別為

(10)

(11)

聯(lián)立式(10)和式(11)可得左譜線與右譜線幅值之比為

(12)

因?yàn)棣摹蔥-0.5,0.5)，則式(12)等號(hào)右邊去掉絕對值符號(hào)后可表示為

(13)

當(dāng)N足夠大，根據(jù)無窮小關(guān)系：sin[(δ-0.5)π/N] ～(δ-0.5)π/N，sin[(δ+0.5)π/N]～(δ+0.5)π/N，則式(13)可近似表示為

由表5可知，安徽省區(qū)域旅游發(fā)展水平差異明顯．黃山市與池州市旅游區(qū)位熵得分較高，旅游產(chǎn)業(yè)發(fā)展水平在省內(nèi)處于領(lǐng)先地位．黃山市、池州市、安慶市、宣城市、合肥市、蕪湖市、六安市旅游區(qū)位熵均大于0.74，旅游業(yè)發(fā)展?fàn)顩r較好，可劃分為旅游發(fā)展的優(yōu)勢區(qū)域;其他城市均小于0.74，旅游業(yè)發(fā)展?fàn)顩r一般，可劃分為一般旅游區(qū)域．

(14)

利用上式可解得頻偏估計(jì)值

(15)

2.2 迭代估計(jì)過程

從圖3可看出，只在初始化階段做了1次apFFT譜分析，而且迭代過程中也只涉及兩個(gè)特殊頻點(diǎn)ωL、ωR的apDTFT值的更新，從而繞開了apDTFT連續(xù)譜的頻率掃描過程，故圖3流程復(fù)雜度低。

2.3 迭代趨勢分析

搜索出apFFT峰值譜位置后，在第1次迭代中，兩個(gè)頻點(diǎn)ωL、ωR上的apDTFT譜不外乎兩種情況：|Y(jωL)|≤|Y(jωR)|和|Y(jωR)|<|Y(jωL)|，其分布分別如圖4(a)、4(b)所示。

圖3 單頻迭代估計(jì)流程圖

(a) |Y(jωL)|<|Y(jωR)|

(b) |Y(jωL)|>|Y(jωR)|

理想情況下，由式可推知|Y(jωL)|=|Y(jωR)|成立。故對于圖3迭代流程，可采用兩者差異的比例作為迭代終止條件。經(jīng)驗(yàn)表明，將閾值取為ε∈(0.004，0.3)，只需2次迭代即可獲得高精度頻率估計(jì)。

2.4 多頻成分的頻率估計(jì)

需指出，與單頻成分不同的是，多頻成分的頻率估計(jì)誤差來源除了外加噪聲和估計(jì)器本身的系統(tǒng)誤差(如式(14)中等價(jià)無窮小引入的誤差)外，各成分的譜間干擾也是誤差源的主要方面。由于apFFT譜分析具有比傳統(tǒng)FFT譜分析更優(yōu)良的抑制譜泄漏性能，故譜間干擾更小，因而其多頻估計(jì)器精度高于文獻(xiàn)[6-18]中的估計(jì)器。以上結(jié)論將通過第3節(jié)中的實(shí)驗(yàn)來證明。

2.5 多頻信號(hào)重構(gòu)

(16)

(17)

式中,ang(·)表示取相角操作。

從而，包含多個(gè)頻率成分的實(shí)信號(hào)的重構(gòu)式為

n∈[-N+1,N-1]

(18)

3 實(shí) 驗(yàn)

例1令N=64,M=3,ω1=7.1Δω,ω2=9.35Δω,ω3=12.2Δω,對信號(hào)x(n)=ejω1n+3ejω2n+ejω3n,n=-N+1,…,N-1，分別進(jìn)行64階apFFT、apDTFT及其127點(diǎn)的FFT和DTFT，其譜圖分別如圖5(a)、5(b)所示。并且分別用本文的頻率估計(jì)器和文獻(xiàn)[18]的Tsui估計(jì)器(選該估計(jì)器的原因是：其誤差方差逼近CRB，幾乎是精度最高的FFT估計(jì)器)分別進(jìn)行頻率估計(jì)，表1給出兩個(gè)估計(jì)器的相對誤差，其算式為

(19)

從圖5明顯可看出，apFFT、apDTFT的譜泄漏遠(yuǎn)小于FFT、DTFT的譜泄漏(即除3個(gè)主峰外，前者的帶外起伏比后者小得多)，這意味著對于3個(gè)成分間的譜干擾，前者也比后者小得多。從表1數(shù)據(jù)也可看出，本文估計(jì)器的相對誤差也遠(yuǎn)小于Tsui估計(jì)器情況。

(a) apFFT、apDTFT譜圖

(b) FFT、DTFT譜圖

η1/%η2/%η3/%本文估計(jì)器0.120.010.06Tsui估計(jì)器0.330.030.36

例2調(diào)研含噪情況下，不同頻率間距的各成分的譜間干擾對頻率估計(jì)器的影響。不妨令N=64,M=2,輸入信號(hào)為x(n)=ejω1n+ejω2n+ξ(n)，其中ξ(n)為加性高斯復(fù)噪聲。固定ω1=5.25Δω，ω2分別取值為9.15Δω，8.15Δω，7.15Δω，6.15Δω。圖6～圖9給出不同信噪比(Signal to Noise Ratio, SNR)條件下，本文估計(jì)器、Tsui估計(jì)器的均方根誤差以及CRB的平方根曲線。

從圖6～圖9的RMSE曲線，可得出如下結(jié)論：

(1) 在低信噪比區(qū)域(SNR<10 dB)，本文估計(jì)器的SNR閾值普遍比Tsui估計(jì)器低，這意味著其抗噪魯棒性強(qiáng)于Tsui估計(jì)器。這是因?yàn)閷τ赥sui估計(jì)器，多頻信號(hào)較大的譜間干擾削弱了抵御強(qiáng)噪聲的能力。

(2) 在中、高信噪比區(qū)域(SNR>10 dB)，對于Tsui估計(jì)器，其RMSE呈現(xiàn)出平坦的形狀，與CRB曲線偏離非常遠(yuǎn)，而本文估計(jì)器僅在兩頻率成分間距較小時(shí)(即圖8、圖9的情況)，才會(huì)與CRB曲線產(chǎn)生較大的偏離。這說明，在中、高信噪比情況，影響估計(jì)器精度的主導(dǎo)因素，既不是外加噪聲，也不是系統(tǒng)誤差，而是各成分間的干擾，而全相位譜分析的譜間干擾遠(yuǎn)小于傳統(tǒng)譜分析情況，故其RMSE曲線更貼近于理論下限。

(3) 當(dāng)頻率間距過小時(shí)(低于1Δω，即對應(yīng)圖9的情況)，兩個(gè)估計(jì)器的RMSE都變得很大，這是因?yàn)闊o論全相位譜，還是傳統(tǒng)FFT譜，其各成分譜的主瓣都相互重疊了，故頻率估計(jì)失效。

(a) ω1=5.25Δω

(b) ω2=9.15Δω

(a) ω1=5.25Δω

(b) ω2=8.15Δω

(a) ω1=5.25Δω

(b) ω2=7.15Δω

(a) ω1=5.25Δω

(b) ω2=6.15Δω

γ=10lg(Px/Pe)

(20)

式中，Px為原信號(hào)x(n)的平均功率，Pe為e(n)的平均功率。表2給出了兩種估計(jì)器的重構(gòu)。

表2 機(jī)械故障信號(hào)的兩種估計(jì)器重構(gòu)信噪比

以測量面缺陷信號(hào)為例，圖10(a)、(b)分別給出了用全相位FFT的振幅譜|Y(k)|、傳統(tǒng)FFT的振幅譜|X(k)|，圖10(c)給出了原始波形、全相位重構(gòu)波形和Tsui重構(gòu)波形。

(a) apFFT振幅譜

(b) FFT振幅譜

(c) 原始波形與兩種重構(gòu)波形

可明顯看出，圖10(a)的振幅譜|Y(k)|的譜間干擾程度比圖10(b)的振幅譜|X(k)|的譜小很多，這保證了前者的各成分的頻率、幅值、相位參數(shù)估計(jì)精度高于后者，從而獲得更高的重構(gòu)信噪比(如表2所示，前者比后者高出2.7 dB)。該結(jié)論對于其他兩種故障信號(hào)的重構(gòu)也成立。

4 結(jié) 論

本文提出基于全相位譜分析的多頻內(nèi)插迭代頻率估計(jì)器，利用全相位譜分析方法有效地抑制譜泄漏，降低多頻成分間的干擾；再結(jié)合內(nèi)插迭代提高多頻信號(hào)的頻率估計(jì)精度。由于本文方法在有外界干擾的低信噪比環(huán)境仍能高精度地做頻率估計(jì)，且仿真實(shí)驗(yàn)和實(shí)際信號(hào)重構(gòu)實(shí)驗(yàn)均驗(yàn)證了本文方法具有更高的參數(shù)估計(jì)精度，故在多頻信號(hào)估計(jì)領(lǐng)域有望得到更廣泛的應(yīng)用。