時(shí)頻分析方法綜述

1、幾種時(shí)頻分析方法簡介1 .傅里葉變換(Fourier Transform)FT:H(f) h(t) e2 ftdtIFT:h(t)離散化(離散取樣) 周期化(時(shí)頻域截?cái)?H(f) e2 ftdtDFT:H需)IDFT:h( nT)N 1T h(kT)e j2 nk/N01 N 1kII I / k、j2 nk/NH (一)eN 0NT2 .小波變換(Wavelet Transforma.由傅里葉變換到窗口傅里葉變換( Gabor Tran sform(Short Time FourierTran sform)/ )從傅里葉變換的定義可知，時(shí)域函數(shù)h(t)的傅里葉變換H(f)只能反映其在整個(gè)實(shí)軸

2、的性態(tài)，不能反映h (t)在特定時(shí)間區(qū)段內(nèi)的頻率變化情況。如果要考察h(t)在特定時(shí)域區(qū)間(比如：t a,b)1(t)1'ta'b0,t/ a,b內(nèi)的頻率成分，很直觀的做法是將h(t)在區(qū)間t a,b與函數(shù)，然后考察h(t) 1(t)傅里葉變換。但是由于1(t)在t= a,b處突然截?cái)啵瑢?dǎo)致中h(t)1(t)出現(xiàn)了原來h(t)中不存在的不連續(xù)，這樣會使得h(t) i(t)的傅里葉變化中附件新的高頻成分。為克服這一缺點(diǎn),D.Gabor在1944年引入了 "窗g(t),稱為窗口函數(shù)，它在有限口傅里葉變換的概念，他的做法是，取一個(gè)光滑的函數(shù)h(t)相乘得到的短時(shí)時(shí)域函數(shù)進(jìn)的

3、區(qū)間外等于0或者很快地趨于 0,然后將窗口函數(shù)與行FT變換以考察h(t)在特定時(shí)域內(nèi)的頻域情況。STFT： Gf(f,) h(t)g(t )e2 ftdtISTFT： h(t)df g(t)Gf(f, )e2 ftdwindow pnH-dLU 反cos(2cos(2cos(2cos(210t)Ost5s25t)5st10s50t)10st15sloot)15st20s仙：IJb ru 1J U IE 1> V E hlICDjmib fl 幾 I 卅 N> k*忖|»a-u2&rT» wndftwl2&rTBwndcw吟 14 fW-T

4、4; HXUiMim 37s ruJOia 通 1： M u ifi 曲 i I可Ftttfmnr nnuJ圖：四個(gè)余弦分量的 STFT« Id 說 =l*(-ab.窗口傅里葉變換（Gabor ）到小波變換（Wavelet Transform*apnl-dLU 電-VI Wavelet I Transform ICUTimeWavelet Analysis圖：小波變換定義滿足條件:假定:t |dt? 0 =0t dt=O基本小波或小波母函數(shù)。的平方可積函數(shù)e (t)(即e (t) L2 (soe0.2書癥Il叼 hjndw ptii0.60.4HaarD00 20 4 O.e OSW

5、wiW ft>kctiDni|pai450 50 OF 04 O.e D»小波函數(shù)db3小波函數(shù)11r11O.Bo.b041 101<f.2”1 Ll_10hO.G10.2 J-11scaling hndlon ffiir/avtML itindion psiedb5小波函數(shù)mexh小波函數(shù)圖：幾種常用的小波函數(shù)ab(t)g , a、b為實(shí)數(shù)，且a丸, a稱e ab為由母函數(shù)生成的有賴于參數(shù)a,b的連續(xù)小波函數(shù)。設(shè) f(t) L2 (8,+ m),定義其小波變換為:Wf a,b f, ab)dta與Fourier類似，小波變化也具有反演公式:Wf a,babdadb2a以

6、及Parseval等式:dadbWf a,b 叫 a,b a2 dadb2aWf a,b小波變換雖然具有頻率愈高相應(yīng)時(shí)間或空間分辨率愈高的優(yōu)點(diǎn),f,g),2dt.但其在頻率域上的分辨率卻相應(yīng)降低。這是小波變換的弱點(diǎn)，使它只能部分地克服Fourier變換的局限性。小波包變換將在一定程度上彌補(bǔ)小波變換的這一缺陷。TimeTime Doman (Shannoini'i&LJmnbm2lAmplitudeFrequency Domaini i(Fourier)limeSIFT (Ga.bof)TimeWa-tfeletA n 且 I聲汩圖：FT變換、STFT變換及Wavelet An

7、alysis比較Signal «n4 ApprwjntfiQfHSIs*1 L心a05a03%1WOOhnjjjfiiiNJchyid.心1-0 2*0 5030,10-0 140060030Q1000圖：Wavelet應(yīng)用1 探測數(shù)據(jù)突變點(diǎn)Slgrul圖：Wavelet應(yīng)用1探測數(shù)據(jù)突變點(diǎn)（樹狀顯示）0 6324圖：Wavelet應(yīng)用2 探測數(shù)據(jù)整體變化趨勢aaard Appn)«linfll»nisF0恥0.030.02Coch, Signal and Dptall(s)1.-iiHniNiniiufflifmuuumumHiuuui.WwSt'fif

8、仲IftililU'l"llllr-3.2 心圖：Wavelet應(yīng)用cfs2U0*05 eod WD IODO2探測數(shù)據(jù)中的頻率成分-M14D1»專 yiA* 事抽 4hTiM*H m IM 2Wu -VfdrHl 會e «fff at 14 r.爭釦灣価ft =-.-_I 、-i hrwR 列PCMVfMrtWI Itwid 3圖：Wavelet應(yīng)用3壓縮數(shù)據(jù)« IDD IMno眄血如IT fl iML aklwtbKMpMHM M IMII 刖圖：Wavelet應(yīng)用3 壓縮數(shù)據(jù)3 .希爾伯特一黃變換(Hilbert-Hua ng Tran s

9、form3.1 希爾伯特與瞬時(shí)頻率(Hilbert Tran sform and in sta ntan eous freque ncy對于任意一個(gè)時(shí)間序列x(t)，它的希爾伯特變換具有如下形式:Y(t)=- P,t-其中，P 積分的柯西主值;希爾伯特變換對于任何屬于Lp空間中的函數(shù)都存立，即上式中X(t) Lp (8, + 8)。通過上述定義，X(t)和Y(t)成為一組復(fù)共軛對，同時(shí)能夠構(gòu)造一個(gè)實(shí)部和虛部分為X(t)和丫的解析信號(Analytic Signal)Z(t)，Z(t)表示為:Z(t)=X(t) iY(t)=a t其中,2 2 1/2a t = X2(t)+Y2(t)t arct

10、an 丫X(t)理論上講有無數(shù)種方式去定義虛部，但是希爾伯特變換是唯一能夠得到解析信號結(jié)果的方法。X(t)的Hilbert變換實(shí)質(zhì)上是將 X(t)與函數(shù)1/t在時(shí)域上做卷積，這就決定了通過X(t)的Hilbert變換能夠考察其局部特性。得到X(t)的瞬時(shí)相位函數(shù)后，其瞬時(shí)頻率為:wtdte m圖：原始信號(三個(gè)正弦波)Real圖：Hilbert變換后解析信號的復(fù)平面圖圖：三個(gè)正弦信號的瞬時(shí)頻率3.2經(jīng)驗(yàn)?zāi)B(tài)分解與固有模態(tài)函數(shù)（Empirical mode decomposition/EMD andIntrin sic mode fun cti on/IMF )固有模態(tài)函數(shù)需要滿足兩個(gè)條件：（1

11、）極值與零點(diǎn)的數(shù)量必須相等或最多相差一個(gè);(2)由局部極大值包絡(luò)和局部極小值包絡(luò)定義的平均包絡(luò)曲線上任何一點(diǎn)的值為A、EMD 篩選過程(Sifting process)x(t)h1 ,hk 1m2mkh2 ,hk .2MTirme . secondhkc1.圖：原始數(shù)據(jù)- =- J MOYMonEiT/«top«t and th* Mean : data2M3S»4&OMO560SOOData400Tiime . secondEnwiopeEfl曲ope Mn-an .m 1圖：極值包絡(luò)與均值m12M圖：h1與原始數(shù)據(jù)M7Arrplitude壽 oMOAr

12、rplitudaWOEnv«k>p4:nd 5a S4B3n oaf J- mss-CQa nn 帯-gn 54Bn .3fc-JT/KsT_ime -secondAsoEnv0lap4* Mnd 5a S4D3、%./ws300Kd*s2me -secondxssowseon函-h3 M 34o _&Env«k>p4ft and ttw Mean M450560400Tirme . secondSOO圖：h4與m5x(t)c1r1riC2r2，rn 1Cnrn .x(t)Cjj 1rn3.3 Hilbert 譜與 Hilbert 邊際譜經(jīng)過篩選過程后,

13、X(t)可以表示為IMF與殘差量的和:X(t)n 1CjrnX2(t)Cj2(t)j 11Cj(t)Ck(t)1IOn 1 n 1Cj(t)Ck(t)/X2(t)j 1 k 1X2(t)1Cj2(t)1對X(t)的每一個(gè)IMF進(jìn)行Hilbert變換可以得到X(t)的Hilbert 譜:HHT : Cj(t)aj(t)ei ti : t dtaj(t)eX(t)nCj (t)j 1i : t dte丿Hilbert SpectrumnHj( ,t)Hilb ertSpectrumi . tFT : X(t)aj(t)e jj 1得到Hilbert譜后可以進(jìn)一步定義Hilbert 邊際譜:Hilb

14、ertTh( )0 H(Magri nal Sp ectrum,t )dt算例1 :一個(gè)有跳變的余弦信號cos(6 t) t 10s y 5 cos(6 t) t 10s號信始原10時(shí)間/s121416182002468101214161820時(shí)間/s105R002468101214161820時(shí)間/s-5圖1:跳變信號及其分量400300位相時(shí)瞬200100002468101214161820時(shí)間/S1法方值數(shù)-率頻時(shí)瞬300rIL11 1-X:6Y:18.59-'亠 -III1LrrJ1 11020010002468101214161820時(shí)間/S-1002法方值數(shù)-率頻時(shí)瞬300

15、11p1 1LL11 11X6Y:19.021111rr11 11200100002468101214161820時(shí)間/S-100圖2:跳變信號EMD分量的瞬時(shí)相位與頻率算例2 :頻率發(fā)生改變的余弦信號cos(6 t) t 10s y cos(4 t) t 10sO號信始原CR圖3:頻率改變余弦信號及其 EMD分解分量時(shí)間/s118.9212.581法方值數(shù)-率頻時(shí)瞬位相時(shí)瞬2法方值數(shù)-率頻時(shí)瞬C CO O2 1圖4:頻率改變余弦信號IMF分量瞬時(shí)相位與瞬時(shí)頻率算例3 :余弦掃頻信號y (1 0.2t)cos(4 t2) 0 t 10s0123789456時(shí)間/S4號信始原10圖5 :余弦掃頻

16、信號及其 EMD分解分量位相時(shí)瞬1法方值數(shù)-率頻時(shí)瞬圖6 :余弦掃頻信號IMF分量瞬時(shí)相位與瞬時(shí)頻率算例4 :兩個(gè)不同頻率的正弦信號的疊加y si n(10t) sin (5t) 0 t 10s號信始原0123678945時(shí)間/S10-10訓(xùn)wwwwwv血12345678910C0.5-0.5045時(shí)間/S10圖7:兩個(gè)不同頻率疊加的正弦信號及其IMF分量OO O O O0 5 51 -位相時(shí)瞬1法方值數(shù)-率頻時(shí)瞬法方值數(shù)-率頻時(shí)瞬時(shí)間Is圖9 :兩個(gè)不同頻率疊加的正弦信號IMF2分量瞬時(shí)相位與瞬時(shí)頻率非線性問題求解Duffing equati ond2xNxdt21d2 x"dF&

17、quot;cos t .0.10.04 HzInitial condition :x(o),x'(0)1 ,1Dumrifl Equation : ODE?3tBIMF Duffngi Equahon : ODE23rB1 ft o & 16 0 -O OS) 20 X 1040SOSOy/100120HO1602040SOSO too 120T【hme tec endHO160laolao200200-2Duffing Equation : Hilbert SpectrumlU :ng EquMlon oO£±ftT A 4丄D4專5A»u

18、87;o« o4熟悉 NCU Matlab HHT 程序:data(n,k)其中n為數(shù)據(jù)長度，k為Function fa.mInputfa(data,dt,ifmethod, no rmmethod, nfilter);IMF個(gè)數(shù)。Out put freq,am;freq ,am 均為n xk矩陣The sp ecificati ons of the calculat ing methods of the in sta ntan eous freque ncyifmetho dCalculati ng methodsFun cti on fileThe no rmalizati on

19、 of input data hilbertHilbert tran sformFAhilbert.mRecomme ndedNot required hilbtmHilbert tran sformFAim philbert.mRecomme ndedNot requiredacos 'Arcos methodFAcos.mRequired zc'Gen eralizedzero-cross ing methodFAzc.mNot recomme ndedFAquadrature.'quad 'Quadrature methodmRequired cosfor

20、Cosi ne formulamethodFAcosfor.mRequiredThe no rmalized methods op ti ons normmethod 'Normalization methodsFun cti on fileRecomme nd how to useReason'none 'NoneNoneFor ' zc'op ti onP ossibleSp li nesp li neno rmalize.Not for en semble EMD spline 'overshono rmalizati onmmethodt

21、Sp li neFor ' hilbert ' or ' acos 'no rmalizati onsp li neno rmalizeeop tionP ossible'splineEP 'with severalp.mNot for en semble EMDovershoend p rocessmethodtHilberthilbert no rmalizeWhen using En semble'hilbert 'amp litude.mEMD methoddefaultno rmalizati onlinearno rm

22、alize.Lin earWhen using En semble'linear 'no rmaliztio nmEMD methodCubic hermiteWhen using En semble'pchip 'sp li nepchipno rmalize.mEMD methodno rmalizati onBlockblock no rmalize.m'block 'no rmalizati onNot to use算例 1 :(參見：ex2012104.m)32x(t) exp cos256cos64320.3si n 32工 320 t 1024s512理論解推導(dǎo)過程如下:解析信號X(t)A t cost iAt sin t對比可知:AM (amplitude modulationP hase an gle :(t)64FM(freque ncy modulati on):d t(t) *16384tx(t)：Aexp2560.3si