第二篇：傅里叶变换与短时傅里叶变换

在上一篇文章中，我们简要介绍了傅里叶变换的一些背景、作用及存在的问题，并通过一些例子直观地感受了一下它的效用。其中一个例子表明，尽管两个原始信号并不相同，但经过傅里叶变换后，它们在频域上的分布却非常相似。

为什么会出现这样的情况呢？

本篇文章就试图回答这个问题，以及上一篇文章遗留的两个问题：

1. 信号的观察窗口多小叫足够小？
2. 观察窗口内的信号如果不是平稳的怎么办？

1、傅里叶变换：信号处理的里程碑

在19世纪初叶，法国数学家J. Fourier发现，任何周期函数都可以表示为复指数函数的无穷和。后来又经过许多数学家许多年的工作，慢慢形成了现在所为人熟知的傅里叶变换。

一般情况下，若傅里叶变换不加任何限定语，则指的是连续傅里叶变换（连续函数的傅里叶变换）。

定义傅里叶变换有许多种不同的形式，这里采用如下的定义：傅里叶变换将可积函数$x$表示成复指数函数的积分形式。公式表示如下：

$$X(f)={\int }_{-\infty }^{\infty }x(t)\cdot e^{-2\pi ift}dt$$

$f$表示频率（单位是Hz），为任意实数；$t$表示时间（单位是s）。

在适当的条件下，傅里叶变换是可逆的：

$$x(t)={\int }_{-\infty }^{\infty }X(f)\cdot e^{2\pi ift}df$$

注意，$x$和$X$分别代表信号在时域和频域中的表示。

第一个公式可以被理解为如下的含义：时域信号$x(t)$乘上一个由指定频率组成的指数项（$e^{-2\pi ift}$），然后在整个时间轴上进行积分，得到了频域信号。

对第一个公式中的指数项应用欧拉公式，得到如下表示：

$$e^{-2\pi ift}=cos(2\pi ft)-i\cdot sin(2\pi ft)$$

所以，傅里叶变换的实质是将原始的时域信号乘上一个由固定频率（$f$）构成的余弦与正弦函数后，再在整个时间轴上进行积分。

如果积分的结果是一个比较大的值，则说明频率$f$是信号$x(t)$的一个主要组成频率；如果积分的结果是一个比较小的值，那么就说频率$f$不是信号$x(t)$的主要组成频率；如果积分值是0，则说明信号$x(t)$中不存在频率$f$。

注意到，积分是在整个时间轴上进行的。这也就是说，对于任一频率$f$，无论它是只出现在某个特定的时间段，还是存在于信号的整个生命周期中，其对最终结果的影响是一样的。

这就回答了本文提出的第一个问题（不同的信号的频域表现是相似的）。也就是说，傅里叶变换仅能够识别出信号中的频谱分量有哪些，但无法进一步判断其出现的时间信息。

数学家们当然意识到了这个问题，并且也提出了一些解决方案，短时傅里叶变换就是其中的一种。

2、短时傅里叶变换（The Short Time Fourier Transform, STFT）

上一篇文章中已经提到，「短时」其实就是选择一个较短的时间窗口对信号进行观察，这个窗口是如此的短以至于我们可以认为观察到的信号是平稳的。

信号平稳了，就可以对其应用傅里叶变换了。

图2.1

上图的信号是通过程序生成的，这个信号的频率每250ms变化一次，在各个区间内都是平稳的。如果选择的观察窗口不大于250ms，那么将窗口沿时间轴从0开始向右移动时，在某些情况下窗口内观察到的信号一定会出现平稳的情况。

有读者可能要问了，无论选择什么长度的时间窗口，将它向右移动时，总会跨过频率的突变点，此时覆盖的信号就是不平稳的，这该怎么办呢？答案是没有任何办法，只能对这部分信号也应用傅里叶变换。这里需要注意的是，窗口信号的平稳性是应用短时傅里叶变换的前提假设，如果不满足这个假设，那么得到的变换结果可能没有参考价值。

这其实也就回答了第一篇遗留的两个问题：如果窗口信号不平稳，那也只能照样对其进行傅里叶变换；要求窗口“尽可能的小”其实没有一个统一的标准，这要视信号的频率而定，一般在应用中，往往会尝试多次不同的窗口来确定合适的取值。

相较于基础的FT，短时傅里叶变换的公式中会多一个表示信号观察时间窗口函数$w(\cdot )$，$w(\cdot )$的取值必须（尽可能）保证在其区间内的信号是平稳的。

由于「时间窗口」是其显著的特征，因此有人把STFT称为「窗口傅里叶变换」。

3、从例子入手

用到了python(3)的第三方库包括：numpy、matplotlib、scipy，需要自行安装。

3.1 对典型的非平稳信号进行STFT

现在，我们给出图2.1所示信号的生成代码：

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS'] 

x = np.linspace(0, 1000, 1000)
y1 = np.sin(2*np.pi*50*x[:250])
y2 = np.sin(2*np.pi*30*x[250:250*2])
y3 = np.sin(2*np.pi*10*x[250*2:250*3])
y4 = np.sin(2*np.pi*5*x[250*3:])
y = np.concatenate((y1, y2, y3, y4))

plt.plot(x, y)
plt.xlabel('时间[ms]')
plt.ylabel('振幅')
plt.show()

由此易知，它的频率每250ms变化一次，从0时刻开始，其频率分别为50Hz、30Hz、10Hz和5Hz。

接下来，我们对这个信号应用STFT并对结果进行可视化：

import scipy
# 应用STFT
f, t, Zxx = scipy.signal.stft(y, fs=1000, nperseg=250)
# 对结果进行可视化
plt.pcolormesh(t*1000, f, np.abs(Zxx), vmin=0, shading='gouraud')
plt.ylim(top=100)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [ms]')
plt.show()

效果如下：

图2.2

尽管这还是一个二维图像，但是通过颜色的深浅读者仍可以了解到第三个纬度的信息。从图中可以看出，信号的频率非常明显地分为4个部分，且各个频率分量在时间上的持续情况也被区分了出来。

读者需要注意，上图中的坐标轴上的数值并非是时间和频率的真实数值，而是经过了一定比例的缩放（通过调整scipy.signal.stft中的fs参数的值，你可以看到纵轴的变化）。

3.2 对典型的平稳信号进行STFT

在本小节，我们构造一个平稳的信号，并对它应用STFT以观察效果。代码如下：

x = np.linspace(0, 1000, 1000)
y1 = np.sin(2*np.pi*50*x)
y2 = np.sin(2*np.pi*30*x)
y3 = np.sin(2*np.pi*10*x)
y4 = np.sin(2*np.pi*5*x)
y = y1 + y2 + y3 + y4
# 原始信号图
plt.plot(x, y)
plt.show()
# 频域图
f, t, Zxx = scipy.signal.stft(y, fs=1000, nperseg=250)
plt.pcolormesh(t*1000, f, np.abs(Zxx), vmin=0, shading='gouraud')
plt.ylim(top=100)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [ms]')
plt.show()

结果如下：

图2.3

图2.4

从频域图上可以看出，这个信号在整个时间轴上都有4个频率值。

3.3 窗口对结果的影响

接下来，我们将使用图2.1中的信号来观察一下不同的窗口宽度对变换结果的影响。

本节选用了高斯窗口函数，在调用时需要额外传入其标准差以控制窗口宽度：

f, t, Zxx = scipy.signal.stft(y, window=('gaussian', 80), fs=1000)

我们知道，标准差越大，整个高斯函数的图像越扁平，对应到STFT中，则体现为在时间轴上的跨度越大。按照上面的猜想，时间跨度越大，则频率分辨率越好、时间分辨率越差。为此，我们设置候选的标准差为[10,40,80]，分别查看一下其转换结果：

f, t, Zxx = scipy.signal.stft(y, window=('gaussian', 10), fs=1000)  # 将此处的10分别换成40、80即可调整标准差
plt.pcolormesh(t*1000, f, np.abs(Zxx), vmin=0, shading='gouraud')
plt.ylim(top=100)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [ms]')

图2.5，标准差为10

图2.6，标准差为40

图2.7，标准差为80

从上面三个图可以看出，随着标准差的增大（窗口逐渐变宽），各个频率分量在频率轴上逐渐「聚焦」，即分辨频率的能力增加，不确定性降低。

读者可以自行尝试继续增加窗口宽度，会发现频谱中各条频率“横线”逐渐横向拉长，彼此交叉。这说明了随着窗口的变宽，覆盖的时间范围增加，导致其在时间上的区分度降低，即时间分辨率降低。

上述例子直观地展示了在使用STFT时所面临的选择困境，即：到底如何选择窗口函数？较窄的窗口函数具有较好的时间分辨率但较差的频率分辨率；较宽的窗口函数则恰恰相反。此外，较宽的窗口函数还会面临截取的信号不平稳的问题。

因此，如何选择窗口函数就成了一个具体问题具体分析、见仁见智的问题了。至此，终于可以引出我们的小波变换了——它在一定程度上就是来解决这个选择困境的。