时频分析方法综述

..几种时频分析方法简介傅里叶变换〔FourierTransform小波变换〔WaveletTransform由傅里叶变换到窗口傅里叶变换〔GaborTransform<ShortTimeFourierTransform>/从傅里叶变换的定义可知,时域函数h<t>的傅里叶变换H<f>只能反映其在整个实轴的性态,不能反映h〔t在特定时间区段内的频率变化情况。如果要考察h<t>在特定时域区间〔比如：t∈[a,b]内的频率成分,很直观的做法是将h<t>在区间t∈[a,b]与函数,然后考察傅里叶变换。但是由于在t=a,b处突然截断,导致中出现了原来h〔t中不存在的不连续,这样会使得的傅里叶变化中附件新的高频成分。为克服这一缺点,D.Gabor在1944年引入了"窗口"傅里叶变换的概念,他的做法是,取一个光滑的函数g<t>,称为窗口函数,它在有限的区间外等于0或者很快地趋于0,然后将窗口函数与h<t>相乘得到的短时时域函数进行FT变换以考察h<t>在特定时域内的频域情况。图：STFT示意图STFT算例图：四个余弦分量的STFT窗口傅里叶变换〔Gabor到小波变换〔WaveletTransform图：小波变换定义满足条件：的平方可积函数ψ<t><即ψ<t>∈L2〔—∞,+∞>为——基本小波或小波母函数。Haar小波函数db3小波函数db4小波函数db5小波函数mexh小波函数图：几种常用的小波函数令,a、b为实数,且a≠0,称ψab为由母函数生成的有赖于参数a,b的连续小波函数。设f<t>∈L2〔—∞,+∞,定义其小波变换为：与Fourier类似,小波变化也具有反演公式：,以及Parseval等式：小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点,但其在频率域上的分辨率却相应降低。这是小波变换的弱点,使它只能部分地克服Fourier变换的局限性。小波包变换将在一定程度上弥补小波变换的这一缺陷。图：FT变换、STFT变换及WaveletAnalysis比较图：Wavelet应用1——探测数据突变点图：Wavelet应用1——探测数据突变点〔树状显示图：Wavelet应用2——探测数据整体变化趋势图：Wavelet应用2——探测数据中的频率成分图：Wavelet应用3——压缩数据图：Wavelet应用3——压缩数据希尔伯特—黄变换〔Hilbert-HuangTransform3.1希尔伯特与瞬时频率〔HilbertTransformandinstantaneousfrequency对于任意一个时间序列X<t>,它的希尔伯特变换具有如下形式：其中,P——积分的柯西主值；希尔伯特变换对于任何属于Lp空间中的函数都存立,即上式中X<t>∈Lp〔—∞,+∞。通过上述定义,X<t>和Y<t>成为一组复共轭对,同时能够构造一个实部和虚部分为X<t>和Y<t>的解析信号<AnalyticSignal>Z<t>,Z<t>表示为：其中,理论上讲有无数种方式去定义虚部,但是希尔伯特变换是唯一能够得到解析信号结果的方法。X<t>的Hilbert变换实质上是将X<t>与函数1/t在时域上做卷积,这就决定了通过X<t>的Hilbert变换能够考察其局部特性。得到X<t>的瞬时相位函数后,其瞬时频率为：图：原始信号〔三个正弦波图：Hilbert变换后解析信号的复平面图图：三个正弦信号的瞬时频率3.2经验模态分解与固有模态函数〔Empiricalmodedecomposition/EMDandIntrinsicmodefunction/IMF固有模态函数需要满足两个条件：〔1极值与零点的数量必须相等或最多相差一个；〔2由局部极大值包络和局部极小值包络定义的平均包络曲线上任何一点的值为0；EMD—筛选过程〔Siftingprocess图：原始数据图：极值包络与均值m1图：h1与原始数据图：h1与m2图：h3与m4图：h4与m53.3Hilbert谱与Hilbert边际谱经过筛选过程后,X<t>可以表示为IMF与残差量的和：对X<t>的每一个IMF进行Hilbert变换可以得到X<t>的Hilbert谱：得到Hilbert谱后可以进一步定义Hilbert边际谱：算例1：一个有跳变的余弦信号图1:跳变信号及其分量图2:跳变信号EMD分量的瞬时相位与频率算例2：频率发生改变的余弦信号图3:频率改变余弦信号及其EMD分解分量图4:频率改变余弦信号IMF分量瞬时相位与瞬时频率算例3：余弦扫频信号图5：余弦扫频信号及其EMD分解分量图6：余弦扫频信号IMF分量瞬时相位与瞬时频率算例4：两个不同频率的正弦信号的叠加图7：两个不同频率叠加的正弦信号及其IMF分量图8：两个不同频率叠加的正弦信号IMF1分量瞬时相位与瞬时频率图9：两个不同频率叠加的正弦信号IMF2分量瞬时相位与瞬时频率非线性问题求解Duffingequation熟悉NCUMatlabHHT程序：Functionfa.mInputfa<data,dt,ifmethod,normmethod,nfilter>;data<n,k>其中n为数据长度,k为IMF个数。Output[freq,am];freq,am均为n×k矩阵ThespecificationsofthecalculatingmethodsoftheinstantaneousfrequencyifmethodCalculatingmethodsFunctionfileThenormalizationofinputdata‘hilbert’HilberttransformFAhilbert.mRecommendedNotrequired‘hilbtm’HilberttransformFAimphilbert.mRecommendedNotrequired‘acos’ArcosmethodFAcos.mRequired‘zc’Generalizedzero-crossingmethodFAzc.mNotrecommended‘quad’QuadraturemethodFAquadrature.mRequired‘cosfor’CosineformulamethodFAcosfor.mRequiredThenormalizedmethodsoptions‘normmethod’NormalizationmethodsFunctionfileRecommendhowtouseReason‘none’NoneNoneFor‘zc’option‘spline’Splinenormalizationsplinenormalize.mNotforensembleEMDmethodPossibleovershot‘splineEP’Splinenormalizationwithseveralendprocesssplinenormalizeep.mFor‘hilbert’or‘acos’optionNotforensembleEMDmethodPossibleovershot‘hilbert’Hilbertamplitudenormalizationhilbertnormalize.mWhenusingEnsembleEMDmethoddefault‘linear’Linearnormaliztionlinearnormalize.mWhenusingEnsembleEMDmethod‘pchip’Cubichermitesplinenormalizationpchipnormalize.mWhenusingEnsembleEMDmethod‘block’Blocknormalizationblocknormalize.mNottouse算例1：〔参见：ex2012104.m理论解推导过程如下：解析信号对比可知：AM〔amplitudemodulation：Phaseangle：FM<frequencymodulation>:图：原始信号图：各种方法得到的解析信号与理论解析信号的复平面对比图：三种不同方法得到的瞬时频率〔IF与理论瞬时频率对比图：三种不同方法得到的瞬时频率〔FM与理论瞬时频率对比〔细节图图：三种不同方法得到的解析信号虚部值与理论虚部分值对比图：三种不同方法得到的解析信号包络值〔AM与理论包络值对比结论：计算信号IM