转载时频分析之STFT:短时傅里叶变换的原理与代码实现(非调用Matlab API)。
这篇文章介绍了FFT的缺点、STFT的原理、STFT在matlab上运用以及STFT的缺点,写的非常非常好。下面这篇文章主要对原理、STFT后生成的图进行个人注释。
原理
补充一张原理图(图片来源https://www.jianshu.com/p/7e160442830f)。
图像分析
close all; clear; clc;
fs = 1000;
t = 0:1/fs:4-1/fs;
L=4000;
% 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)];
% fplot(@(t) 10 * cos(2 * pi * 10 * t),[0 1],'g')
% hold on
% fplot(@(t) 20 * cos(2 * pi * 20 * t),[1 2],'g')
% hold on
% fplot(@(t) 30 * cos(2 * pi * 30 * t),[2 3],'g')
% hold on
% fplot(@(t) 40 * cos(2 * pi * 40 * t),[3 4],'g')
% hold off
% grid on
signal=fenduan(t);
plot(t,signal);
f = (0:(L-1))*(fs/L);
Y = fft(signal);
figure
subplot(2,1,1);
plot(f,abs(Y*2/L))
xlim([0,100])
% 窗口大小,推荐取 2 的幂次
wlen = 256;
% hop size 即移动步长,一般要取一个小于 wlen 的数,推荐取 2 的幂次
hop = wlen/4;
% FFT 点数,理论上应该不小于wlen,推荐取 2 的幂次
nfft = 256;
win = blackman(wlen, 'periodic');
[S, F, T] = spectrogram(signal, win, wlen - hop, nfft, fs);
subplot(2,1,2);
PlotSTFT(T, F, S,win);
title('spectrogram');
其中
function m=fenduan(t)
m=10*cos(2*pi*20*t).*(t>=0 & t<1)+80*cos(2*pi*80*(t-1)).*(t>1 & t<=2)+60*cos(2*pi*60*(t-2)).*(t>2 & t<=3)+40*cos(2*pi*40*(t-3)).*(t>3 & t<=4);
end
signal的时域图
它的频域图和经STFT处理后的图形
这张图就可以表面STFT的两个优点。
1.傅里叶变换只能获取一段信号总体上包含哪些频率的成分,而STFT处理可以对各成分出现的时刻加以展现。
2.对于信号中的突变,傅里叶变换很难及时捕捉。而STFT可以捕捉信号的突变。
最后放两张STFT处理后信号图像的三维旋转图。