时频分析方法综述.pdf

-几种时频分析方法简介几种时频分析方法简介1 1 傅里叶变换（傅里叶变换（Fourier TransformFourier Transform）N1n j2nk/N2ftDFT:H()Th(kT)eFT:H(f)h(t)edtNT0离散化(离散取样)周期化（时频域截断）N12ft1kj2nk/NIFT:h(t)H(f)edtIDFT:h(nT)H()eNTN02 2 小波变换（小波变换（Wavelet TransformWavelet Transform）a.a.由傅里叶变换到窗口傅里叶变换（由傅里叶变换到窗口傅里叶变换（Gabor Transform(Short Time Fourier Transform)/Gabor Transform(Short Time Fourier Transform)/）从傅里叶变换的定义可知，时域函数h(t)的傅里叶变换 H(f)只能反映其在整个实轴的性态，不能反映 h（t）在特定时间区段内的频率变化情况。如果要考察h(t)在特定时域区间（比如：ta,b）内的频率成分，很直观的做法是将h(t)在区间 ta,b与函数1(t)1,ta,b，然后考察h(t)1(t)傅里叶变换。但是由于1(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)相乘得到的短时时域函数进行FT变换以考察 h(t)在特定时域内的频域情况。STFT：Gf(f,)h(t)g(t)e2ftdtISTFT：h(t)dfg(t)Gf(f,)e2ftd图：STFT 示意图STFTSTFT 算例算例 cos(210t)0s t 5scos(225t)5s t 10sx(t)=cos(250t)10s t 15scos(2100t)15s t 20s-图：四个余弦分量的 STFTb.b.窗口傅里叶变换（窗口傅里叶变换（GaborGabor）到小波变换（）到小波变换（Wavelet TransformWavelet Transform）图：小波变换-定义满足条件：ff2tdt0=0 tdt=0df 假定：的平方可积函数(t)(即(t)L2（，+）)为基本小波或小波母函数。基本小波或小波母函数。Haar 小波函数db3 小波函数db4 小波函数db5 小波函数-mexh 小波函数图：几种常用的小波函数令ab(t)1at b，a、b 为实数，且 a0，a称ab为由母函数母函数生成的有赖于参数 a,b 的连续小波函数。连续小波函数。设 f(t)L2（，+），定义其小波变换为：Wfa,b f,ab与 Fourier 类似，小波变化也具有反演公式：1at bftdtaft以及 Parseval 等式：1C Wfa,babtdadb，2a Wfa,bWga,b1C Wfa,b2dadbCf,g,2a2dadbftdt.a2小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点，但其在频率域上的分辨率小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点，但其在频率域上的分辨率却相应降低。却相应降低。这是小波变换的弱点，使它只能部分地克服Fourier 变换的局限性。小波包变换将在一定程度上弥补小波变换的这一缺陷。-图：FT 变换、STFT 变换及 Wavelet Analysis比较图：Wavelet应用 1探测数据突变点-图：Wavelet应用 1探测数据突变点（树状显示）图：Wavelet应用 2探测数据整体变化趋势-图：Wavelet应用 2探测数据中的频率成分图：Wavelet应用 3压缩数据-图：Wavelet应用 3压缩数据3 3 希尔伯特黄变换（希尔伯特黄变换（HilbertHilbert-Huang TransformHuang Transform）3.13.1 希尔伯特与瞬时频率（希尔伯特与瞬时频率（Hilbert Transform and instantaneous frequencyHilbert Transform and instantaneous frequency）对于任意一个时间序列X(t)，它的希尔伯特变换具有如下形式：Y(t)=1P-X()d,t-其中，P积分的柯西主值；希尔伯特变换对于任何属于Lp空间中的函数都存立，即上式中X(t)Lp（，+）。通过上述定义，X(t)和 Y(t)成为一组复共轭对，同时能够构造一个实部和虚部分为X(t)和 Y(t)的解析信号解析信号(Analytic Signal)(Analytic Signal)Z(t)，Z(t)表示为：Z(t)=X(t)iY(t)=ate其中，it,1/2Y(t)22at=X(t)+Y(t),t arctan.X(t)理论上讲有无数种方式去定义虚部，但是希尔伯特变换是唯一能够得到解析信号结果的方法。X(t)的Hilbert变换实质上是将X(t)与函数1/t在时域上做卷积，这就决定了通过X(t)的Hilbert变换能够考察其局部特性。得到X(t)的瞬时相位函数后，其瞬时频率为：d(t)wt.dt-32.521.51Amplitude0.50-0.5-1sinx0.5+sinx1.5+sinx1020304050time/s60708090100-1.5-20图：原始信号（三个正弦波）1.510.5Imaginary0-0.5-1-1.5-1-0.500.5Real11.522.5图：Hilbert 变换后解析信号的复平面图-4sinx0.5+sinx1.5+sinx32Istantaneousfrequency10-1-2-301020304050time/s60708090100图：三个正弦信号的瞬时频率3.23.2 经验模态分解与固有模态函数经验模态分解与固有模态函数（Empirical mode decomposition/Empirical mode decomposition/EMDEMD and Intrinsic modeand Intrinsic modefunction/function/IMFIMF）固有模态函数需要满足两个条件：（1）极值与零点的数量必须相等或最多相差一个；（2）由局部极大值包络和局部极小值包络定义的平均包络曲线上任何一点的值为0；A、EMD筛选过程（Sifting process）x(x(t t)m m1 1 h h1 1,h h1 1 m m2 2 h h2 2,.h hk k 1 1 m mk k h hk k.h hk k c c1 1.-图：原始数据图：极值包络与均值 m1-图：h1 与原始数据图：h1 与 m2-图：h3 与 m4图：h4 与 m5x(x(t t)c c1 1 r r1 1,r r1 1 c c2 2 r r2 2,.r rn n 1 1 c cn n r rn n.x(x(t t)c cj j 1 1n nj j r rn n.3.3 Hilbert3.3 Hilbert 谱与谱与 HilbertHilbert 边际谱边际谱经过筛选过程后，X(t)可以表示为 IMF 与残差量的和：-X(X(t t)C Cj j r rn n X X(t t)C C(t t)2 2 2 22 2j jj j 1 1T Tj j 1 1n nn n 1 1n n 1 1 n n 1 1j j 1 1 k k 1 1 C C(t t)C)Cj jk k(t t)IOIO t t 0 0 j j 1 1n n 1 1 n n 1 1 2 22 22 2C C(t t)C)C(t t)/X X(t t)0 0 X X(t t)C Cj j(t t)j jk kk k 1 1j j 1 1 n n 1 1对 X(t)的每一个 IMF 进行 Hilbert 变换可以得到 X(t)的 Hilbert 谱：HHTHHT:C Cj j(t t)a aj j(t t)e)ei i t t a aj j(t t)e)e i i j j t t dt dt X X(t t)C Cj j(t t)n nn na aj j(t t)e)e i i j j t t dt dt H Hj j(,t,t)n nj j 1 1j j 1 1HilbertHilbert Spectrum Spectrumn nF FT T:X X(t t)a a j jt tj j(t t)e)ei ij j 1 1得到 Hilbert 谱后可以进一步定义 Hilbert 边际谱：h(h()T T0 0H(H(,t,t)dt)dtHilbertHilbert Magrinal Spectrum Magrinal Spectrum算例算例 1 1：一个有跳变的余弦信号y cos(6t)t 10s5cos(6t)t 10s6号4信始2原0-202468101214161820时 间/s421C0-2-402468101214161820时 间/s105R0-502468101214161820时 间/s图 1:跳变信号及其分量-j j 1 1HilbHilbertert Spectru Spectrum m-400瞬时相位30020010000246810时 间/s1214161820瞬时频率数值方法13002001000-1000246810时 间/s1214161820X:6Y:18.59瞬时频率数值方法23002001000-1000246810时 间/s1214161820X:6Y:19.02图 2:跳变信号 EMD 分量的瞬时相位与频率算例算例 2 2：频率发生改变的余弦信号cos(6t)t 10sy cos(4t)t 10s1原始信号0-10246810时 间/s12141618202C10-20 x 10-3246810时 间/s121416182020R-20246810时 间/s1214161820图 3:频率改变余弦信号及其EMD 分解分量-400瞬时相位30020010000246810时 间/s1214161820瞬时频率数值方法120181614120246810时 间/s121412.5816182018.92瞬时频率数值方法220181614120246810时 间/s1212.581416182018.92图 4:频率改变余弦信号 IMF 分量瞬时相位与瞬时频率算例算例 3 3：余弦扫频信号y (10.2t)cos(4t2)0 t 10s420-2-4012345时 间/s678910原始信号C1420-2-4012345时 间/s6789100.050R-0.05012345时 间/s678910图 5：余弦扫频信号及其EMD 分解分量-1500瞬时相位10005000-500012345时 间/s678910瞬时频率数值方法14003002001000-100012345时 间/s678910瞬时频率数值方法24003002001000-100012345时 间/s678910图 6：余弦扫频信号 IMF 分量瞬时相位与瞬时频率算例算例 4 4：两个不同频率的正弦信号的叠加y sin(10t)sin(5t)0 t 10s2原始信号0-2012345时 间/s6789101C10-1012345时 间/s6789101C20-1012345时 间/s6789100.50-0.5R012345时 间/s678910图 7：两个不同频率叠加的正弦信号及其IMF 分量-100瞬时相位500-50012345时 间/s678910瞬时频率数值方法180604020001234X:4.26Y:10.155时 间/s678910瞬时频率数值方法2806040200012345时 间/s678910图 8：两个不同频率叠加的正弦信号IMF1 分量瞬时相位与瞬时频率6040200-20012345时 间/s678910瞬时频率数值方法1瞬时相位20151050012345时 间/s678910X:4.42Y:5.129瞬时频率数值方法220151050012345时 间/s678910图 9：两个不同频率叠加的正弦信号IMF2 分量瞬时相位与瞬时频率非线性问题求解非线性问题求解Duffing equationDuffing equation-d d2 2x xd d2 2x x3 3 x x x x 2 2 x x 1 1 x x2 2 cocos s t t.2 2d dt tdt dt 1 1 0.10.1 0.040.04 HzHzInitialInitial conditioncondition:x(x(o o),x x(0 0)11,1 1-熟悉 NCU Matlab HHT 程序：Function fa.mFunction fa.mInputfa(data,dt,ifmethod,normmethod,nfilter);data(n,k)其中 n 为数据长度，k 为 IMF 个数。Output freq,am;freq,am 均为 nk 矩阵The specifications of the calculating methods of the instantaneous frequencyThe specifications of the calculating methods of the instantaneous frequencyifmethodifmethodhilberthilbtmacoszcquadcosfornormmethodnormmethodnonesplineCalculating methodsCalculating methodsHilbert transformHilbert transformHilbert transformHilbert transformArcos methodArcos methodGeneralized zeroGeneralized zero-crossingcrossingmethodmethodQuadrature methodQuadrature methodCosine formula methodCosine formula methodNormalizationNormalizationmethodsmethodsNoneNoneSplineSplinenormalizationnormalizationFunction fileFunction fileFAhilbert.mFAhilbert.mFAimphilbert.mFAimphilbert.mFAcos.mFAcos.mFAzc.mFAzc.mFAquadrature.mFAquadrature.mFAcosfor.mFAcosfor.mThe normalization of input dataThe normalization of input dataRecommendedRecommendedRecommendedRecommendedNot requiredNot requiredNot requiredNot requiredRequiredRequiredNot recommendedNot recommendedRequiredRequiredRequiredRequiredThe normalized methods optionsThe normalized methods optionsFunction fileFunction fileNoneNonesplinenormalize.msplinenormalize.mRecommend how to useRecommend how to useFor zc optionFor zc optionNot for ensemble EMDNot for ensemble EMDmethodmethodFor hilbert or acos optionFor hilbert or acos optionNot for ensemble EMDNot for ensemble EMDmethodmethodWhen using Ensemble EMDWhen using Ensemble EMDmethodmethodReasonReasonPossiblePossibleovershotovershotPossiblePossibleovershotovershotdefaultdefaultsplineEPSplineSplinenormalizationnormalizationsplinenormalizeep.msplinenormalizeep.mwith several endwith several endprocessprocessHilbertHilbertamplitudeamplitudehilbertnormalizehilbertnormalize.m.mhilbert-normalizationnormalizationlinearLinearLinearnormaliztionnormaliztionCubic hermiteCubic hermitesplinesplinenormalizationnormalizationBlockBlocknormalizationnormalizationlinearnormalize.mlinearnormalize.mpchipnormalize.mpchipnormalize.mblocknormalize.mblocknormalize.mWhen using Ensemble EMDWhen using Ensemble EMDmethodmethodWhen using Ensemble EMDWhen using Ensemble EMDmethodmethodNot to useNot to usepchipblock-算例算例 1 1：（参见：参见：ex2012104.mex2012104.m）t2 t2tx(t)expcoscos32 0.3sin32 0 t 1024s2565123251264理论解推导过程如下：解析信号X(t)Ateit AtcostiAtsint x(t)ix(t)对比可知：AM（amplitude modulation）：At expt256 t2(t)32 0.3sin32Phase angle：6451232512 t2FM(frequency modulation):t2dtt0.3t(t)cos32dt16384819232 512f(t)t2DATA:Damped Chirp Duffing Model0.60.40.2Signalintensity0-0.2-0.4-0.6-0.80100200300400500Time/s6007008009001000图：原始信号-1TrueHTNHTDQ0.5Imaginary0-0.5-1-0.5Real00.51图：各种方法得到的解析信号与理论解析信号的复平面对比0.08TrueHTNHTDQ0.060.040.02IF0-0.02-0.04-0.06-0.08-0.10100200300400500Time/s6007008009001000图：三种不同方法得到的瞬时频率（IF）与理论瞬时频率对比-0.04TrueHTNHTDQ0.030.020.010-0.01520540560580600Time/s620640660680700trueHTNHTDQ图：三种不同方法得到的瞬时频率（FM）与理论瞬时频率对比（细节图）10.80.60.4Imaginatypart0.20-0.2-0.4-0.60100200300400500Time/s6007008009001000图：三种不同方法得到的解析信号虚部值与理论虚部分值对比-1TrueHTNHTDQ0.90.80.70.6AM0.50.40.30.20.100100200300400500Time/s6007008009001000图：三种不同方法得到的解析信号包络值（AM）与理论包络值对比结论：计算信号 IMF 分量的瞬时频率（FM）和包络（AM）采用 DQ 法效果最好。但是在使用DQ（direct quadrature method）法之前需要对信号的每一个固有模态分量（IMF component）进行两步关键的操作。第一步是将每一个固有模态分量进行标准化处理（得到nimf），然后在第二步再对其进行 AM-FM 分解。一个经过标准化的 IMF 分量进行 AM-FM 分解后满足：nimf AtFt Atcost假定：2Ft cost sin t 1 Ft上式正负号怎么确定呢？上式正负号怎么确定呢？-0.050.50.0450.450.040.40.0350.350.30.250.20.150.10.050.030.0250.020.0150.010.0050010020030040050060070080090010002.521.510.50-0.5-1-1.5-2-2.501020304050607080-0 80001600024000010002000300040005000600070008000 x 1012-41086420010002000300040005000600070008000-x 101.2-310.80.60.40.20010002000300040005000600070000.1120.09100.080.0780.060.050.0460.0340.0220.01001020304050607080-0.1120.09110.08100.0790.0680.0570.0460.0350.0240.013010203040506070802-