时频分析方法综述

1、* *几种时频分析方法简介1 .傅里叶变换(Fourier Transform)FT:H(f)IFT:h(t)h(t) e2 ftdtH(f) e2 ftdt离散化(离散取样) 周期化(时频域截断)DFT：H如IDFT:h(nT)N 1T h(kT)e0j2 nk/Nk、j2 nk/Nn 0 H(NT)e2 .小波变换(Wavelet Transform)a.由傅里叶变换到窗口傅里叶变换( Gabor Tran sform(Short Time FourierTran sform)/ )从傅里叶变换的定义可知，时域函数h(t)的傅里叶变换 H(f)只能反映其在整个实轴的性态，不能反映h (t)

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

3、T变换以考察h(t)在特定时域内的频域情况。STFT: Gf(f, ) h(t)g(t )e2 ftdtISTFT: h(t) df g(t)Gf(f, )e2 ftd* *windowmpn 兰 dETime图：STFT示意图STFT算例cos(210t)Ost5scos(225t)5st10sx(t)=cos(250t)10st15scos(2100t)15st20s1 i W0 M It ISA14parEiOT-MT T J% heSswtftOiaw! wi* T 13 vni26to windowMl12& wndwQ輩*Vi vh T 375 皿馆 I* 力吐冲 Wi nW T

4、1DOD rmft W 0 H 16 阳 r|v|BE K- dm:, tt- _ff o o Q DJQ图：四个余弦分量的 STFT* *b.窗口傅里叶变换(Gabor )到小波变换(Wavelet Transform)apnl_-dLUWI Wavelet BTransform图：小波变换定义满足条件：2? f假定：t|dt,df1? 0 =0t dt=0的平方可积函数” (t)(即“ (t) L2 (a,+ I )为一一基本小波或小波母函数。Selling FlindiiOn ptii0a.204 g 0Haar小波函数db3小波函数db4小波函数皿师lt funcliori 口并mex

5、h小波函数图：几种常用的小波函数ab(t)f , a、b为实数，且 a丸,.a a称“ ab为由母函数生成的有赖于参数a,b的连续小波函数。设 f(t) L2 (a,+ R),定义其小波变换为:Wf a,bdt与Fourier类似，小波变化也具有反演公式:Wf a,babdadb2-， a以及Parseval等式:dadbWf a,b W a,b 厂 ga22 dadb2aWf a,b小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点,相应降低。这是小波变换的弱点，使它只能部分地克服:fg,2dt.但其在频率域上的分辨率却Fourier变换的局限性。小波包变换Epru=duJp and t

6、he Man h3Deta . h3 Enwlope EnvelopeMan m4300碍0400450Tirpe %ond500SSOGOO图：h3与m4图：h4与m5x(t) Ciri ,r 1C2r2 ,rn 1Cnrn nx(t)Cj rn3.3 Hilbert 谱与 Hilbert 边际谱经过筛选过程后,X(t)可以表示为IMF与残差量的和:X(t)Cjn 12 2X2(t)Cj2(t)j 11Cj(t)Ck(t)1IOn 1 n 12Cj(t)Ck(t)/X2(t)j 1 k 1X2(t)1Cj2(t)1对X(t)的每一个IMF进行Hilbert变换可以得到X(t)的Hilbert

7、 谱:HHT : Cj(t)i ti j t dtaj(t)eaj(t)eX(t)nCjj 1(t)i : t dta(t)eHilbert SpectrumnHj( ,t) sJ Hilb ert Spectrumi tFT : X(t) a(t)e jj 1得到Hilbert谱后可以进一步定义Hilbert边际谱:Th( )0 H( ,t )dtHilbert Magrinal Spectrum算例1 :一个有跳变的余弦信号cos(6 t) t 10s5 cos(6 t) t 10s* *号信始原101214161820时间/s2C o-2-42468101214时间/s16 18 200

8、105R0-52468101214时间/s16 18 200图1:跳变信号及其分量位相时瞬004003002001002468101214时间/s1618 20法方值数一率频时瞬3002001000-1000X：；8.598101214时间/s1618 202法方值数一率频时瞬11IILLii-X6Y:19.02-11irni3002001000-1000246168101214时间/s18 20图2:跳变信号EMD分量的瞬时相位与频率算例2 :频率发生改变的余弦信号cos(6 t) t 10scos(4 t) t 10s* *位相时瞬20号信始原2 di416181CR时间/s图3:频率改变

9、余弦信号及其 EMD分解分量oo002468121416 18 203000200di10时间/s时间/s2法 方 值 数率 频 时 瞬20181614122468121416 18 20010时间/s图4:频率改变余弦信号IMF分量瞬时相位与瞬时频率算例3 :余弦扫频信号y (1 0.2t)cos(4 t2) 0 t 10sR图6 :余弦扫频信号IMF分量瞬时相位与瞬时频率图5 :余弦扫频信号及其 EMD分解分量* *算例4 :两个不同频率的正弦信号的叠加y si n(10t) si n(5t) 0 t 10soo号信始原o1 o1o06oR图7:两个不同频率叠加的正弦信号及其IMF分量00

10、位相时瞬50O50-法方值数一率频时瞬60-1.iLLLX: 4.26Y: 10.15rirrrriJlL604020123456789时间/s01040O时间/s图8 :两个不同频率叠加的正弦信号IMF1分量瞬时相位与瞬时频率时间/s时间/s时间/s图9 :两个不同频率叠加的正弦信号IMF2分量瞬时相位与瞬时频率非线性问题求解* *Dufing equati ond2xdt2x3d2 xdt2cos t .0.10.04 HzInitial condition :x(o),x(0)1 ,1PufTiingi ouatian : ODEZTBTinrte secondIMF Duffiirg

11、Equation : ODE23TB0M斗 060MJ1001201401601*02ODTlrr*e teeond* *Duffing Equation : Hilbert Spectrum* *熟悉 NCU Matlab HHT 程序:data(n,k)其中n为数据长度，k为Function fa.mIn putfa(data,dt,ifmethod, no rmmethod, nfilter);IMF个数。Output freq,am;freq ,am 均为n xk矩阵The specificati ons of the calculat ing methods of the in st

12、a ntan eous freque ncyifmetho dCalculati ng methodsFun cti on fileThe no rmalizati on of in put data hilbertHilbert tran sformFAhilbert.mRecomme ndedNot requiredhilbtmHilbert tran sformFAimphilbert.mRecomme ndedNot requiredacos Arcos methodFAcos.mRequired zcGen eralizedzero-cross ing methodFAzc.mNot

13、 recomme ndedFAquadrature.quad Quadrature methodmRequired cosforCosi ne formulamethodFAcosfor.mRequiredThe no rmalized methods opti ons normmethod Normalization methodsFun cti on fileRecomme nd how to useReasonnone NoneNoneFor zc opti onPossibleSpli nespli neno rmalize.Not for en semble EMD spline o

14、vershono rmalizati onmmethodtSpli neFor hilbert or acos no rmalizati onspli neno rmalizeeoptionPossiblesplineEP with severalp.mNot for en semble EMDovershoend processmethodtHilberthilbert no rmalizeWhen using En semblehilbert amplitude.mEMD methoddefaultno rmalizati onlinearno rmalize.Lin earWhen us

15、ing En semblelinear no rmaliztio nmEMD methodCubic hermiteWhen using En semblepchip spli nepchip normalize.mEMD methodno rmalizati onBlockblock no rmalize.mblock no rmalizati onNot to use* *算例 1 :(参见：ex2012104.m)x(t) exptt2cos cos -25664512320.3sin32 51232t 1024s2理论解推导过程如下: 解析信号X(t) A t e1 t A t cos

16、 t iA t sin t x(t) iX(t)对比可知：):A t expt256t2Phase an gle :(t)64 512 32AM ( amplitude modulation0.3si n FM(freque ncy modulati on):t232(t)d tdt32 512t0.3 tt2cos3216384819232512f(t)* *VFrs nnr anaps0.5-0.5DATA:Damped Chirp Duffi ng Model4 2 a ao2 a-4 6 a a - -0.801002003004005006007008009001000Time/s图

17、：原始信号-1-0.500.5Real图：各种方法得到的解析信号与理论解析信号的复平面对比0.080.060.040.02-0.02-0.04-0.06-0.08-0.10100 2003004005006007008009001000Time/s图：三种不同方法得到的瞬时频率（IF）与理论瞬时频率对比0.040.030.020.01-0.01520540560580600620640660680700Time/s图：三种不同方法得到的瞬时频率（FM ）与理论瞬时频率对比（细节图）ILrap ytanTQa 卩图：三种不同方法得到的解析信号虚部值与理论虚部分值对比图：三种不同方法得到的解析信号包络值（AM ）与理论包络值对比结论：计算信号IMF分量的瞬时频率（FM ）和包络（AM ）采用DQ法效果最好。但是在使用DQ （ direct quadrature method）法之前需要对信号的每一个固有模态分量（ IMFcomp onent ）进行两步关键的操作。第一步是将每一个固有模态分量进行标准化处理（得到nimf ），然后在第二步再对其进行 AM-FM 分解。一个经过标准化的IMF分量