强大的MATLAB2019——短时傅里叶变换stft函数

MATLAB最新版本2019b，推出了一个超好用的函数——stft。它能非常方便地为你实现加窗傅里叶变换。

什么是加窗傅里叶变换呢，简单来说就是把一段时间信号切成一段一段的，再给每一段分别进行傅里叶变换，可以想见，只要你切得足够细，就能得到频率随时间的变化规律。

加窗傅里叶变换的动画演示可以查看本人视频

当然这里有些讲究。首先切信号不能直接像切面条似的把信号截出来，而是应该用两端为0的时窗去截，否则会产生吉布斯现象，即在截断附近出现大的信号跳变，造成原始信号畸变。

其次呢，我们不能像切面条一样把信号切成一段一段不重叠的段，而是应该互相有所重叠，否则的话，万一你切段的位置正好切在信号上怎么办呢。

明确这两点，我们就来看一下MATLAB的stft函数吧！

例子1：

产生两秒钟的压控振荡器输出，该输出由以10 kHz采样的正弦波控制。

fs = 10e3; 

t = 0:1/fs:2; 

x = vco(sin(2*pi*t),[0.1 0.4]*fs,fs);

计算并绘制信号的STFT。使用长度为256且形状参数β = 5的Kaiser窗口。指定重叠的长度为220个样本，DFT的长度为512点。用默认的颜色图和视图绘制STFT。

stft(x,fs,'Window',kaiser(256,5),'OverlapLength',220,'FFTLength',512);

看一下各个参数，

Window后面可以选择时窗，这个时窗的定义呢就和滤波器里时窗的定义方式一样~

OverlapLength就是重叠宽度，也就是相邻两个的重叠宽度。

FFTLength就是每个小段的长度。

也可选用另一种Hamming时窗。

stft(x,fs,'Window',hamming(128,'periodic'),'OverlapLength',50);

如果我们想看震撼的立体效果怎么办呢，我们只要将视角调成俯视就好啦！

view(-45,65) colormap jet

接下来用我自己生成的do re mi fa so la xi(40号音 到 51号音)七个音去算短时距傅里叶，得到很清晰的结果~

注意到我的音色函数是写了泛音列（基音的整数倍）呢。泛音列的原理可以看我的视频呦 

那么如何根据stft函数得到每一段的频率数值呢，也很简单。

令s=stft(参数)，得到一个时窗长度×n的矩阵，每一列就是每个时窗里的fft变换结果。

和由FFT结果画幅值频谱图一样，取出某一列画图就行了。

abs(s[:,m]);

就能生成这个时间段的幅值频谱图啦！

参考资料：

https://ww2.mathworks.cn/help/signal/ref/stft.html