时频分析方式综述

几种时频分析方式简介
1． 傅里叶变换（Fourier Transform ）
1
2/201
22/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 χ的傅里叶转变中附件新的高频成份。

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

22(,)()()()()(,)ft f ft
f STFT ISTF G f h t
g t e dt
h 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 20s
t t x t t ππππ≤≤⎧⎪≤≤⎪⎨≤≤⎪⎪≤≤⎩
图：四个余弦分量的STFT
b. 窗口傅里叶变换（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 、b 为实数，且a ≠0， 称ψab 为由母函数生成的有赖于参数a,b 的持续小波函数。

设f(t)∈L 2
（—∞，+
∞）
，概念其小波变换为：
(
)(),,f ab t b W a b f f t dt a ψψ+∞
-∞
-⎛⎫
==
⎪⎝⎭
与Fourier 类似，小波转变也具有反演公式：
()()()
2
1
,f ab dadb
f t W a b t C a ψ
ψ+∞+∞
-∞
-∞
=
⎰⎰
， 和Parseval 等式：
(
)()
()
()2
22
2,,,,1,.f g f dadb
W a b W a b C f g a dadb
W 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）
希尔伯特与瞬时频率（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/2
22
(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
θ=
图：Hilbert 变换后解析信号的复平面图
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与m5
11122n 1n n n
j
n j 1
x(t )c r ,r c r ,x(t )c
r ...
r c r ..
-=-=-=-
⇒=-=∑
Hilbert 谱与Hilbert 边际谱
通过挑选进程后，X(t)能够表示为IMF 与残差量的和：
n n 1n 1n 12
2j n j
j
k
j 1j 1
j 1k 1
T
n 1
n 1
2j k t 0j 1k 1n 12
2j 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 n n
i t dt
j j j 1
j 1
Hilbert Spectrum
Hilb n i t dt i t j j j j j 1
n
i ert Spectru t
j j 1
m
C (t )a (t )e
a (HHT :a (t )e X (t )C (t )t )e H (,t )
X (t )a t T )e
F :(ωωθωω====⎰==⇒====⎰∑∑∑∑
取得Hilbert 谱后能够进一步概念Hilbert 边际谱：
Hilbert Magrinal Spectrum
T
h()
H(,t )dt ωω=⎰
算例1：一个有跳变的余弦信号
cos(6) 10
t t s y π≤⎧
=⎨
算例2：频率发生改变的余弦信号
cos(6)10
t t s
π≤⎧
算例3：余弦扫频信号
算例4：两个不同频率的正弦信号的叠加
非线性问题求解
Duffing equation
()
22
32
22
.
1
0.1
0.04Hz
Initial condition:
[x(o),x
d x d x
x x x1x co
'(0)][1,1
s t d dt
]
t
ε
γ
ω
εεγω++=++
-
=
=
=
=
=
Function
Input fa(data,dt,ifmethod,normmethod,nfilter); data(n,k)其中n为数据长度，k为IMF个数。

Output [freq,am]; freq ,am均为n×k矩阵
The specifications of the calculating methods of the instantaneous frequency
The normalized methods options
算例1：（参见：）
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 ：2
2()320.3sin 32 6451232512
t t t ππθ⎛⎫
⎛⎫⎛⎫=+++ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭
FM(frequency modulation):
()2
0.3()cos 32 16384819232512
d t t t
t 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.005
0.01
0.015
0.02
0.0250.03
0.035
0.04
0.045
0.050.05
0.1
0.150.20.250.30.35
0.4
0.45
0.5
010002000300040005000600070008000
2
4
6
8
10
12
x 10x 10-3
01020304050607080
00.01
0.02
0.030.040.05
0.06
0.07
0.08
0.09
0.12
46
8
1012 2
3
456789101112。