..幾種時頻分析方法簡介傅里葉變換〔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