一種多信號(hào)頻帶的瞬時(shí)頻率估計(jì)方法

一種多信號(hào)頻帶的瞬時(shí)頻率估計(jì)方法

非穩(wěn)定信號(hào)的時(shí)差頻率是很多實(shí)際處理中非常重要的工作。Boashash等系統(tǒng)地闡述了瞬時(shí)頻率的概念和一些瞬時(shí)頻率估計(jì)的方法,如中心有限差分法,相位建模法,時(shí)頻分布法等。在這些方法中,非時(shí)頻方法對(duì)噪聲比較敏感,所以應(yīng)用場合受到很大限制。時(shí)頻方法可以在時(shí)頻二維平面上給出信號(hào)的時(shí)變頻譜,能夠在時(shí)頻平面用遮隔去掉只在頻域?yàn)V波無法濾除的干擾和噪聲,使瞬時(shí)頻率估計(jì)精度顯著提高。然而,時(shí)頻方法運(yùn)算過程復(fù)雜,而且運(yùn)算量大,僅Wigner-Ville分布的計(jì)算復(fù)雜度就達(dá)到O(N2logN)。另外,對(duì)于信號(hào)和噪聲存在時(shí)頻耦合,即信號(hào)和噪聲的投影在時(shí)域和頻域都有重疊的情況下,時(shí)頻方法無法將噪聲遮隔掉。近年來,由于分?jǐn)?shù)階Fourier變換(fractionalFouriertransform,FRFT)的獨(dú)特優(yōu)點(diǎn),FRFT已經(jīng)迅速發(fā)展為信號(hào)處理的重要工具。例如,對(duì)上述時(shí)頻耦合的情況,只要信號(hào)和噪聲線性可分,就可以在合適的FRFT域取得好的濾波和干擾分離效果。同時(shí)也出現(xiàn)了多種利用信號(hào)的FRFT研究瞬時(shí)頻率估計(jì)的方法。W.X.Cong等利用信號(hào)模值的平方和Fresnel譜,提出的瞬時(shí)頻率估計(jì)和相位提取方法、,由于采用的是遞歸算法,使得運(yùn)算量和時(shí)間開銷都很大,激發(fā)了人們尋找非遞歸算法的興趣。文獻(xiàn)等利用FRFT的Wigner-Ville分布在FRFT域上的旋轉(zhuǎn)進(jìn)行相位提取,需要在[0,π]的大角度范圍內(nèi)進(jìn)行FRFT運(yùn)算,計(jì)算量非常大并且有些情況下無法實(shí)現(xiàn)。本文中利用任意確定階數(shù)的分?jǐn)?shù)階Fourier譜和信號(hào)相位微分的關(guān)系,提出了一種估計(jì)信號(hào)瞬時(shí)頻率的新方法,該方法不需要迭代運(yùn)算并且運(yùn)算量低,估計(jì)精度高。1frft反變換的核信號(hào)x(t)在α=pπ2上的FRFT定義為Xα(u)=Fα[x(t)](u)=∫∞-∞Κ(α,t,u)x(t)dt(1)其中Xα(u)表示信號(hào)x(t)的FRFT變換,Fα表示FRFT算子,α=pπ2,0＜|p|＜2。變換核K(α,t,u)為Κ(α,t,u)={√1-jcotα2πexp(jt2+u22cotα-jutcscα),若a≠nπδ(t-u),若α=2nπδ(t+u),若α=(2n+1)π(2)當(dāng)α=0和α=π時(shí),K(α,t,u)退化成δ(t-u)和δ(t+u)。從(2)式還可看出,FRFT反變換的核K(-α,t,u)=K*(α,t,u).由于分?jǐn)?shù)階Fourier變換對(duì)應(yīng)著信號(hào)的Wigner-Ville分布Wx(t,f)=∫∞-∞x(t+12τ)x*(t-12τ)exp(-j2πfτ)dτ(3)和模糊函數(shù)A(τ,v)=∫∞-∞x(t+12τ)x*(t-12τ)exp(-j2πvt)dτ(4)在時(shí)頻平面的旋轉(zhuǎn),所以FRFT可以看成將信號(hào)在時(shí)間軸上逆時(shí)針旋轉(zhuǎn)角度α到u軸上的表示,因而FRFT也被認(rèn)為是一種時(shí)頻分析方法。當(dāng)α=0時(shí),FRFT就是時(shí)域信號(hào)自身,當(dāng)α=π/2和α=3π/2時(shí),FRFT對(duì)應(yīng)著經(jīng)典的Fourier變換Xπ/2(u)和Fourier反變換X-π/2(u)。FRFT算子一個(gè)重要性質(zhì)是具有相加性,即Fα+β=FαFβ,另一個(gè)重要性質(zhì)是,分?jǐn)?shù)階功率譜|Xα(u)|2,即信號(hào)的FRFT模值的平方,等于信號(hào)的Wigner-Ville分布Wx(t,f)在u軸上的投影|Xα(u)|2=∫∞-∞Wx(ucosα-fsinα,usinα+fcosα)df(5)關(guān)于FRFT的其它性質(zhì)可參見文獻(xiàn)。2變量范圍估計(jì)方法FRFT是經(jīng)典Fourier變換的廣義形式,在信號(hào)處理中有著廣泛的應(yīng)用,一般認(rèn)為,能使用Fourier變換的場合就能使用FRFT。設(shè)信號(hào)x(t)的Fourier變換為Xπ/2(u),則Fourier反變換為x(t)=12π∫∞-∞Xπ/2(u)ejutdu(6)在上式中兩邊對(duì)t進(jìn)行微分有dx(t)dt=12π∫∞-∞juXπ/2(u)ejutdu(7)整理得F-π/2[Xπ/2(u)u](t)=-idx(t)dt(8)考慮到FRFT的可相加性以及x(t)=X0(t),不難得出F-π/2[Xα+π/2(u)u](t)=-idXα(t)dt(9)這就把β=α+π/2階的FRFT域的濾波表示為α階FRFT的微分。由FRFT的性質(zhì),Fα[x(u)u](t)=(tcosα+isinαD)Fα[x(u)](t)(10)可得F-α[Xα(u)u](t)=(tcosα-isinαD)x(t)(11)其中D=d/dt為微分算子。顯然,從(11)和(8)有F?α[X±α(u)u](t)=cosαtx(t)±isinαdx(t)dt=cosαF0[X0(u)u](t)±sinαF-π/2[Xπ/2(u)u](t)(12)把階數(shù)+α和-α分別代入(12)式,聯(lián)立后解方程組可解得信號(hào)的微分dx(t)dt=iF-π/2[Xπ/2(u)u](t)和tx(t)=F0[X0(u)u](t)如下-idx(t)dt=12sinα{F-α[Xα(u)u](t)-Fα[X-α(u)u](t)}tx(t)=12cosα{F-α[Xα(u)u](t)+Fα[X-α(u)u](t)}(13)設(shè)x(t)=a(t)ejφ(t),φ(t)是信號(hào)x(t)的相位,則從(12)式可以求得|F?α[X±α(u)u](t)|2=sin2α(da(t)dt)2+|a(t)|2(tcosα±sinαdφ(t)dt)2注意到|a(t)|2=|x(t)|2,則兩FRFT功率譜之差為|F-α[Xα(u)u](t)|2-|Fα[X-α(u)u](t)|2=2(sin2α)t|x(t)|2dφ(t)dt(14)即dφ(t)dt=12(sin2α)t|x(t)|2{|F-α[Xα(u)u](t)|2-|Fα[X-α(u)u](t)|2}(15)其中α是定值。本文中提出的方法就是利用(15)式來進(jìn)行瞬時(shí)頻率估計(jì),由兩個(gè)分?jǐn)?shù)階功率譜和信號(hào)模值的平方|x(t)|2就可以重構(gòu)瞬時(shí)頻率fΙF(t)=12πdφ(t)dt。一般來說信號(hào)x(t)可以由它的幅度和瞬時(shí)頻率重構(gòu)出來,則本文方法說明由兩個(gè)分?jǐn)?shù)階功率譜和信號(hào)模值就可以重構(gòu)出相位不為常量的信號(hào)。由于FRFT是一種一維的線性算子,其快速算法可用FFT實(shí)現(xiàn),計(jì)算復(fù)雜度為O(NlogN),遠(yuǎn)低于Wigner-Ville分布的O(N2logN),因此該方法運(yùn)算量低。作為一種不需要迭代運(yùn)算的瞬時(shí)頻率估計(jì)方法,該方法需要的數(shù)據(jù)量少,適用范圍廣,實(shí)現(xiàn)過程簡單。另外,在某些場合下,FRFT域的濾波及信噪分離等更為有效,可以更方便地應(yīng)用該方法。3frft信號(hào)轉(zhuǎn)換與光譜圖2假設(shè)信號(hào)是均勻采樣的,采樣間隔為T,并假設(shè)數(shù)據(jù)的FRFT功率譜|F-α[Xα(u)u](nT)|2和|Fα[X-α(u)u](nT)|2是已知的。這兩個(gè)分?jǐn)?shù)階功率譜可以通過計(jì)算信號(hào)FRFT模值的平方得到,也可通過Wigner-Ville變換與FRFT的關(guān)系計(jì)算得到,在光學(xué)中還可以直接測(cè)量得到。將信號(hào)瞬時(shí)頻率的估計(jì)值記為?fΙF(t),則用于估計(jì)信號(hào)瞬時(shí)頻率的(15)式的離散形式為?f(nΤ)={|F-α[Xα(u)u](nΤ)|2-|Fα[X-α(u)u](nΤ)|2}2(sin2α)(nΤ)|x(nΤ)|2(16)4仿真實(shí)驗(yàn)和結(jié)果下面通過兩個(gè)例子,分別針對(duì)不含噪聲和含噪聲的兩種仿真信號(hào)進(jìn)行瞬時(shí)頻率估計(jì),驗(yàn)證了本方法的有效性。另外,為了演示算法的估計(jì)精度,還將估計(jì)出的瞬時(shí)頻率用于信號(hào)的重構(gòu),并與原信號(hào)進(jìn)行了對(duì)比。應(yīng)用瞬時(shí)頻率重構(gòu)信號(hào)的方法為x^(nΤ)=|x(nΤ)|exp[j∑m=-Μn2πf^ΙF(mΤ)Τ]?Μ的取值滿足n<-M時(shí)fIF(nT)=0。對(duì)于噪聲中的信號(hào),由于無法直接進(jìn)行波形對(duì)比,因而給出了原信號(hào)的和重構(gòu)后信號(hào)的pseudoWigner-Ville分布(PWVD)的圖形。對(duì)FRFT的運(yùn)算采用的是文獻(xiàn)中的算法。(1)考慮信號(hào)x(t)=exp(-5t2)exp[i∫-∞tsin(5πt)dt]?t∈[-0.5,0.5],采樣點(diǎn)數(shù)為256,采樣率為100Hz。圖1(a)給出了原始瞬時(shí)頻率fIF(nT)的和由(16)式重構(gòu)的f^ΙF(nΤ)的圖形,圖1(b)給出了原始和根據(jù)瞬時(shí)頻率重構(gòu)的相位圖形。圖1(c)給出了原始信號(hào)和重構(gòu)信號(hào)的圖形。從圖中可見,本文算法不僅有效而且具有很高的精度,重構(gòu)出的信號(hào)與原始的非常接近。(2)考慮噪聲中的信號(hào):x(t)=exp[-(2.25t)8]×{Aexp[i∫-∞t[800arctan(25t)+400π]dt]+v(t)},其中v(t)為高斯復(fù)值白噪聲,信噪比為SNR=20log(A/σv),t∈[-0.5,0.5]。圖2給出了SNR=10dB時(shí)的仿真結(jié)果。圖2(a)給出了原始瞬時(shí)頻率fIF(nT)的和由(16)式重構(gòu)的f^(nΤ)的圖形,圖2(b)給出原始相位和重構(gòu)相位的圖形。圖2(c)給出了原始信號(hào)的PWVD圖形,圖2(d)給出了重構(gòu)信號(hào)的PWVD圖形。從圖2的PWVD圖形可以看出,對(duì)于噪聲中的信號(hào),本文方法同樣能夠有效地提取瞬時(shí)頻率,由于噪聲的原因,估計(jì)精度有所下降。隨著信噪比的降低,瞬時(shí)頻率估計(jì)的質(zhì)量逐步降低。通過多次實(shí)驗(yàn)發(fā)現(xiàn),當(dāng)信噪比高于3dB時(shí),估計(jì)結(jié)果能夠滿足要求;當(dāng)信噪比低于3dB時(shí),本方法不再有效。另外,隨著α取值的不同,估計(jì)的結(jié)果會(huì)有少許不同。相