时频分析

在信号处理中，时频分析（time–frequency analysis）是指同时在时域和频域对信号进行研究的技术，其使用各种时频表示（representations）。

时频分布是一项让我们能够同时观察一个信号的时域和频域资讯的工具，而时频分析就是在分析时频分布。传统上，我们常用傅里叶变换来观察一个信号的频谱。然而，这样的方法不适合用来分析一个频率会随着时间而改变的信号，由于傅里叶变换只分析了一维的信号分布，而时频分析却能分析二维(时域跟频域)的信号分布，因此在信号处理中更常被运用。

时频分析也可以说是傅里叶分析的一般化，通常用于频率特性会随时间而变化的信号上，而在日常生活中符合符合此特性的信号非常多，像是演讲、音乐、影像、医学信号等，因此能应用的领域相当广泛。

另外，更实际应用时频分析的动机为传统傅里叶分析假设信号在时域是无限长或是周期性出现的，然而在现实中许多信号都只有短暂的存在，而且在信号持续期间可能有相当大的变化。举例来说，传统的音乐乐器不会持续产生无限长的正弦波，反而可能突然有一巨声，然后渐渐减弱。因此时频分析的研究势不可挡。

让我们看看以下这个频率会随时间变化的信号例子：

$x(t)={\begin{cases}\cos(\pi t);&t<10\\\cos(3\pi t);&10\leq t<20\\\cos(2\pi t);&t>20\end{cases}}$

一旦这样的数学式成立，便可利用时频分析的各种技术，萃取信号中的各种有用资讯，并分离噪音或干扰。

历史

最早的时频分析方法应见于Alfréd Haar提出的哈尔小波转换(1909)，然而在当时因时频分析所需的运算量仍是个无法忽视的议题，因此并未广泛应用于信号处理。而后更多的贡献来自于加博尔·德奈什，像是小波前身Gabor原子(1947)，以及加伯转换和改进型的短时距傅里叶变换。维格纳准几率分布(Ville 1948)也是一个重要的开端。

特别在1930年代及1940年代，早期的时频分析方法恰好与量子力学的发展一致，这反映了位置-动量平面及时域-频域平面的数学机制有些共通性，像是海森堡不确定性原理(量子力学)与加伯限制(时频分析)最终都得出了扭对称几何结构。

常见的时频分布函数

常见的时频分布函数有短时距傅里叶变换（包含加伯转换）、科恩分布函数（包含韦格纳分布）、改进型韦格纳分布 ，以及加伯-韦格纳分布（Gabor-Wigner distribution function）函数及S转换等。

而这些看似不同的时频分析函数，其数学公式的由来都有些相关性，若想对时频分析的了解更加透彻，应在学习时将它们一起理解，而非都视为单一函数，像是做1/4次傅里叶变换可以解读成傅里叶变换在时频分析平面上转90°，而这个旋转做了4次后就会回到原本的函数，只做2次时则会视反转的图形。

理想的时频分布函数

一个理想的时频分布函数有助于我们做时频分析，而它大致上具有以下四种性质：

1. “高清晰度”：可让我们分析更容易。
2. “没有交叉项（cross-term）”：可避免我们把信号和噪声混淆。
3. “好的数学性质”：有利于我们在许多方面的应用。
4. “较低的运算复杂度”：使得我们分析的速度变快。

在这里我们比较几个较常用的时频分析之优劣度。

	清晰度	交叉项	好的数学性质	运算复杂度
加伯转换	较差	无	较差	低
韦格纳分布函数	最好	有	最好	高
加伯-韦格纳分布函数	好	几乎可以消除	好	高
锥状分布	好	无	好	中

为了能顺利的分析各信号之时频分布，选择适当的时频分布函数是很重要的。而至于要如何选择时频分布函数呢？这端看于我们所要应用它的地方在哪边。韦格纳分布之定义中的自相关函数是一把双面刃，它让韦格纳分布函数拥有高的清晰度，然而，它也同时让它产生了交叉项的问题。

因此，如果我们想要分析一个只有单一项的信号，此时不会有交叉项的产生，因此我们通常选择韦格纳分布函数来获得高清晰度；另一方面，如果我们要分析的信号是由很多个项所组成的，此时若用韦格纳分布会有交叉项产生，所以我们可能选择用加伯转换或是加伯-韦格纳分布函数会比较好。

应用

在接下来即将介绍的应用中，我们除了需要时频分布函数，还需要搭配其他的运算才能达到目的，而著名的线性标准转换（linear canonical transform，LCT）可以帮助我们。我们可以利用线性标准转换来任意的改变一个信号在时频分布平面上面的形状和位置，像是水平以及垂直的移动、扩大、扭曲（shearing），以及旋转（用分数傅里叶变换，fractional Fourier transform, FRFT）等。由此可见，线性标准转换让我们对于时频分布的处理更灵活。

这边我们列举一些时频分布之应用的例子。

找出瞬间频率

瞬间频率的定义是 ${\frac {1}{2\pi }}{\frac {d}{dt}}\phi (t)$，其中$\phi (t)$ 是信号的瞬时相位。我们可以直接由时频分布的图形中看出每个时刻的瞬时频率是多少，不过前提是这个时频分布的图形要够清晰，因此，我们经常选用韦格纳分布函数来做进一步的分析。

滤波器设计

滤波器的目的就是要移除我们不要的部分，并保留我们要的部分。在没有应用时频分布之前，我们只能分别在时域跟频域上面来做过滤的动作，如下所示。

像上面这样只能分别在时域或频域上过滤的方式，并不适合处理每一种信号。如果信号在时域上或在频域上有重叠的话，这时候使用时频分布函数来做分析过滤，并搭配线性完整转换的操作，就可以做出更有效且灵活的滤波器。让我们看看以下的例子。

而在滤波器设计的应用中，时频分布通常处理的信号是由很多个项所组成的，因此若用韦格纳分布来做分析的话，将会产生交叉项的问题。或许加伯转换、加伯-韦格纳分布函数，亦或科恩类分布函数会是比较好的选择。

信号分解

信号分解的概念就跟滤波器设计很类似。

取样定理

由 奈奎斯特－香农采样定理且经过一番推导，我们大致上可以说一个信号经过取样后而不产生失真（aliasing）的最低取样点数，会和这信号在时频平面上图形的面积相等（事实上，没有一个信号在时频平面上的面积有限的，因此我们省略了一些精确度）。接下来，让我们看看传统取样定理跟结合了时频分析以后的取样定理之差异。

若浅绿色的部分是我们取样的涵盖范围，则我们可以很明显的看出使用时频分析后，所需取样的点数会比之前少了许多，因此加快了我们的运算。当我们使用韦格纳分布函数时，可能会产生交叉项；另一方面，若使用加伯转换做分析的话，又可能会因为清晰度不佳而让所需要取样的面积又变大了。因此，选用哪个函数要视信号的情形而定，如果信号是单一项组成的，那么就使用韦格纳分布函数；然而，如果信号是由多项组成的，则用加伯转换、加伯-韦格纳分布函数，或是科恩类分布函数。

取样方法

取样点数(sampling points) = 时频分析面积的总和 + 其余额外参数

• 如何使时频分布的面积更小?

1. 将原本信号切割成数个部分。
2. 使用啁啾乘法(chirp multiplications)、啁啾卷积(chirp convolutions)、分数傅里叶转换或线性标准转换减小面积。

Step 1. 解析信号转换

转换信号到座标轴的同一边(一般是取该信号的实数区)

$x(t)\rightarrow {\begin{alignedat}{2}x_{a}(t)&=x(t)+jx_{H}(t)\\x_{a}(f)&=x(f)+jH(f)x(f)\\\end{alignedat}},\ x_{H}(t):{\text{Hilbert transform of }}x(t)$

$H(f)={\begin{cases}-j,&{\text{for }}f>0\\j,&{\text{for }}f<0\end{cases}},x_{a}(f)={\begin{cases}2x(f),&{\text{for }}f>0\\0,&{\text{for }}f<0\end{cases}}$

x(f) v.s. x_a(f) analytic signal conversion first step

Step 2. 信号拆解

使用短时距傅里叶变换(因为信号包含许多不同成分)来拆解信号成许多部分。

separate the frequency components

Step 3. 使用斜推(shearing)或旋转(rotation)使各个部分减少到最小"面积"

使用韦格纳分布方程(Wigner distribution function, WDF; 因为此时信号为单一成分且属随机程序)来斜推和翻转各个部分。

Step 4. 使用传统采样理论采样各个成分

• 传统的取样方式

$x_{d}[n]=x(n\vartriangle _{t}),\vartriangle _{t}<1/F$

• 重建

$x(t)=\sum _{n}x_{d}[n]sinc({\frac {t}{\vartriangle _{t}}}-n)$

• 新的取样方式

1. $x(t)\rightarrow x_{a}(t)=x(t)+jx_{H}(t),\ x_{H}(t):{\text{Hilbert transform of }}x(t)$
2. $x(t)\rightarrow x_{a}(t)=x_{1}(t)+x_{2}(t)+...+x_{K}(t)$
3. $y_{k}(t)=exp(j2\pi a_{k}t^{2})x_{k}(t),k=1,2,...,K$
4. ${\begin{aligned}x_{d,k}[n]&=y_{k}(n\vartriangle _{t,k})\\&=exp(j2\pi a_{k}n^{2}\vartriangle _{t,k}^{2})x_{k}(n\vartriangle _{t,k}),k=1,2,...,K\\\end{aligned}}$

• 重建

1. $y_{k}(t)=\sum _{n}x_{d,k}[n]sinc({\frac {t}{\vartriangle _{t,k}}}-n)$
2. $x_{k}(t)=exp(-j2\pi a_{k}t^{2})y_{k}(t)$
3. $x_{a}(t)=x_{1}(t)+x_{2}(t)+...+x_{K}(t)$
4. $x(t)=\Re [x_{a}(t)]$

严格来说，没有任何一个信号的时频分布"面积"是有限的，但是我们可以选择一个阈值Δ，使时频分析 $\left\vert X(t,f)\right\vert >\bigtriangleup$或是分布的面积是有限的。但若以"面积"来讨论取样点数，也会牺牲一些精确度。

• 理论

如果$x(t)$是时间上有限的$(x(t)=0{\text{ for }}t<t_{1}{\text{ and }}t>t_{2})$, 则频率上则不可能是有限的。

如果$x(t)$是频率上有限的$(X(f)=0{\text{ for }}f<f_{1}{\text{ and }}f>f_{2})$,则时间上则不可能是有限的。

只取$t\in [t_{1},t_{2}]$和$f\in [f_{1},f_{2}]$牺牲的能量所占的比例，$err={\frac {\int _{-\infty }^{t_{1}}|x(t)|^{2}dt+\int _{t_{2}}^{\infty }|x(t)|^{2}dt+\int _{-\infty }^{f_{1}}|X_{1}(f)|^{2}df+\int _{f_{2}}^{\infty }|X_{1}(f)|^{2}df}{\int _{-\infty }^{\infty }|x(t)|^{2}dt}}$

$X_{1}(f)=FT[x_{1}(t)],x_{1}(t)=x(t){\text{ for }}t\in [t_{1},t_{2}],x_{1}(t)=0{\text{ otherwise }}$

• 韦格纳分布方程 (Wigner distribution function, WDF)

$|x(t)|^{2}=\int _{-\infty }^{\infty }W_{x}(t,f)df,|X(f)|^{2}=\int _{-\infty }^{\infty }W_{x}(t,f)dt$

$\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }W_{x}(t,f)dfdt=\int _{-\infty }^{\infty }|x(t)|^{2}dt={\text{ energy of }}x(t)$

${\begin{aligned}&\int _{-\infty }^{t_{1}}|x(t)|^{2}dt+\int _{t_{2}}^{\infty }|x(t)|^{2}dt+\int _{-\infty }^{f_{1}}|X_{1}(f)|^{2}df+\int _{f_{2}}^{\infty }|X_{1}(f)|^{2}df\\&=\int _{-\infty }^{t_{1}}\int _{-\infty }^{\infty }W_{x}(t,f)dfdt+\int _{t_{2}}^{\infty }\int _{-\infty }^{\infty }W_{x}(t,f)dfdt+\int _{-\infty }^{\infty }\int _{-\infty }^{f_{1}}W_{x_{1}}(t,f)dfdt+\int _{-\infty }^{\infty }\int _{f_{2}}^{\infty }W_{x_{1}}(t,f)dfdt\\&=\int _{-\infty }^{t_{1}}\int _{-\infty }^{\infty }W_{x}(t,f)dfdt+\int _{t_{2}}^{\infty }\int _{-\infty }^{\infty }W_{x}(t,f)dfdt+\int _{t_{1}}^{t_{2}}\int _{-\infty }^{f_{1}}W_{x_{1}}(t,f)dfdt+\int _{t_{1}}^{t_{2}}\int _{f_{2}}^{\infty }W_{x_{1}}(t,f)dfdt\\&\cong \int _{-\infty }^{t_{1}}\int _{-\infty }^{\infty }W_{x}(t,f)dfdt{\text{ (A) }}+\int _{t_{2}}^{\infty }\int _{-\infty }^{\infty }W_{x}(t,f)dfdt{\text{ (B) }}+\int _{t_{1}}^{t_{2}}\int _{-\infty }^{f_{1}}W_{x}(t,f)dfdt{\text{ (C) }}+\int _{t_{1}}^{t_{2}}\int _{f_{2}}^{\infty }W_{x}(t,f)dfdt{\text{ (D) }}\end{aligned}}$

$err\cong 1-{\frac {\int _{t_{1}}^{t_{2}}\int _{f_{1}}^{f_{2}}W_{x}(t,f)dfdt}{\int _{-\infty }{\infty }|x(t)|^{2}dt}}$

summation of area of Wigner distribution function on time-frequency analysis

调变与多工

传统上，调变（modulation）与 多工（multiplexing）都只有分别在时域及频域上下功夫，也就是尽量塞满时域及频域上的空间，这都是一维的操作。如果我们利用时频分布函数，就可以将调变与多工的触角延伸至二维的时频平面上，所要做的就是塞满整个时频平面，做最有效的利用。由以下例子可以让我们更了解。

由上例可知，使用韦格纳分布来分析会有严重的交叉项问题，这非常不利于调变与多工的作业，因此不能选择它来做这种应用。

通常可借由以下方法来实行调变和多工，

(1)加伯转换(Gabor transform)或加伯–韦格纳转换(the Gabor-Wigner transform)

(2)水平和垂直的位移(horizontal and vertical shifting)、扩张(dilation)、斜推(shearing)、广义斜推(generalized shearing)和旋转(rotation)

• 传统调变方法

信号$x_{1}(t),x_{2}(t),x_{3}(t),......,x_{K}(t)$可以被成功传递，如果${\begin{aligned}&{\text{Allowed Bandwidth}}\geq \sum _{k=1}^{K}B_{k}\\&B_{k}{\text{:the bandwidth (including the negative frequency part) of }}x_{k}(t)\\\end{aligned}}$

• 基于时频分析的调变方法

信号$x_{1}(t),x_{2}(t),x_{3}(t),......,x_{K}(t)$可以被成功传递，如果${\begin{aligned}&{\text{Allowed Time duration }}\times {\text{ Allowed Bandwidth }}\geq \sum _{k=1}^{K}A_{k}\\&A_{k}{\text{:the area of the time-frequency distribution of }}x_{k}(t)\\\end{aligned}}$

电磁波的传递

应用时频分析的观念，我们可以将一个电磁波表示成一个2x1的矩阵${\begin{bmatrix}x\\y\end{bmatrix}}$。而当电磁波经过一段自由空间时，著名的 菲涅耳衍射就产生了。菲涅耳衍射可以用线性完整转换的参数矩阵${\begin{bmatrix}a&b\\c&d\end{bmatrix}}={\begin{bmatrix}1&\lambda z\\0&1\end{bmatrix}}$来表达，其中z是电磁波在自由空间中传递的距离，而$\lambda$则是波长。

光学

光也是一种电磁波，所以在光学上的应用就跟电磁波传递很类似。

如果电磁波通过一片球面透镜片或是经过一个碟型面的反射，则线性完整转换的参数矩阵可分别表示为${\begin{bmatrix}a&b\\c&d\end{bmatrix}}={\begin{bmatrix}1&0\\{\frac {-1}{\lambda f}}&1\end{bmatrix}}$ 和 ${\begin{bmatrix}a&b\\c&d\end{bmatrix}}={\begin{bmatrix}1&0\\{\frac {1}{\lambda R}}&1\end{bmatrix}}$ ，其中f是球面透镜的焦距，而R是碟型面的半径。

用 LCT 来分析光学系统的好处是只需要用到2x2的矩阵运算 ${\begin{bmatrix}a&b\\c&d\end{bmatrix}}{\begin{bmatrix}x\\y\end{bmatrix}}$，避免了复杂的物理理论和数学积分，但只有在“近轴”的情形下才能准确得到。

光学也可以用在无线通讯(wireless communication), 光学系统分析(optical system analysis), 激光(laser)和雷达系统分析(radar system analysis)。

• 光在真空中传输 (菲涅尔转换，LCT的特例): 啁啾卷积(chirp convolution)
• 光在透镜或碟盘传输: 啁啾相乘(chirp multiplication)

信号鉴别

以下两个信号无法经由单纯的傅里叶分析分辨出来，它们的频谱都长得一样。
$x_{1}(t)={\begin{cases}\cos(\pi t);&t<10\\\cos(3\pi t);&10\leq t<20\\\cos(2\pi t);&t>20\end{cases}}$

$x_{2}(t)={\begin{cases}\cos(\pi t);&t<10\\\cos(2\pi t);&10\leq t<20\\\cos(3\pi t);&t>20\end{cases}}$
不过幸亏有时频分布函数，我们可以看出随时间改变之频率的起落，进而鉴别信号。这个想法也可以应用至模式识别。

语音

语音信号的特性就是它的频率随着时间剧烈变化。因为语音信号所涵盖的资讯非常的多，所以相对的计算时间会是很重要的考量。

根据奈奎斯特采样频率及人耳可听的的频率上限约为20000Hz这两个条件，因此语音信号的取样频率需为40000Hz左右。 然而我们对于时频分析的输出在时间轴的分辨率要求往往不会到这么高，加上时频分析出来的结果为原本输入信号维度的两倍，为减少运算时间，我们会降低输出的取样频率，如100Hz。

以短时距傅里叶变换为例

${X}\left({t,f}\right)=\int _{-\infty }^{\infty }{w\left({t-\tau }\right)\cdot }{x}\left({\tau }\right)\,{e^{-j2\pi \,f\tau }}\cdot d\tau$

可改写为

${X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=n-Q}^{n+Q}{w\left({(nS-p){\Delta _{\tau }}}\right){x}\left({p{\Delta _{\tau }}}\right)}{e^{-j2\pi \,mp{\Delta _{\tau }}{\Delta _{f}}}}{\Delta _{\tau }}$

其中${\Delta _{\tau }}$为输入信号的取样间隔， ${\Delta _{t}}$为输出信号的取样间隔

$S={\frac {\Delta _{t}}{\Delta _{\tau }}},\qquad {\Delta _{t}}\neq {\Delta _{\tau }},\qquad B=Q{\Delta _{\tau }}$

$\left|t\right|>B,w(t)\cong 0\qquad {\frac {B}{\Delta _{t}}}=Q$

其中${\Delta _{\tau }}$需满足下列条件

(1)${\Delta _{\tau }}{\Delta _{f}}={\frac {1}{N}}$ N为一整数

(2)$N>2Q+1$

(3)${\frac {1}{2\Omega }}>{\Delta _{\tau }}\qquad |X(t,f)|\approx 0{\mbox{when}}|f|>\Omega$

${X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)=\sum \limits _{p=nS-Q}^{nS+Q}{w\left({(nS-p){\Delta _{\tau }}}\right){x}\left({p{\Delta _{\tau }}}\right)}{e^{-j{\textstyle {{2\pi \,pm} \over N}}}}{\Delta _{\tau }}$

令$q=p-(nS+Q)\to p=(nS-Q)+q$

${X}\left({n{\Delta _{t}},m{\Delta _{f}}}\right)={\Delta _{\tau }}{e^{j{\textstyle {{2\pi \,(Q-nS)m} \over N}}}}\sum \limits _{q=0}^{N-1}{x_{1}\left({q}\right){e^{-j{\textstyle {{2\pi \,qm} \over N}}}}}$

其中 ${\begin{cases}{x_{1}}\left(q\right)=w\left({(Q-q){\Delta _{\tau }}}\right)x\left({(nS-Q+q){\Delta _{\tau }}}\right),&{\mbox{for}}{\rm {0}}\leq q\leq 2{\rm {Q}}\\{x_{1}}\left(q\right)=0,&{\mbox{for}}{\rm {2}}Q\leq q\leq N\end{cases}}$

时频分析在语音领域同时包含音乐信号、声音信号以及声纹辨识。

例如语音信号:

• 声纹:不同人说话的声音频谱(声纹)不同
• 同一个人发不同的声音也拥有不同的声音频谱
• 语调的不同也会使频谱的变化情况不同
• 同一个字音，其中子音与母音的频谱亦不相同
• 双母音也有不同频谱的变化

例如研究指出在声纹(中文注音第一、二、三、四声和轻声)当中的语调会使时间和瞬时频率的关系有所不同。

生医工程

时频分析在生医工程上的应用几乎都是用以分析生理信号，如肌电图(EMG)、心电图(ECG)等等。

肌电图(EMG)信号处理

其中肌电图(EMG)是肌肉收缩时的电位变化和肌纤维震动的变化所产生的生理信号，故常用以探讨肌肉收缩力量程度大小或用来判定肌肉是否产生疲劳的工具，进而推估身体状态的生理信号。EMG信号为非周期性且随机的信号，所以若将时域信号进行快速傅里叶转换(FFT)并不适宜，取而代之的是对EMG信号进行短时傅里叶转换(STFT)，获得EMG信号的功率频谱密度函数(power spectral density function，PSDF)，其反映了EMG信号频率随时间而变化。

短时傅里叶转换属于时频分析的一种，是在傅里叶转换中加入一个移动的视窗函数(window function) w(t)，用来对输入的信号做切割，在对视窗内的信号做傅里叶转换，产生一个二维的时间频率分布图。视窗选择方面大约包含了Hamming、Hanning、Gaussian 等等，而视窗选择对于STFT有很大的影响，若为较小的window function可得到较佳的时域分辨率，但会牺牲频域分辨率；反之，若选择较大的window function将会得到较佳的频域分辨率但较差的时域分辨率，借由适当window size的设定，就能观察出肌肉收缩力量程度大小随着时间而改变。

HRV analysis of ECG signal in frequency domain

心电图(ECG)信号处理

心电图也是常见的电生理信号之一，记录了心脏整个活动过程，心脏收缩前必先产生电气活动(electrical activity)，此电气活动是由许多心肌的兴奋波所组成。兴奋波起源于心脏的窦房节(SA Node)，并经由心肌特化的传导系统(conducting system)将此兴奋波传遍整个心脏，心电图即是源于心脏的兴奋区与未兴奋区之间的电位差。而在心电图(ECG)的应用最常见的例子即为心率变异度分析（heart rate variability, HRV），将ECG信号使用时频分析可用以检测自主神经系统活性，亦可进行进行个人压力与情绪分析。

为了评估自主神经系统功能及对心血管活动的影响，亦可对ECG信号进行希尔伯特-黄变换（Hilbert-Huang transform,HHT）以获得Hilbert时频谱，利用Hilbert-Huang时频来做时频特征提取和分析。依据短时程HRV信号的线性频域分析指标，得到不同生理频带的Hilbert能量图，提取总能量，各生理频带的能量和其归一化能量以及生理频段的能量比值作为评价心率变异性的时频特征。基于Hilbert谱的时频特征的区分性能好，有较清晰的生理意义，能反映人的生理病理变化，为短时程HRV信号分析提供了一种有效方法。

以往大多数的心率变异度分析都是以快速傅里叶转换频谱分析为基础，来探讨研究现象之心率变异频谱特性与自主神经活性之间关联性。但是对于快速变化的生理信号，如ECG是一个非稳定(non-stationary)信号，使得傅里叶转换无法达到有效的处理，使用快速傅里叶转换频谱分析有其局限性。另一种较为合适的方法是以平滑式魏格纳－韦立分布时频分析法(smoothed pseudo Wigner-Ville distribution; SPWVD)为基础，分析受试者之心率变异度之频谱特性，来检视受试者心率变异度中之低频频谱能量(LF)与高频频谱能量(HF)的增减现象，其比值为自主神经平衡指标(LF/HF)，其频率范围定义如下所列。

1. 极低频范围的正常心跳间期的变异数(VLF) - 截取频率范围为0 - 0.04 Hz

2. 低频范围的正常心跳间期的变异数(LF) - 截取频率范围为0.04 - 0.15 Hz，为交感神经活性指标

3. 高频范围的正常心跳间期的变异数(HF) - 截取频率范围为0.15 - 0.4 Hz，为副交感神经活性指标

心音(PCG)信号处理

PCG信号含有预知身体状况的资讯，所以可说是有临床诊断的价值。

在一个心动周期当中，有两个主要的声音，也就是第一心音及第二心音。

第一心音是由四个成分组成，是因为主动脉膜的闭合所产生。

第二心音是由两个成分组成，是因为肺动脉瓣膜的闭合所产生。

虽然心音信号可以说是一个不段重复的信号，但是我们不将它定义为一个简单的周期信号。

这是因为，心音信号具有非平稳性，而且心音信号是由许多复杂的事件所产生的信号所叠加的多成分信号。

而传统的频谱分析方法，只能反映信号的静态频谱特征，但是，能够反映心脏及心血管系统病理特征的心音信号，一般来说是非平稳时变特征。

所以，使用频谱分析方法，就无法反映心音信号随时间的变化。

因此，心音信号的时频分析用于心脏方面的研究和临床诊断，有非常好的分析效果。

脑电图(EEG)信号处理

EEG分析，是一种研究大脑电信号的方法，并且是一种无创伤的方法。

EEG可以分成三种，也就是非瞬时自发性EEG、瞬时自发性EEG和诱发性EEG。

傅里叶转换和各种衍生的方法，可以很好地用在非瞬时自发性EEG。

但是，对于瞬时自发性EEG和诱发性EEG，用傅里叶转换的方式来分析显然无法给予充分的资讯，因为这两种信号都是非平稳信号，因此就会用到时频分析的方法。

另外，EEG信号中，有时会出现一些短瞬时脉冲，这些脉冲可能含有病理资讯，也可能仅仅只是干扰信号。不论属于哪种类型，对这些脉冲信号做检测和分析都十分重要，因此将时频分析方法应用在EEG信号处理中，是许多学者的研究方向。

机械设备故障诊断

一个机械设备的故障诊断可以概括为以下五个部分：

(1) 采集信号

(2) 从采集的信号中，用信号处理的方法，提取出能够反映机械设备状态的特征。

(3) 以某些状态识别方法，监测机械设备并判断机械设备的状态，检查是否有故障的状况。

(4) 分析与诊断机械设备的状态。若机械设备发生故障，要分析故障的类型、性质、发生部分等。

(5) 最后，根据机械设备的状态及可能的发展趋势做出决策，例如维修。

而对于搜集来的信号，不仅只会搜集到能够反映出机械设备某特定部分的工作状态资讯，也会搜集到其他零件的资料，这些资料对于研究该特定部分的状态而言，是一种背景噪声。

另外，一般而言，这样的背景噪声会比状态资讯来的更大，而由于大多数的机械信号都是非平稳信号，特征时频谱范围又广，成分非常复杂，并且又如前述参杂许多背景噪声。

所以，要怎样把所需的状态资讯提取出来，提高讯杂比，并且将故障特征资讯放大，是一个重要的机械故障诊断的研究方向。

时频分析则为机械故障诊断带来了一线曙光。作为一种联合时间与频率的分析方法，它能够很好的分析非平稳信号，并且将信号与噪声做分离，找出所谓的故障信号。

目前，几种主要的时频分析方法已被提出并应用于机械设备故障诊断，主要包括短时傅里叶变换（STFT）、小波变换（WT）和希尔伯特-黄变换（HHT）。

滚动轴承故障诊断

滚动轴承是机械设备中至关重要的部件，其运行状态直接影响整个机械系统的可靠性和效率。滚动轴承故障通常以振动信号中的周期性冲击表现出来，这些冲击具有强烈的非平稳性。因此，传统的振动分析方法如时域分析、频域分析往往无法有效捕捉这类故障的特征。为了克服这一挑战，时频分析方法如小波变换（Wavelet Transform, WT）在滚动轴承故障诊断中得到了广泛应用。

小波变换具有多分辨率分析的特性，能够同时处理信号的高频和低频成分，并且能在时间和频率上自适应地选择分解层级。因此，它能够有效捕捉到滚动轴承故障引起的瞬时冲击信号。此外，结合深度学习技术，这些时频特征可以进一步用于故障分类。比如，Zhang等人（2020年）提出了一种基于时频分析和深度学习的滚动轴承故障诊断方法。该方法首先利用连续小波变换方法提取滚动轴承的故障特征，然后将这些特征输入深度学习模型进行故障分类，结果显示该方法能够显著提高故障诊断的准确性，尤其是在处理复杂故障模式时。

齿轮故障诊断

齿轮故障是机械设备中常见的问题，故障会引起机械系统的振动信号变化。由于齿轮故障引起的振动信号通常较为微弱，传统的频域方法难以有效检测。为了更准确地识别这些微弱信号，Liu等人（2022年）提出了一种基于连续小波变换（CWT）的齿轮箱故障诊断方法。该方法通过CWT对齿轮箱故障振动信号进行时频分析，构建其时频图，捕捉到信号在时间和频率域的变化特征。CWT能够有效地揭示齿轮故障引起的微弱振动，并提供更为灵敏的故障识别能力。这一方法克服了传统频域方法的局限，实现了对早期齿轮故障的敏感检测，并显著提高了诊断准确性

引擎故障诊断

引擎是机械系统中一个至关重要的部件，其运行状态直接影响整体性能和效率。在可变气门正时系统（VVT）控制的引擎中，机械部件的故障可能会表现为特定的声音特征，这些声音随着引擎运行状态的变化而改变。因此，声音信号的分析成为引擎故障诊断的一个重要方法。传统的诊断方法通常依赖于振动或其他物理量的测量，而利用时频分析方法来识别引擎中各个机械部件的故障则能提供更多的诊断信息。

时频分析方法，特别是希尔伯特黄变换（HHT），在这一领域中发挥了重要作用。HHT能够有效地处理非平稳信号，并且在时间和频率域中自适应地分析信号的特征。通过对引擎声音信号进行HHT分析，可以提取出不同机械部件的声音能量分布，从而识别哪些部件处于良好工作状态，哪些部件可能存在故障。例如，研究中使用HHT分析引擎声音信号，发现引擎的声音能量分布与引擎输出扭矩之间存在关联，这一发现有助于早期识别故障并提高故障诊断的准确性。

近岸高频测流雷达

运用无线电受到海面反射来的回波频谱来估算表层海流。原理包括了几种物理现象，布拉格散射、多普勒效应、深水波假设。

海洋的表面是粗糙的，其中包含了各种波长的波浪，当一个近岸的测站发射无线电波，电波碰撞到海面，因为布拉格散射，波浪波长为二分之一的无线电波长的波浪会造成强的回波。

若是打向一静止的粗糙面，回波信号的频率应该与发射频率相符，但海面上的波浪是行进的，对于电波接收站而言，其所收到的信号，是一个移动中的信号来源，所以观察到的回波频率因为多普勒效应，回波信号频率不再是原发射频率，而会发生在发射频率加上一偏移频率。

因为波浪没有特定的方向，可说是在海面往四面传递，对于无线电侧站位置而言，远离的波会造成负的偏移频率，即回波频率降低，反之，靠近测站的波浪会形成一偏移频率为正的回波。

其偏移频率的大小值与波速相关，波浪波速根据深水波的假设，波速是波长之函数，因为先前布拉格散射，当无线电波频率给定，其回波信号主要是二分之一无线电波长的波浪，故此偏移频率可以估计出来。

由于表层海流载着海表面的波浪走，所以观测到的无线电频谱，和预测上的频谱会有所不同，而这中间不同的量可以来估算出表层海流。

通常无线电波站都有Ｘ, Y, Z轴三个方向的天线，借由三个方向的接收信号，雷达可以分辨出回波的方位，借由到达回波在时间轴上可以分辨出回波的距离。

以下举一个虚构的例子来解说。

An trivial example of received signal spectrum

无线电测站发出一频率为Fc，也就是波长为${\frac {c}{Fc}}$的无线电波，因为1) 海表面有波长不一的波浪2) 因为布拉格散射，波长为${\frac {c}{2Fc}}$的波浪会产生较大的回波能量 3)这些波浪四面传递，且波浪波速在深水波假设中，可以近似成$c_{w}ave={\sqrt {\frac {g\lambda }{2\pi }}}$。

回波的产生源，波浪，因为以波浪波速移动，多普勒现象的缘故，回波信号会有一个频率位移Δf，远离测站靠近测站的波浪分别造成Fc-Δf， Fc+Δf的回波（如虚线所示）。

假设又一个靠岸（流向测站）表层海流，海流带着波浪行进，其结果是离岸的浪减速，靠岸的流加速，其在频率域相当于右偏（频率增加）。借着观察回波频率极强值（实线）和无流情况推算出来的频率强值之间的差距，来推算出迳向上流速（远离或靠近测站方向上的流速）。

鲸豚研究

鲸豚的研究常常仰赖着研究者现场目视观察，由于人力物力有限，只能在小范围内，海况佳，日间的情况下来研究。鲸豚的一大习性是发声，透过被动地监听鲸豚，研究者可以更进一步的来研究。以海豚的声音为例，最粗浅地看至少有两种，即click 和 whistle。二种声音从时频图来看是一目了然，二种在时频图有截然不同的特性，click是短时间(50-150 μs)宽频的信号(65 – 100kHz)，whistle是频率调变的信号（2 – 20kHz）。一般来说，click信号和海豚用来定位导航，猎食有关，而whistle和海豚的社交，沟通有关。也有可能同一时间发出两种声音。这当中许多学问，包含各种声音的功能，发声的机制，有待厘清。

资料压缩

资料的压缩包含图像、影像及语音的压缩，其中最广为利用的时频分析方法应为小波分析。以影像来说，是用二维离散的小波转换进行压缩，图一为二维离散小波转换的结构图，所得的四个结果分别为影像中的低频及不同范围的高频成分，$x_{1,L}[m,n]$为影像中的低频成分，$x_{1,H1}[m,n]$为水平方向的边缘，$x_{1,H2}[m,n]$为垂直方向的边缘，$x_{1,H3}[m,n]$为图像中的角落。而保留低频成分$x_{1,L}[m,n]$，舍弃其他高频部分，再进行数次的二维离散的小波转换，可得粗略但很接近原图的缩图，以达到缩图的效果，每进行一次二维离散的小波转换可将资料量缩为原图的1/4倍左右，其中小波分析的影像压缩以JPEG 2000为主。

图一

卫星信号

台湾的第二颗卫星FORMOSAT-2(FS-2)是具高分辨率的遥感卫星，而FORMOSAT-3(FS-3)/COSMIC(Formosa 3号卫星，气象、电离层和气候星座观测系统）包括6个低地球轨道(Low-Earth-Orbit, LEO)卫星是第一个演示近实时数值天气的星座，使用来自全球定位系统(Global Positioning System, GPS)卫星的无线电信号进行预测(Numerical Weather Prediction, NWP)。资料中指出，卫星任务每天都会受到自动重配置命令 (Automatic Reconfiguration Order, ARO)的干扰，FS-2记录着相当多ARO事件，另外在另一颗FS-3上也记录了许多计算机异常事件(卫星重启或重置)，同时也有相关纪录在FS-2上。而这些ARO事件也大多归因于发现这些异常当中有相当star tracker data多应归因于单事件失败(Single Event Upset, SEU)，也因此我们必须了解更多SEU的起因来避免ARO的事件一再发生而使卫星受到干扰。在研究中发现FS-2的ARO和星体追踪仪数据的遗失位元(Lost Bytes, LB)与太空天气成正向相关，而几个太空气象的主要参数正是影响卫星任务的原因:例如地球磁场的Kp指数，质子密度，电子密度和10.7 cm radio flux (RF)。另外也发现FS-3的电脑重置也与地球磁场的Kp指数，质子密度，10.7 cm radio flux (RF)和X射线相关。因此若能使用时频分析当中的希尔伯特-黄转换和其他方法来分析此非线性且不稳定的太空数据，可以结合这些太空天气的相关参数，进而避免SEU甚至其余ARO事件的发生几率，也可使得卫星任务更为稳定。

加速规分析

加速规在各种应用中扮演着重要角色，尤其是在动态测量、健康监测等。透过时频分析，我们可以深入了解加速规数据的时间变化及其频谱特征，特别适用于处理非稳定信号和动态事件。以下将介绍几个加速规数据利用时频分析的应用范畴。

时频分析的核心在于能够同时观察信号在时间域与频域上的特征，这在处理如加速规这类动态信号时格外有用。例如，当我们面对加速规数据时，往往需要了解信号中包含的频率成分及其随时间的变化。这正是时频分析所提供的价值，它能够在时间与频率上提供综合视角。

1. 跌倒检测

在跌倒检测领域，加速规被广泛用于监测使用者的活动，并区分跌倒与其他日常活动。透过时频分析，可以有效地将跌倒事件与非跌倒事件区分开来。例如，研究显示，跌倒事件在特定时间窗内会出现高能量的频率成分，这些频率通常集中在 2 到 3.5 Hz 之间，而跑步、跳跃等活动则表现出不同的频率范围。经过时频图的分析，能够实现高灵敏度和高特异性，显示出时频分析在动态信号识别中的强大能力。

2. 人体健康监测

对于患有帕金森氏症的人，加速度计数据的时频分析可以用来检测如“冻结步态”（freezing of gait）等症状，利用加速规数据进行步态冻结的实时检测，并使用时频分析来分析加速规信号的特征。短时傅里叶变换（STFT）来进行时频分析，这使得他们能够捕捉到步态运动中的频率变化，尤其是在帕金森病患者的步态冻结事件发生时。

时频分析能够将加速规信号的时间变化和频率成分同时显示出来，帮助识别信号中的异常模式。例如，当患者的步态发生冻结时，步态运动的频率成分会显著改变，这一变化可以通过时频分析中的频率图清晰地观察到。该方法提供了精确的时间-频率定位，能够及时检测到步态冻结，进而实现对帕金森病患者的即时预警。

经济资料分析

传统上，经济数据被视为不稳定且具嘈噪声号的时间序列资讯，且被正统理论认为是随机程序，基本上不存在规律性。经济或金融资料属于时变的资讯，而其在此领域中大多使用统计方法做分析，例如自回归模型(AR)，经济学家大多使用经验数据来手动拟合人工线性模型，而此随机模型的拟合通常只能计算均值和方差，对于大量的金融经济应用层面并不足够。因为经济金融的时序统计资料在频域上是相对未开发的领域。在这种情况下，其频谱资讯会随着时间而一直改变，以传统的傅里叶转换是不足以充分描述该资讯的循环特性，因此后人研发了联合时频表示法(joint time-frequency representation, JTFR)来克服这个问题，因为它能同时在时域和频域表现其时间序列和做分析。使用JTFA，我们不仅可以知道此信号(资讯)当中存在哪种类型的周期，而且知道它们何时发生以及持续多长时间。当有持续一定时间的周期集中在联合时频域中时，随机噪声趋于均匀地扩展到整个时频平面。因此，JTFA具有更高的信噪比(SNR)。例如部分研究中，会在JTFA当中先用加伯转换(Gabor transform)和短时傅里叶转换两个线性模型分析原始的ISE资料，以及两个二次模型(韦格纳分布和Page distribution)做时频分析。通过特殊合成的时间序列可以借此评估每个TFR在检测和解码原始ISE数据中可能存在的特性，也可以对它们的趋势或是周期分量进行分析和计算。此方法可以大致模拟股票指数系列(例如原始ISE系列)的模式，并以此做为比较绩效分析的基准。接着比较不同时频分析方法以取最优，并再使用过去指数数据重复做分析，以察看结果是否对从新兴市场到成熟市场的股票数据来源变化是否可靠。另外也可以使用希尔伯特-黄转换，不是先进行预白化，而是找到一条平滑的趋势曲线以拟合经验数据，以便差值包含尽可能多的有意义的周期;尽管平滑曲线呈现长期趋势，但差值可进一步用于分析短期行为。尽管频谱表明了残差所包含的频率成分，但它并不能告诉它们何时出现或持续多长时间。因为只有那些持续一定时间的频率才被认为是有意义的周期，所以无法仅从频谱上确定经济周期。

使用JTFA或其余时频分析方法进行经济数据分析的发展是指日可待的，这种发展趋势可以使金融界发生分析及预测上大革命。

地震波分析

以信号处理的方式解析地震波有助于地质学家更准确的掌握地底构造，透过观察地震波在地底中传播的情形，辨识岩层、矿物等物质。除了自然发生的地震，人为爆破而生成的震波，搭配信号处理技术的使用，也可以应用于石油探测及矿石开采等实务面。

近代随着科技进步，电脑辅助运算成为必要之工具，地震观测主要为数位地震仪系统，这样一来，可获取大量取样点及高分辨率、高带宽、足够的动态范围的数码化地震波形资料。由于地震观测系统得动态范围和频率范围皆为有限值，并且记录地面的活动有着不同程度的变形；因此不同型态的外界干扰和系统内部噪声会对于地震波纪录造成干扰。在资料记录过程中，会出现倾斜、平移、波形变形等现象，因此在对地震波进行时频分析前，必须进行修正，调整仪器内部影响，以及过滤噪声。

以时频分析而言，对于地震波的处理，可以分为以下几种，分别为使用短时傅里叶转换(Short Time Fourier Transform), 希尔伯特-黄转换(Hilbert-Huang Transform)以及小波转换(Wavelet Transform)。

短时傅里叶转换(Short Time Fourier Transform, STFT)

$X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }\,d\tau$

窗函数(window function)若为矩形，则写为 $w(t)=1,t<|T|$， $o.w.,w(t)=0$

计算速度快，呈现的图形直观且方便，可以于野外观测现场进行。

分析过程中，主要目的为确定各窗函数(window function)中的主频，因此窗函数大都选择为矩形，矩形窗(rectangular window)的宽度 T 一般选择大于地震波的周期，透过窗函数的滑动，得到地震波随时间变化的频率表现。

希尔伯特-黄转换(Hilbert-Huang Transform, HHT)

适合用来对非线性及非稳定的信号做时频分析。

透过经验模态分解及希尔伯特转换，得到包含时间、频率、及振幅三种资讯的结果。由于经验模态分解的效率高，加上此时频分析方法不依赖于傅里叶转换的使用，因此能有效反映信号内部特征，达到去除噪声并保有信号的非线性及非稳定特征，符合地震波的信号分析。

小波转换(Wavelet Transform)

小波转换分析具有良好的时频局部化性质，优于传统的傅里叶分析。

透过时间窗(time window)及频率窗(frequency window)的改变，能在时间轴及频率轴上表达局部特征，因此，适合用于观察一正常信号中的反常结果，在分析信号的过程中具有放大镜的功能。

其它时频分析方法

另外，ZAM(Zhao-Atlas-Marks distribution )方法为非线性时频分析方法，其母函数为锥形函数，研究显示，此类时频分析方法适合应用像地震波等非平稳信号之研究，能有效抑制时频分析结果中的交互干扰项，同时使时频分析之结果密集度较高。

以下举例以时频分析对地震波之分析结果:

对自然发生之地震及人为爆破事件进行ZAM时频分析，观察结果得知

1.包括天然地震和人为爆破在内的数位地震波形中，S波之较高能量密度区域处于低频区域，其频带比P波之高能量密度区域频带来得窄，且其能量密度比P波能量密度更大；

2.天然地震的高能量密度分布较为离散，而人为爆破的高能量密度分布区域较为集中。

延伸阅读

• Jian-Jiun Ding, Time frequency analysis and wavelet transform class note, the Department of Electrical Engineering, National Taiwan University (NTU), Taipei, Taiwan, 2018.
• 葛哲学, and 陈仲生. "Matlab 时频分析技术及其应用." 人民邮电出版社, pp10-15 (2006).