时频分析之STFT：短时傅里叶变换的原理与代码实现（非调用Matlab API）

1. 引言

在信号分析中，傅里叶变换可称得上是神器。但在实际应用中，人们发现它还是存在一些不可忽视的缺陷。

为了便于叙述考察以下两种情形：

Case 1

考察这样一个函数：

fs = 1000;
t = 0:1/fs:1 - 1/fs;
x = [10 * cos(2 * pi * 10 * t), 20 * cos(2 * pi * 20 * t),...
        30 * cos(2 * pi * 30 * t), 40 * cos(2 * pi * 40 * t)];

绘制这个函数的时域图像和经过傅立叶变换后的频谱图像，长这个样子：

现在把信号反转过来：

x = [10 * cos(2 * pi * 10 * t), 20 * cos(2 * pi * 20 * t),...
        30 * cos(2 * pi * 30 * t), 40 * cos(2 * pi * 40 * t)];
x = x(end:-1:1);

再次绘制时域和频域的图像，它长这样：

不难发现，尽管这两个信号的时域分布完全相反，但是它们的频谱图是完全一致的。显然，FFT无法捕捉到信号在时域分布上的不同。

Case 2

考察一个普普通通的信号:

fs = 1000;
t = 0:1/fs:1-1/fs;

x = 2 * cos(2 * 10 * t) + 4 * sin(2 * 30 * t);

同样绘制它的时域以及频域图像：

现在给信号加入一个高频突变：

sharp = zeros(1, length(x));
% 给信号中间加一个突变
sharp(501:510) = 5 * cos(2 * pi * 100 * linspace(0, 1, 10));
x = x + sharp;

然后绘图：

对比两个信号的时域图，我们能很明显发现在第二个信号中央的部分出现了一个突变扰动。然而在频域图中，这样的变化并没有很好的被捕捉到。注意到红框中部分，显然傅里叶变换把突变解释为了一系列低成分的高频信号的叠加，并没有很好的反应突变扰动给信号带来的变化。

为什么我们需要时频分析

通过以上的两个例子，我们不难发现傅立叶变换的缺陷。

第一个例子告诉我们，傅里叶变换只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。

第二个例子告诉我们，对于信号中的突变，傅里叶变换很难及时捕捉。而在有些场合，这样的突变往往是十分重要的。

当然如果非要硬杠，也不是完全没办法——这就需要需分析相位谱了，但在实际应用中，有谁会不嫌麻烦地去看相位谱呢？

总而言之，傅里叶变换非常擅长分析那些频率特征均一稳定的平稳信号。但是对于非平稳信号，傅立叶变换只能告诉我们信号当中有哪些频率成分——而这对我们来讲显然是不够的。我们还想知道各个成分出现的时间。知道信号频率随时间变化的情况，各个时刻的瞬时频率及其幅值——这也就是时频分析（引用自知乎）。

所谓时频分析，就是既要考虑到频率特征，又要考虑到时间序列变化。常用的有两种方法：短时傅里叶变化，以及小波变换。本文我们只介绍短时傅里叶变换

2. 短时傅里叶变换原理

短时傅里叶变换的思路非常直观：既然对整个序列做FFT会丢失时间信息，那我一段一段地做FFT不就行了嘛！这也正是短时傅里叶变换名称的来源，Short Time Fourier Transorm，这里的 Short Time 就是指对一小段序列做 FFT。

那么怎么一段一段处理呢？直接截取信号的一段来做 FFT 吗？一般我们通过加窗的方法来截取信号的片段。定义一个窗函数 $\text{w}(t)$，比如这样。

将窗函数位移到某一中心点 $\tau$，再将窗函数和原始信号相乘就可以得到截取后的信号 y(t)。

$$y(t)=x(t)\cdot \text{w}(t-\tau )$$

前面提到的直接截取的方法其实就是对信号加一个矩形窗，不过一般我们很少选用矩形窗，因为矩形窗简单粗暴的截断方法会产生的频谱泄露以及吉布斯现象，不利于频谱分析。更多关于窗函数的内容，可以看这里：加窗法。

对原始信号 $x(t)$ 做 STFT 的步骤如下。

首先将将窗口移动到信号的开端位置，此时窗函数的中心位置在 $t={\tau }_0$处，对信号加窗处理

$$y(t)=x$$