8.3 短时傅里叶变换(STFT)
（1）基本概念
傅里叶变换只反映出信号在频域的特性，无法在时域内对信号进行分析。为了将时域和频域相联系，Gabor于1946年提出了短时傅里叶变换(short-time Fourier transform
，STFT)，其实质是加窗的傅里叶变换。STFT的过程是：在信号做傅里叶变换之前乘一个时间有限的窗函数 h(t)，并假定非平稳信号在分析窗的短时间隔内是平稳的，通过窗函数
h(t)在时间轴上的移动，对信号进行逐段分析得到信号的一组局部“频谱”。
傅里叶变换非常擅长分析那些频率特征均一稳定的平稳信号。但是对于非平稳信号，傅立叶变换只能告诉我们信号当中有哪些频率成分——而这对我们来讲显然是不够的。我们还想知
道各个成分出现的时间。知道信号频率随时间变化的情况，各个时刻的瞬时频率及其幅值——这也就是时频分析。
（2）短时傅里叶变换原理介绍
短时傅里叶变换的核心思想是对信号进行分段处理从而达到对信号的时频分析，如何对信号进行分段处理呢?一般我们通过加窗的方法来截取信号的片段。定义一个窗函数 \(\textrm{w}(t)\) ，比如这样。
将窗函数位移到某一中心点 \(\tau\)，再将窗函数和原始信号相乘就可以得到截取后的信号 y(t)
\[\mathrm{y}(\mathrm{t})=\mathrm{x}(\mathrm{t}) \cdot \mathrm{w}(\mathrm{t}-\tau)\]
一般我们很少选用矩形窗，因为矩形窗简单粗暴的截断方法会产生的频谱泄露以及吉布斯现象，不利于频谱分析。
对原始信号 x(t) 做 STFT 的步骤如下。首先将窗口移动到信号的开端位置，此时窗函数的中心位置在 \(t = \tau_0\) 处，对信号加窗处理
\[\mathrm{y}(\mathrm{t})=\mathrm{x}(\mathrm{t}) \cdot \mathrm{w}(\mathrm{t}-\tau_0)\]
然后进行傅里叶变换：
\[\mathrm{X}(\omega)=\mathcal{F}(\mathrm{y}(\mathrm{t}))=\int_{\infty}^{+\infty} \mathrm{x}(\mathrm{t}) \cdot \mathrm{w}\left(\mathrm{t}-\tau_{0}\right) \mathrm{e}^{-\mathrm{j} \omega \mathrm{t}} \mathrm{dt}\]
由此得到第一个分段序列的频谱分布 \(X(\omega)\) 。在现实应用中，由于信号是离散的点序列，所以我们得到的是频谱序列 X[N]。
为了便于表示，我们在这里定义函数 \(S(\omega, \tau)\) ，它表示，在窗函数中心为 \(\tauτ\) 时，对原函数进行变换后的频谱结果 \(X(\omega)\) ，即：
\[\mathrm{S}(\omega, \tau)=\mathcal{F}(\mathrm{x}(\mathrm{t}) \cdot \mathrm{w}(\mathrm{t}-\tau))=\int_{\infty}^{+\infty} \mathrm{x}(\mathrm{t}) \cdot \mathrm{w}(\mathrm{t}-\tau) \mathrm{e}^{-\mathrm{j} \omega \mathrm{t}} \mathrm{dt}\]
对应到离散场景中，\(S[\omega, \tau]\) 就是一个二维矩阵，每一列代表了在不同位置对信号加窗，对得到的分段进行傅里叶变换后的结果序列。
完成了对第一个分段的FFT操作后，移动窗函数到 \(\tau_1\) 。把窗体移动的距离称为 Hop Size。移动距离一般小于窗口的宽度，从而保证前后两个窗口之间存在一定重叠部分，我们管这个重叠叫 Overlap
重复以上操作，不断滑动窗口、FFT，最终得到从 \(\tau_0 \sim \tau_N\)
最终我们得到的 S ，就是 STFT 变换后的结果。
（3）matlab中的应用
在进行短时傅里叶变换matlab当中的应用分析前，我们再次明确一下目标，即为什么要进行时频分析？先说一下傅里叶变换的缺点：
1.傅里叶变换只能获取一段信号总体上包含哪些频率的成分，但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。
2.对于信号中的突变，傅里叶变换很难及时捕捉。而在有些场合，这样的突变往往是十分重要的。
即傅里叶变换非常擅长分析那些频率特征均一稳定的平稳信号。但是对于非平稳信号，傅立叶变换只能告诉我们信号当中有哪些频率成分——而这对我们来讲
显然是不够的。我们还想知道各个成分出现的时间。知道信号频率随时间变化的情况，各个时刻的瞬时频率及其幅值——这也就是时频分析
代码示例
close all; clear; clc;
fs = 1000;
t = 0:1/fs:1 - 1/fs;
% 窗口大小，推荐取 2 的幂次
wlen = 256;
% hop size 即移动步长，一般要取一个小于 wlen 的数，推荐取 2 的幂次
hop = wlen/4;
% FFT 点数，理论上应该不小于wlen，推荐取 2 的幂次
nfft = 256;
%示例信号
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)];
figure;
subplot(2, 2, 1);
plot(x);
% 随便选的一个窗函数
win = blackman(wlen, 'periodic');
[S, f, t] = mystft(x, win, hop, nfft, fs);
subplot(2, 2, 2);
PlotSTFT(t,f,S);
%对上述生成的x信号进行倒叙排列
x = x(end:-1:1);
subplot(2, 2, 3);
plot(x);
win = blackman(wlen, 'periodic');
[S, f, t] = mystft(x, win, hop, nfft, fs);
subplot(2, 2, 4);
PlotSTFT(t,f,S);
可以看到，在 FFT 中无法区分的频谱图像在 STFT 中区分就非常明显，可以看出按照不同的时间分段，频谱分布的变化
为了更好地理解，将右上角的图做一次三维旋转：
可以非常清晰地看出频率分布随时间的变换。注意到分界线处存在异常的高频成分（就是 STFT 图像中那三条竖线），这是因为时域信号突变导致的高频成分。
spectrogram函数与stft函数的对比
在老的版本中，Matlab 中 STFT 的函数名为 spectrogram，而在新版本中，引入了新的函数 stft，用法和我上面的实现的程序基本一致
close all; clear; clc;
fs = 1000;
t = 0:1/fs:1 - 1/fs;
% 窗口大小，推荐取 2 的幂次
wlen = 256;
% hop size 即移动步长，一般要取一个小于 wlen 的数，推荐取 2 的幂次
hop = wlen/4;
% FFT 点数，理论上应该不小于wlen，推荐取 2 的幂次
nfft = 256;
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)];
figure;
subplot(1, 3, 1);
win = blackman(wlen, 'periodic');
[S, f, t] = mystft(x, win, hop, nfft, fs);
PlotSTFT(t,f,S);
title('My STFT');
subplot(1, 3, 2);
[S1, f1, t1] = spectrogram(x, win, wlen - hop, nfft, fs);
PlotSTFT(t1, f1, S1);
title('spectrogram');
subplot(1, 3, 3);
[S2, f2, t2] = stft(x, fs, 'Window', win, 'OverlapLength',wlen - hop,'FFTLength',nfft);
PlotSTFT(t2, f2, S2);
title('stft');
生成图像如下：
要注意的是spectrogram 输出的是单边谱，而 stft 输出的是双边谱，spectrogram 还可以输出功率谱，而 stft 不行。
STFT变换的缺点
简单来说stft变换就是加窗的傅里叶变换，它克服了福利叶变换局部分析的能力不足的缺点。但其窗口一经确定，便不能改变，所以其不能根据待分析信号频率的改变而改变分辨率。
窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差；
窗太宽，时域上又不够精细，时间分辨率低。
从定量的角度来看，STFT的时间分辨率取决于滑移宽度 H ，而频率分辨率则取决于 \(\frac{F_s}{H}\) ，显然，一方的增加必然意味着另一方的减小。这就是所谓的时频测不准原理，具体关系为：
\[\Delta \mathrm{t} \cdot \Delta \mathrm{f} \geqslant \frac{1}{4 \pi}\]
另外，固定的窗口大小过于死板。对低频信号而言，有可能连一个周期都不能覆盖；对高频信号而言，可能覆盖过多周期，不能反映信号变化。
参考：
时频分析之STFT.