语音信号处理 | 傅里叶变换、短时傅里叶变换、小波变换、希尔伯特变换、希尔伯特黄变换

在信号处理领域，存在诸多变换，比如标题中的五个变换。本文将对这五个变换进行介绍和比较。在开始之前，我们需要先理清什么是平稳信号，什么是非平稳信号。

我们知道，自然界中几乎所有信号都是非平稳信号，比如我们的语音信号就是典型的非平稳信号。那么何谓平稳信号和非平稳信号呢？一个通俗的理解即，平稳信号在不同时间得到的采样值的统计特性(比如期望、方差等)是相同的，非平稳信号则与之相反，其特性会随时间变化。在信号处理中，这个特性通常指频率。

通常傅里叶变换只适合处理平稳信号，对于非平稳信号，由于频率特性会随时间变化，为了捕获这一时变特性，我们需要对信号进行时频分析，就包括短时傅里叶变换、小波变换、希尔伯特变换、希尔伯特黄变换这几种变换。以下逐一进行分析介绍。

傅里叶变换(Fourier Transform, FT)

首先考虑一个连续信号f(t)f(t)f(t)的傅里叶变换和它的反变换，如下
F(ω)=F[f(t)]=∫−∞+∞f(t)e−jωtdtF(\omega)=F[f(t)]=\int_{-\infty}^{+\infty}f(t)e^{-j\omega t}dt F(ω)=F[f(t)]=∫−∞+∞f(t)e−jωtdt

f(t)=F−1[F(ω)]=12π∫−∞+∞F(ω)ejωtdωf(t)=F^{-1}[F(\omega)]=\frac{1}{2\pi}\int_{-\infty}^{+\infty}F(\omega)e^{j\omega t}d\omega f(t)=F−1[F(ω)]=2π1∫−∞+∞F(ω)ejωtdω

在实际应用中，计算机只能处理离散信号，所以对连续信号x(t)x(t)x(t)进行时域采样，得到一组离散样本x(n)x(n)x(n)，对它进行傅里叶变换得到
X(ω)=∑n=−∞∞x(n)e−jωnX(\omega)=\sum\limits_{n=-\infty}^\infty x(n)e^{-j\omega n} X(ω)=n=−∞∑∞x(n)e−jωn
上式即为离散时间傅里叶变换(DTFT)，由于变换后得到的频域值仍然是连续的，我们继续对频域进行采样，得到
X(k)=∑n=0N−1x(n)e−j2πknNX(k)=\sum\limits_{n=0}^{N-1}x(n)e^{-j\frac{2πkn}{N}} X(k)=n=0∑N−1x(n)e−jN2πkn
上式就是离散傅里叶变换(DFT)，当前计算机中常用的快速傅里叶变换(FFT)，就是DFT的快速算法。如下即为信号x(t)=cos(200πt)+0.5cos(400πt)x(t)=cos(200\pi t)+0.5cos(400\pi t)x(t)=cos(200πt)+0.5cos(400πt)的采样信号和经过FFT得到的幅度谱，其中采样频率fs=500Hzf_s=500Hzfs=500Hz

可以看到，该信号由两个频率分别为100Hz和200Hz的分量构成。利用频谱分析，我们可以发现一些时域中看不到的信息，比如信号的频率分量组成、信号能量在频域上的分布等。

然而，傅里叶变换是一种全局的变换，时域信号经过傅里叶变换后，就变成了频域信号，从频域是无法看到时域信息的，我们可以从上节中的傅里叶变换和反变换公式进行解释，进行正变换时，积分区间为整个时域，所以变换结果将不包含时域信息，反变换同理。

如果信号的频率特性在任何时间都不发生改变(即该信号是平稳信号)的话，使用傅里叶变换是没有问题的，然而如果该信号是非平稳信号，这时候时域信息就相当重要了。举个栗子，以下分别是0_{100Hz线性递增扫频信号和100}0Hz线性递减扫频信号的幅度谱，这两个信号都是非平稳信号，可以看到它们的幅度谱是相同的，很显然我们无法知道每一个频率分量出现的时间，即无法知道扫频信号的模式。

上述例子可能比较简单，可以直接从时域信号看出扫频模式，但是对于一些复杂的信号和它的幅度谱，比如以下的两路扫频信号叠加后的信号和它的幅度谱

仅仅通过时域信号或者幅度谱，我们是很难分析这段非平稳信号的特征的。下面我们将引入一个时频分析( Time-Frequency Analysis)方法——短时傅里叶变换(STFT)。

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

短时傅里叶变换定义为
X(n,ω)=∑m=−∞∞x(m)w(n−m)e−jωmX(n,\omega)=\sum_{m=-\infty}^\infty x(m)w(n-m)e^{-j\omega m} X(n,ω)=m=−∞∑∞x(m)w(n−m)e−jωm
其中x(m)x(m)x(m)为输入信号，w(m)w(m)w(m)是窗函数，它在时间上反转并且有n个样本的偏移量。X(n,ω)X(n,\omega)X(n,ω)是时间nnn和频率ω\omegaω的二维函数，它将信号的时域和频域联系起来，我们可以据此对信号进行时频分析，比如S(n,ω)=∣X(n,ω)∣2S(n,\omega)=|X(n,\omega)|^2S(n,ω)=∣X(n,ω)∣2就是语音信号所谓的语谱图(Spectrogram)。

画出上节中两路扫频信号叠加后的信号的语谱图，如下图

可见该信号是由一个0 ~ 250Hz二次递增的扫频信号和一个250 ~ 0Hz二次递减的扫频信号的叠加。通过STFT，我们可以很容易地得出非平稳信号的时变特性。

继续对STFT进行分析。

计算语谱S(n,ω)S(n,\omega)S(n,ω)时采用不同窗长度，可以得到两种语谱图，即窄带和宽带语谱图。长时窗（至少两个基音周期）常被用于计算窄带语谱图，短窗则用于计算宽带语谱图。窄带语谱图具有较高的频率分辨率和较低的时间分辨率，良好的频率分辨率可以让语音的每个谐波分量更容易被辨别，在语谱图上显示为水平条纹。相反宽带语谱图具有较高的时间分辨率和较低的频率分辨率，低频率分辨率只能得到谱包络，良好的时间分辨率适合用于分析和检验英语语音的发音。

如下图所示，分别为一段语音的帧长为128和512的语谱图。

可见，对于帧长固定的短时傅里叶变换，在全局范围内的时间分辨率和频率分辨率是固定的。如果我们想要在低频区域具有高频率分辨率，在高频区域具有高时间分辨率，显然STFT是不能满足要求的。我们继续引入另一种时频分析方法——小波变换。

小波变换(Wavelet Transform, WT)

对于任意能量有限信号f(t)f(t)f(t)，其连续小波变换(CWT)定义为
Wf(a,b)=1a∫−∞+∞f(t)ψ∗(t−ba)dtW_f(a,b)=\frac{1}{\sqrt a}\int_{-\infty}^{+\infty}f(t)\psi^*(\frac{t-b}{a})dt Wf(a,b)=a1∫−∞+∞f(t)ψ∗(at−b)dt
其中ψ(t)\psi(t)ψ(t)是母小波或者基本小波，它满足ψ(±∞)=0\psi(±\infty)=0ψ(±∞)=0，ψ(0)=0\psi(0)=0ψ(0)=0，∫−∞+∞ψ(t)dt=0\int_{-\infty}^{+\infty}\psi(t)dt=0∫−∞+∞ψ(t)dt=0，前两个条件表明ψ(t)\psi(t)ψ(t)在时域上是一个有限长的函数，最后一个条件则表明ψ(t)\psi(t)ψ(t)必须时正时负地波动，否则它的积分结果不会为零，因此它在频域上也是有限的。所以不同于傅里叶变换的基函数是一个无限长的正弦波，小波变换的基函数是一个经过衰减处理的有限长小波，小波基函数在时域和频域上都是局部化的，如下图分别为傅里叶变换和小波变换的基函数

将ψ(t)\psi(t)ψ(t)进行伸缩和平移得到一族函数ψa,b(t)\psi_{a,b}(t)ψa,b(t)，称为分析小波，这就是小波变换的基函数族，其中aaa为伸缩参数，当a>1a>1a>1时，沿时间轴方向拉伸，因子1/a1/\sqrt a1/a是为了保持伸缩之后能量不变；bbb为平移参数。

使用MATLAB的小波变换工具箱画出上节两路扫频信号叠加信号的小波变换结果，如下图，其中纵轴(频率轴)是对数轴。

可见在低频区域的变换结果具有较高的频率分辨率(频率轴是对数轴，在低频区域跨度较小)，在高频区域具有较高的时间分辨率。

一些总结

我们对上述方法进行总结，可以发现上述方法的区别都可以归结为选取的时频窗尺寸不同。画出时域信号、频域信号、使用STFT得到的时频域信号、使用小波变换得到的时频域信号的时频窗，如下图，每一个小方块表示一个时频窗，沿时间方向的边长表示时间分辨率，沿频率方向的边长表示频率分辨率

由上图可以知道：

对于时域信号，它可以有很高的时间分辨率，然而其频率分辨率为零。
经过傅里叶变换得到的频域信号可以实现很高的频率分辨率，然而其时间分辨率为零。
对于短时傅里叶变换(STFT)，它在时域和频域都有一定的分辨率，并且在全局范围内STFT的时频分辨率都是一样的。但是由于Heisenberg不确定原理(也就是量子力学中的测不准原理)的制约，每一个时频窗的面积都是固定的，即时间分辨率和频率分辨率成反比，所以这两个分辨率不能同时很高。
小波变换在不同时间和频率上具有不同尺寸的时频窗，可以在低频区域实现较高的频率分辨率，然而其仍然受到Heisenberg不确定原理的限制，时间分辨率和频率分辨率不能两全其美。同时小波变换的时频窗并非完全是自适应的，它还需要人为地选择基函数。

上述的方法都会受到Heisenberg不确定原理的限制，而且并不是完全自适应的方法。接下来介绍一种不受Heisenberg不确定原理限制、同时还有更好的自适应性的时频分析方法——希尔伯特黄变换

希尔伯特变换(Hilbert Transform, HT)

在介绍希尔伯特黄变换之前，我们先介绍一下希尔伯特变换。

希尔伯特变换也是傅里叶变换的一种扩展，它常常用于通信系统中的调制解调，当然它也可以用于信号的时频分析。其计算方法为

计算输入信号的FFT，保存为向量F
创建一个向量h，其中
h(i)={1,i=1,n2+12,i=2,3,⋯,n20,i=n2+2,⋯,nh(i)=\begin{cases}1&,i=1,\frac{n}{2}+1 \\2&,i=2,3,\cdots,\frac{n}{2}\\0&,i=\frac{n}{2}+2,\cdots,n\end{cases} h(i)=⎩⎪⎨⎪⎧120,i=1,2n+1,i=2,3,⋯,2n,i=2n+2,⋯,n
计算F与h的内积
计算上步得到的序列的iFFT

稍微介绍下它在通信系统中的应用，利用MATLAB中的hilbert函数对一个扫频信号进行希尔伯特变换，得到的结果是一个解析信号，在同一坐标系画出该解析信号的实数部分和虚数部分，如下图

可见虚数部分较实数部分滞后π/2\pi / 2π/2，我们可以利用这个性质消掉负频率。使用MATLAB中的fft函数对一个信号进行傅里叶变换，结果中的负频率保存在下半轴。如下是一个扫频信号的幅度谱以及使用希尔伯特变换得到的解析信号的幅度谱

由上图可见，解析信号的幅度谱没有负频率，并且各个频率分量的幅度是原来实信号的两倍。

通过希尔伯特变换可以得到没有负频率的解析信号，基于此可以实现信号的单边带调制，从而节省带宽和降低发射功率。感兴趣可以继续查看使用希尔伯特变换进行单边带幅度调制

继续探讨基于希尔伯特变换的时频分析。在时频分析领域，希尔伯特变换主要用于瞬时频率估计。

由上述分析可知使用希尔伯特变换可以得到原始信号的解析信号，假设解析信号为z(t)z(t)z(t)
z(t)=c(t)+jy(t)=aejθ(t)z(t)=c(t)+jy(t)=ae^{j\theta(t)} z(t)=c(t)+jy(t)=aejθ(t)
其中a(t)=c(t)2+y(t)2a(t)=\sqrt{c(t)^2+y(t)^2}a(t)=c(t)2+y(t)2，表示瞬时幅值，θ(t)=arctany(t)c(t)\theta(t)=arctan\frac{y(t)}{c(t)}θ(t)=arctanc(t)y(t)，表示瞬时相位，ω=dθ(t)dt\omega=\frac{d\theta(t)}{dt}ω=dtdθ(t)，表示瞬时频率。由瞬时幅值和瞬时频率可将信号表示为
x(t)=a(t)ej∫ω(t)dtx(t)=a(t)e^{j\int\omega(t)dt} x(t)=a(t)ej∫ω(t)dt
若使用∣a(t)∣2|a(t)|^2∣a(t)∣2表示瞬时能量，则可在时间-频率面上画出信号的瞬时能量分布，这个分布谱图就是Hilbert谱，记为H(ω,t)H(\omega,t)H(ω,t)。如下为一个0~100Hz扫频信号的时域信号和它的Hilbert谱

上图中的信号为单频率成分信号，即同一时刻只有一个频率分量的信号，我们可以由Hilbert谱很好地观察出信号的时频特征，且有很高的的时间分辨率，但是信号边界处的误差往往较大。

如果是对一个多频率成分信号(同一时刻有多个频率分量的信号)进行希尔伯特变换，结果会怎样呢？如下为两路扫频信号叠加信号的Hilbert谱

从上图中的Hilbert谱我们并不能概括出该信号的时频特征，所以对多频率成分信号不能直接进行Hilbert变换，我们还需要对其进行进一步处理，将原始信号分解成单频率信号的叠加，这就要用到希尔伯特黄变换中的EMD分解。

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

相比于HT，HHT就多了一个经验模态分解(Empirical Mode Decomposition, EMD)，EMD就是把复杂信号分解成从高频到低频的若干个固有模态函数(Intrinsic Mode Function, IMF)，IMF需要满足两个条件：

信号极值点的数量与零点数相等或相差为1
信号的由极大值定义的上包络和由极小值定义的下包络的局部均值为0(即包络上下对称)

简单的理解就是，EMD是依次提取信号在每个局部的最高频分量的过程，所以每个IMF实际上是一个单频率分量信号，这样我们就可以对每个IMF分量进行Hilbert变换，从而得到每个分量的Hilbert谱。如下是对两个扫频信号叠加后的信号进行EMD分解得到的IMF分量

对上述的IMF分量进行Hilbert变换，求得Hilbert谱，如下

从上图我们大致可以看出信号的时频特性：该信号是一个0 ~ 250Hz二次递增的扫频信号和一个250 ~ 0Hz二次递减的扫频信号的叠加。

当然HHT并不是完美的，目前对于它的关键步骤EMD分解的研究尚不完善，缺乏一些理论基础。从上图我们也可以看到，HHT在低频区域可能会出现一些不存在的频率分量。

参考

https://blog.csdn.net/Forlogen/article/details/88535027

https://www.sohu.com/a/246972969_607269

唐晓初. 小波分析及其应用[M]. 重庆大学出版社, 2006.

Marple L . Computing the discrete-time “analytic” signal via FFT[J]. IEEE Transactions on Signal Processing, 1999.

罗利春. 用希尔伯特变换构造解析信号进行时频分析[J]. 航天电子对抗, 2003(03):27-30.

Qian S , Chen D . Joint time-frequency analysis[J]. IEEE Signal Processing Magazine, 1999, 16(2): P.52-67.