时频分析方法综述

几种时频分析方法简介

1． 傅里叶变换（Fourier Transform）

N1

nj2nk/N

2ftDFT:H()Th(kT)eFT:H(f)h(t)edtNT0离散化(离散取样)

周期化（时频域截断）N1

2ft1kj2nk/N IFT:h(t)H(f)edtIDFT:h(nT)H()eNTN0

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]与函数

1(t)

1,ta,b，然后考察h(t)1(t)傅里叶变换。但是由于1(t)在t= a,b处突然0,ta,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)e2ftdt





ISTFT：h(t)df







g(t)Gf(f,)e2ftd

图：STFT示意图

STFT算例

cos(210t) 0st5scos(225t) 5st10s

x(t)=

cos(250t) 10st15scos(2100t) 15st20s

图：四个余弦分量的STFT

b. 窗口傅里叶变换（Gabor）到小波变换（Wavelet Transform）

图：小波变换

定义满足条件：





ˆff

2



tdt

ˆ0=0tdt=0

假定：





的平方可积函数ψ(t)(即ψ(t)∈L2（—∞，+∞）)为——基本小波或小波母函数。

Haar小波函数

db3小波函数

db4

小波函数

db5小波函数

mexh小波函数 图：几种常用的小波函数

令

ab(t)



tb

，a、b为实数，且a≠0， a

称ψab为由母函数生成的有赖于参数a,b的连续小波函数。设f(t)∈L2（—∞，+∞），定义其小波变换为：

Wf

a,bf,ab

与Fourier类似，小波变化也具有反演公式：







tbftdt

a

ft

以及Parseval等式：

1

C









Wfa,babt

dadb

， 2

a









Wfa,bWga,b

1C





Wfa,b

2

dadb

Cf,g,2

a

2dadb

ftdt.a2

小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点，但其在频率域上的分辨率

却相应降低。这是小波变换的弱点，使它只能部分地克服Fourier变换的局限性。小波包变换将在一定程度上弥补小波变换的这一缺陷。

图：FT变换、STFT变换及Wavelet Analysis比较

图：Wavelet应用1——探测数据突变点

图：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)，它的希尔伯特变换具有如下形式：

Y(t)=

1



P



-

X()

, t-

其中，P——积分的柯西主值；

希尔伯特变换对于任何属于Lp空间中的函数都存立，即上式中X(t)∈Lp（—∞，+∞）。 通过上述定义，X(t)和Y(t)成为一组复共轭对，同时能够构造一个实部和虚部分为X(t)和Y(t)的解析信号(Analytic Signal)Z(t)，Z(t)表示为：

Z(t)=X(t)iY(t)=ateit,

其中，

1/2Y(t)22

at=X(t)+Y(t),tarctan. X(t)

理论上讲有无数种方式去定义虚部，但是希尔伯特变换是唯一能够得到解析信号结果的方法。

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

wt

d(t)

.

dt

3.2经验模态分解与固有模态函数（Empiricalmode decomposition/EMDand Intrinsic mode function/IMF）

固有模态函数需要满足两个条件：（1）极值与零点的数量必须相等或最多相差一个；（2）由局部极大值包络和局部极小值包络定义的平均包络曲线上任何一点的值为0； A、 EMD—筛选过程（Sifting process）

x(t)m1h1,h1m2h2,..........

hk1mkhk.hkc1.

图：原始数据

图：极值包络与均值m1

图：h1与原始数据

图：h1与

m2

图：h3与m4

图：h4与m5

x(t)c1r1,r1c2r2,...rn1cnrn.x(t)

n

c

j1

j

rn.

3.3Hilbert谱与Hilbert边际谱

经过筛选过程后，X(t)可以表示为IMF与残差量的和：

X(t)CjrnX(t)C(t)2

2

2j

j1T

j1

nn1n1n1

j1k1

C(t)C

j

k

(t)

n1

IO

t0j1



k1

n1

n122

Cj(t)Ck(t)/X(t)0X(t)C2j(t)

j1

对X(t)的每一个IMF进行Hilbert变换可以得到X(t)的Hilbert谱：

HHT:Cj(t)aj(t)e

n

it

aj(t)e

ijtdt

X(t)Cj(t)

j1

nn

aj(t)e

ijtdt

j1

Hilbert Spectrum

Hj(,t)

j1

Hilbert Spectrum

n

FT:X(t)aj(t)e

j1

ijt

得到Hilbert谱后可以进一步定义Hilbert边际谱：

Hilbert Magrinal Spectrum

h()H(,t)dt

T

算例1：一个有跳变的余弦信号

cos(6t) t10

s

y

算例2：频率发生改变的余弦信号

cos(6t)t10

s

算例3：余弦扫频信号

算例4：两个不同频率的正弦信号的叠加

非线性问题求解 Duffing equation

d2xd2x3

xx2x1x2cost.2

dtdt1

0.1

0.04Hz

Initialcondition:

[x(o),x'(0)][1,1]

熟悉NCU Matlab HHT程序： Function fa.m

Input fa(data,dt,ifmethod,normmethod,nfilter); data(n,k)其中n为数据长度，k为IMF个数。 Output [freq,am]; freq,am均为n×k矩阵

算例1：（参见：ex2012104.m）

t2t2tx(t)expcoscos320.3sin32 0t1024s 2565123251264

理论解推导过程如下： 解析信号

X(t)Ate

对比可知：

it

AtcostiAtsint

(t) x(t)ix

t

 256

AM（amplitude modulation）：Atexp

t2

(t)320.3sin32Phase angle：

6451232512

t2

FM(frequency modulation):

t2dtt0.3t

(t)cos32

dt[1**********]512f(t)

t

2

图：原始信号

计算信号IMF分量的瞬时频率（FM）和包络（AM）采用DQ法效果最好。但是在使用DQ（direct quadrature method）法之前需要对信号的每一个固有模态分量（IMF component）进行两步关键的操作。第一步是将每一个固有模态分量进行标准化处理（得到nimf），然后在第二步再对其进行AM-FM分解。一个经过标准化的IMF分量进行AM-FM分解后满足：

nimfAtFtAtcost

假定：

Ftcostsint

上式正负号怎么确定呢？

几种时频分析方法简介

1． 傅里叶变换（Fourier Transform）

N1

nj2nk/N

2ftDFT:H()Th(kT)eFT:H(f)h(t)edtNT0离散化(离散取样)

周期化（时频域截断）N1

2ft1kj2nk/N IFT:h(t)H(f)edtIDFT:h(nT)H()eNTN0

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]与函数

1(t)

1,ta,b，然后考察h(t)1(t)傅里叶变换。但是由于1(t)在t= a,b处突然0,ta,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)e2ftdt





ISTFT：h(t)df







g(t)Gf(f,)e2ftd

图：STFT示意图

STFT算例

cos(210t) 0st5scos(225t) 5st10s

x(t)=

cos(250t) 10st15scos(2100t) 15st20s

图：四个余弦分量的STFT

b. 窗口傅里叶变换（Gabor）到小波变换（Wavelet Transform）

图：小波变换

定义满足条件：





ˆff

2



tdt

ˆ0=0tdt=0

假定：





的平方可积函数ψ(t)(即ψ(t)∈L2（—∞，+∞）)为——基本小波或小波母函数。

Haar小波函数

db3小波函数

db4

小波函数

db5小波函数

mexh小波函数 图：几种常用的小波函数

令

ab(t)



tb

，a、b为实数，且a≠0， a

称ψab为由母函数生成的有赖于参数a,b的连续小波函数。设f(t)∈L2（—∞，+∞），定义其小波变换为：

Wf

a,bf,ab

与Fourier类似，小波变化也具有反演公式：







tbftdt

a

ft

以及Parseval等式：

1

C









Wfa,babt

dadb

， 2

a









Wfa,bWga,b

1C





Wfa,b

2

dadb

Cf,g,2

a

2dadb

ftdt.a2

小波变换虽然具有频率愈高相应时间或空间分辨率愈高的优点，但其在频率域上的分辨率

却相应降低。这是小波变换的弱点，使它只能部分地克服Fourier变换的局限性。小波包变换将在一定程度上弥补小波变换的这一缺陷。

图：FT变换、STFT变换及Wavelet Analysis比较

图：Wavelet应用1——探测数据突变点

图：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)，它的希尔伯特变换具有如下形式：

Y(t)=

1



P



-

X()

, t-

其中，P——积分的柯西主值；

希尔伯特变换对于任何属于Lp空间中的函数都存立，即上式中X(t)∈Lp（—∞，+∞）。 通过上述定义，X(t)和Y(t)成为一组复共轭对，同时能够构造一个实部和虚部分为X(t)和Y(t)的解析信号(Analytic Signal)Z(t)，Z(t)表示为：

Z(t)=X(t)iY(t)=ateit,

其中，

1/2Y(t)22

at=X(t)+Y(t),tarctan. X(t)

理论上讲有无数种方式去定义虚部，但是希尔伯特变换是唯一能够得到解析信号结果的方法。

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

wt

d(t)

.

dt

3.2经验模态分解与固有模态函数（Empiricalmode decomposition/EMDand Intrinsic mode function/IMF）

固有模态函数需要满足两个条件：（1）极值与零点的数量必须相等或最多相差一个；（2）由局部极大值包络和局部极小值包络定义的平均包络曲线上任何一点的值为0； A、 EMD—筛选过程（Sifting process）

x(t)m1h1,h1m2h2,..........

hk1mkhk.hkc1.

图：原始数据

图：极值包络与均值m1

图：h1与原始数据

图：h1与

m2

图：h3与m4

图：h4与m5

x(t)c1r1,r1c2r2,...rn1cnrn.x(t)

n

c

j1

j

rn.

3.3Hilbert谱与Hilbert边际谱

经过筛选过程后，X(t)可以表示为IMF与残差量的和：

X(t)CjrnX(t)C(t)2

2

2j

j1T

j1

nn1n1n1

j1k1

C(t)C

j

k

(t)

n1

IO

t0j1



k1

n1

n122

Cj(t)Ck(t)/X(t)0X(t)C2j(t)

j1

对X(t)的每一个IMF进行Hilbert变换可以得到X(t)的Hilbert谱：

HHT:Cj(t)aj(t)e

n

it

aj(t)e

ijtdt

X(t)Cj(t)

j1

nn

aj(t)e

ijtdt

j1

Hilbert Spectrum

Hj(,t)

j1

Hilbert Spectrum

n

FT:X(t)aj(t)e

j1

ijt

得到Hilbert谱后可以进一步定义Hilbert边际谱：

Hilbert Magrinal Spectrum

h()H(,t)dt

T

算例1：一个有跳变的余弦信号

cos(6t) t10

s

y

算例2：频率发生改变的余弦信号

cos(6t)t10

s

算例3：余弦扫频信号

算例4：两个不同频率的正弦信号的叠加

非线性问题求解 Duffing equation

d2xd2x3

xx2x1x2cost.2

dtdt1

0.1

0.04Hz

Initialcondition:

[x(o),x'(0)][1,1]

熟悉NCU Matlab HHT程序： Function fa.m

Input fa(data,dt,ifmethod,normmethod,nfilter); data(n,k)其中n为数据长度，k为IMF个数。 Output [freq,am]; freq,am均为n×k矩阵

算例1：（参见：ex2012104.m）

t2t2tx(t)expcoscos320.3sin32 0t1024s 2565123251264

理论解推导过程如下： 解析信号

X(t)Ate

对比可知：

it

AtcostiAtsint

(t) x(t)ix

t

 256

AM（amplitude modulation）：Atexp

t2

(t)320.3sin32Phase angle：

6451232512

t2

FM(frequency modulation):

t2dtt0.3t

(t)cos32

dt[1**********]512f(t)

t

2

图：原始信号

计算信号IMF分量的瞬时频率（FM）和包络（AM）采用DQ法效果最好。但是在使用DQ（direct quadrature method）法之前需要对信号的每一个固有模态分量（IMF component）进行两步关键的操作。第一步是将每一个固有模态分量进行标准化处理（得到nimf），然后在第二步再对其进行AM-FM分解。一个经过标准化的IMF分量进行AM-FM分解后满足：

nimfAtFtAtcost

假定：

Ftcostsint

上式正负号怎么确定呢？