时频分析方法综述

1、几种时频分析方法简介1 .傅里叶变换(Fourier Transform)FT:H(f) h(t) e2 ftdtIFT:h(t)离散化(离散取样) 周期化(时频域截断)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)/ )从傅里叶变换的定义可知，时域函数h(t)的傅里叶变换H(f)只能反映其在整个实轴

2、的性态，不能反映h (t)在特定时间区段内的频率变化情况。如果要考察h(t)在特定时域区间(比如：t a,b)1(t)1'ta'b0,t/ a,b内的频率成分，很直观的做法是将h(t)在区间t a,b与函数，然后考察h(t) 1(t)傅里叶变换。但是由于1(t)在t= a,b处突然截断，导致中h(t)1(t)出现了原来h(t)中不存在的不连续，这样会使得h(t) i(t)的傅里叶变化中附件新的高频成分。为克服这一缺点,D.Gabor在1944年引入了 "窗g(t),称为窗口函数，它在有限口傅里叶变换的概念，他的做法是，取一个光滑的函数h(t)相乘得到的短时时域函数进的

3、区间外等于0或者很快地趋于 0,然后将窗口函数与行FT变换以考察h(t)在特定时域内的频域情况。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图：四个余弦分量的 STFT« Id 说 =l*(-ab.窗口傅里叶变换（Gabor ）到小波变换（Wavelet Transform*apnl-dLU 电-VI Wavelet I Transform ICUTimeWavelet Analysis图：小波变换定义满足条件:假定:t |dt? 0 =0t dt=O基本小波或小波母函数。的平方可积函数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»小波函数db3小波函数11r11O.Bo.b041 101<f.2”1 Ll_10hO.G10.2 J-11scaling hndlon ffiir/avtML itindion psiedb5小波函数mexh小波函数图：几种常用的小波函数ab(t)g , a、b为实数，且a丸, a称e ab为由母函数生成的有赖于参数a,b的连续小波函数。设 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小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点,f,g),2dt.但其在频率域上的分辨率却相应降低。这是小波变换的弱点，使它只能部分地克服Fourier变换的局限性。小波包变换将在一定程度上弥补小波变换的这一缺陷。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应用1 探测数据突变点Slgrul图：Wavelet应用1探测数据突变点（树状显示）0 6324图：Wavelet应用2 探测数据整体变化趋势aaard Appn)«linfll»nisF0恥0.030.02Coch, Signal and Dptall(s)1.-iiHniNiniiufflifmuuumumHiuuui.WwSt'fif

8、仲IftililU'l"llllr-3.2 心图：Wavelet应用cfs2U0*05 eod WD IODO2探测数据中的频率成分-M14D1»专 yiA* 事抽 4hTiM*H m IM 2Wu -VfdrHl 會e «fff at 14 r.爭釦灣価ft =-.-_I 、-i hrwR 列PCMVfMrtWI Itwid 3图：Wavelet应用3压缩数据« IDD IMno眄血如IT fl iML aklwtbKMpMHM M IMII 刖图：Wavelet应用3 压缩数据3 .希尔伯特一黄变换(Hilbert-Hua ng Tran s

9、form3.1 希尔伯特与瞬时频率(Hilbert Tran sform and in sta ntan eous freque ncy对于任意一个时间序列x(t)，它的希尔伯特变换具有如下形式:Y(t)=- P,t-其中，P 积分的柯西主值;希尔伯特变换对于任何属于Lp空间中的函数都存立，即上式中X(t) Lp (8, + 8)。通过上述定义，X(t)和Y(t)成为一组复共轭对，同时能够构造一个实部和虚部分为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)理论上讲有无数种方式去定义虚部，但是希尔伯特变换是唯一能够得到解析信号结果的方法。X(t)的Hilbert变换实质上是将 X(t)与函数1/t在时域上做卷积，这就决定了通过X(t)的Hilbert变换能够考察其局部特性。得到X(t)的瞬时相位函数后，其瞬时频率为:wtdte m图：原始信号(三个正弦波)Real图：Hilbert变换后解析信号的复平面图图：三个正弦信号的瞬时频率3.2经验模态分解与固有模态函数（Empirical mode decomposition/EMD andIntrin sic mode fun cti on/IMF )固有模态函数需要满足两个条件：（1

11、）极值与零点的数量必须相等或最多相差一个;(2)由局部极大值包络和局部极小值包络定义的平均包络曲线上任何一点的值为A、EMD 筛选过程(Sifting process)x(t)h1 ,hk 1m2mkh2 ,hk .2MTirme . secondhkc1.图：原始数据- =- J MOYMonEiT/«top«t and th* Mean : data2M3S»4&OMO560SOOData400Tiime . secondEnwiopeEfl曲ope Mn-an .m 1图：极值包络与均值m12M图：h1与原始数据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 边际谱经过筛选过程后,

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)的每一个IMF进行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谱后可以进一步定义Hilbert 边际谱:Hilb

14、ertTh( )0 H(Magri nal Sp ectrum,t )dt算例1 :一个有跳变的余弦信号cos(6 t) t 10s y 5 cos(6 t) t 10s号信始原10时间/s121416182002468101214161820时间/s105R002468101214161820时间/s-5图1:跳变信号及其分量400300位相时瞬200100002468101214161820时间/S1法方值数-率频时瞬300rIL11 1-X:6Y:18.59-'亠 -III1LrrJ1 11020010002468101214161820时间/S-1002法方值数-率频时瞬300

15、11p1 1LL11 11X6Y:19.021111rr11 11200100002468101214161820时间/S-100图2:跳变信号EMD分量的瞬时相位与频率算例2 :频率发生改变的余弦信号cos(6 t) t 10s y cos(4 t) t 10sO号信始原CR图3:频率改变余弦信号及其 EMD分解分量时间/s118.9212.581法方值数-率频时瞬位相时瞬2法方值数-率频时瞬C CO O2 1图4:频率改变余弦信号IMF分量瞬时相位与瞬时频率算例3 :余弦扫频信号y (1 0.2t)cos(4 t2) 0 t 10s0123789456时间/S4号信始原10图5 :余弦扫频

16、信号及其 EMD分解分量位相时瞬1法方值数-率频时瞬图6 :余弦扫频信号IMF分量瞬时相位与瞬时频率算例4 :两个不同频率的正弦信号的叠加y si n(10t) sin (5t) 0 t 10s号信始原0123678945时间/S10-10训wwwwwv血12345678910C0.5-0.5045时间/S10图7:两个不同频率叠加的正弦信号及其IMF分量OO O O O0 5 51 -位相时瞬1法方值数-率频时瞬法方值数-率频时瞬时间Is图9 :两个不同频率叠加的正弦信号IMF2分量瞬时相位与瞬时频率非线性问题求解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为数据长度，k为Function fa.mInputfa(data,dt,ifmethod, no rmmethod, nfilter);IMF个数。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理论解推导过程如下:解析信号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