几种时频分析方法简介1． 傅里叶变换（Fourier Transform ）12/20122/0()()()()1()()()(::::)N j nk N ft N ft j nk N n H T h kT e H f h t e d DFT FT IFT IDFT t NT k h t H f e dt h nT H e N NT ππππ--∞--∞∞--∞⎫=⎫⎪=⋅⎪⎪−−−−−−−→⎬⎬⎪⎪=⋅=⎭⎪⎭∑⎰⎰∑离散化(离散取样)周期化（时频域截断） 2． 小波变换（Wavelet Transform ）a. 由傅里叶变换到窗口傅里叶变换（Gabor Transform(Short Time Fourier Transform)/）从傅里叶变换的定义可知，时域函数h(t)的傅里叶变换H(f )只能反映其在整个实轴的性态，不能反映h （t ）在特定时间区段内的频率变化情况。

如果要考察h(t)在特定时域区间（比如：t ∈[a,b]）内的频率成分，很直观的做法是将h(t)在区间t ∈[a,b]与函数[][]11,t ,()0,t ,a b t a b χ⎧∈⎪=⎨∈⎪⎩，然后考察1()()h t t χ傅里叶变换。

但是由于1()t χ在t= a,b 处突然截断，导致中1()()h t t χ出现了原来h （t ）中不存在的不连续，这样会使得1()()h t t χ的傅里叶变化中附件新的高频成分。

为克服这一缺点，D.Gabor 在1944年引入了“窗口”傅里叶变换的概念，他的做法是，取一个光滑的函数g(t),称为窗口函数，它在有限的区间外等于0或者很快地趋于0，然后将窗口函数与h(t)相乘得到的短时时域函数进行FT 变换以考察h(t)在特定时域内的频域情况。

22(,)()()()()(,)ft f ft f STFT ISTF G f h t g t e dth t df g t G f e d T ππτττττ+∞--∞+∞+∞-∞-∞=-=-⎰⎰⎰：：图：STFT 示意图STFT 算例cos(210) 0s t 5s cos(225) 5s t 10s (t)=cos(250) 10s t 15s cos(2100) 15s t 20st t x t t ππππ≤≤⎧⎪≤≤⎪⎨≤≤⎪⎪≤≤⎩图：四个余弦分量的STFTb. 窗口傅里叶变换（Gabor ）到小波变换（Wavelet Transform ）图：小波变换定义满足条件：()()()()2=ˆ=00ˆ0t dt t dt f df fψψψψ+∞-∞+∞<+∞-∞+∞-∞⎰<+∞−−−−−−→⇔⎰⎰假定：的平方可积函数ψ(t)(即ψ(t)∈L 2（—∞，+∞）)为——基本小波或小波母函数。

Haar 小波函数db3小波函数db4小波函数db5小波函数mexh 小波函数 图：几种常用的小波函数令()ab t b t a aψ-⎛⎫=⎪⎝⎭，a 、b 为实数，且a ≠0， 称ψab 为由母函数生成的有赖于参数a,b 的连续小波函数。

设f(t)∈L 2（—∞，+∞），定义其小波变换为：()(),,f ab t b W a b f f t dt a aψψ+∞-∞-⎛⎫==⎪⎝⎭⎰与Fourier 类似，小波变化也具有反演公式：()()()21,f ab dadbf t W a b t C aψψ+∞+∞-∞-∞=⎰⎰， 以及Parseval 等式：()()()()2222,,,,1,.f g f dadbW a b W a b C f g adadbW a b f t dt C a ψψ+∞+∞-∞-∞+∞+∞+∞-∞-∞-∞==⎰⎰⎰⎰⎰ 小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点，但其在频率域上的分辨率却相应降低。

这是小波变换的弱点，使它只能部分地克服Fourier 变换的局限性。

小波包变换将在一定程度上弥补小波变换的这一缺陷。

图：FT变换、STFT变换及Wavelet Analysis比较图：Wavelet应用1——探测数据突变点图：Wavelet应用1——探测数据突变点（树状显示）图：Wavelet应用2——探测数据整体变化趋势图：Wavelet应用2——探测数据中的频率成分图：Wavelet应用3——压缩数据图：Wavelet 应用3——压缩数据3． 希尔伯特—黄变换（Hilbert-Huang Transform ）3.1希尔伯特与瞬时频率（Hilbert Transform and instantaneous frequency ） 对于任意一个时间序列X(t)，它的希尔伯特变换具有如下形式：-1()(t)=,-X Y P d t ττπτ∞∞⎰其中，P ——积分的柯西主值；希尔伯特变换对于任何属于L p 空间中的函数都存立，即上式中X(t)∈L p （—∞，+∞）。

通过上述定义，X(t)和Y(t)成为一组复共轭对，同时能够构造一个实部和虚部分为X(t)和Y(t)的解析信号(Analytic Signal)Z(t)，Z(t)表示为：()()(t)=(t)(t)=a ,i t Z X iY t eθ+其中，()()1/222(t)a =(t)+(t),arctan .X(t)Y t X Y t θ⎛⎫⎡⎤= ⎪⎣⎦⎝⎭理论上讲有无数种方式去定义虚部，但是希尔伯特变换是唯一能够得到解析信号结果的方法。

X(t)的Hilbert 变换实质上是将X(t)与函数1/t 在时域上做卷积，这就决定了通过X(t)的Hilbert 变换能够考察其局部特性。

得到X(t)的瞬时相位函数后，其瞬时频率为：()()(t).d w t dtθ=3.2经验模态分解与固有模态函数（Empirical mode decomposition/EMD and Intrinsic mode function/IMF ）固有模态函数需要满足两个条件：（1）极值与零点的数量必须相等或最多相差一个；（2）由局部极大值包络和局部极小值包络定义的平均包络曲线上任何一点的值为0； A 、 EMD —筛选过程（Sifting process ）11122k 1k k k 1x(t )m h ,h m h ,..........h m h .h c .--=-=-==⇒图：原始数据图：极值包络与均值m1图：h1与原始数据图：h1与m2图：h3与m4图：h4与m511122n 1n n njn j 1x(t )c r ,r c r ,x(t )cr ...r c r ..-=-=-=-⇒=-=∑3.3 Hilbert 谱与Hilbert 边际谱经过筛选过程后，X(t)可以表示为IMF 与残差量的和：n n 1n 1n 122j n jjkj 1j 1j 1k 1Tn 1n 12j k t 0j 1k 1n 122j j 1X(t )C r X (t )C (t )2C (t )C(t )C (t )C (t )/X (t )IO X (t )C (t )0++++=====++====+⇒=+⎛⎫⇒=≈⇒ ⎪⎝⎭=∑∑∑∑∑∑∑∑对X(t)的每一个IMF 进行Hilbert 变换可以得到X(t)的Hilbert 谱：()()()j j j nni t dtj j j 1j 1Hilbert SpectrumHilb n i t dti t j j j j j 1ni ert Spectru tj j 1mC (t )a (t )ea (HHT :a (t )e X (t )C (t )t )e H (,t )X (t )a t T )eF :(ωωθωω====⎰==⇒====⎰∑∑∑∑得到Hilbert 谱后可以进一步定义Hilbert 边际谱：Hilbert Magrinal SpectrumTh()H(,t )dt ωω=⎰算例1：一个有跳变的余弦信号cos(6) 10t t s y π≤⎧=⎨算例2：频率发生改变的余弦信号cos(6)10t t sπ≤⎧算例3：余弦扫频信号算例4：两个不同频率的正弦信号的叠加非线性问题求解 Duffing equation()223222.10.10.04HzInitial condition :[x(o ),x d x d x x x x 1x co '(0)][1,1s t d dt ]t εγωεεγω++=++-=====熟悉NCU Matlab HHT程序：Function fa.mInput fa(data,dt,ifmethod,normmethod,nfilter); data(n,k)其中n为数据长度，k为IMF个数。

Output [freq,am]; freq ,am均为n×k矩阵算例1：（参见：ex2012104.m ）22()exp cos cos 320.3sin 32 0t 1024s 2566451232512t t t x t ππ⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=-+++≤≤ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭理论解推导过程如下： 解析信号()()()()()()()cos sin ()()i t X t A t e A t t iA t t x t ix t θθθ==+⎡⎤⎡⎤⎣⎦⎣⎦=+对比可知：AM （amplitude modulation ）：()exp 256t A t ⎛⎫=-⎪⎝⎭Phase angle ：22()320.3sin 32 6451232512t t t ππθ⎛⎫⎛⎫⎛⎫=+++ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭FM(frequency modulation):()20.3()cos 32 16384819232512d t t tt t dt θπππω⎛⎫⎛⎫==-++ ⎪ ⎪⎝⎭⎝⎭()()2t f t ωπ=图：原始信号计算信号IMF 分量的瞬时频率（FM ）和包络（AM ）采用DQ 法效果最好。

但是在使用DQ （direct quadrature method ）法之前需要对信号的每一个固有模态分量（IMF component ）进行两步关键的操作。

第一步是将每一个固有模态分量进行标准化处理（得到nimf ），然后在第二步再对其进行AM-FM 分解。

一个经过标准化的IMF 分量进行AM-FM 分解后满足：()()()()cos nimf A t F t A t t θ==⎡⎤⎣⎦假定：()()()cos sin F t t t θθ=⇒=⎡⎤⎡⎤⎣⎦⎣⎦上式正负号怎么确定呢？100200300400500600700800900100000.0050.010.0150.020.0250.030.0350.040.0450.050.050.10.150.20.250.30.350.40.450.501000200030004000500060007000800024681012x 10x 10-3...010203040506070800.010.020.030.040.050.060.070.080.090.12468101223456789101112。