Matlab中短时傅⾥叶变换spectrogram和stft的⽤法
在Matlab中，做短时傅⾥叶变换需要使⽤函数spectrogram，⽽在Matlab2019中，引⼊了⼀个新的函数stft，下⾯我们就来看下这两个函数都如何使⽤。
短时傅⾥叶变换的基本原理就是将数据分段加窗，做fft，在分段时会有overlap，因此⼀个向量的短时傅⾥叶变换结果是⼀个矩阵。了解了这点，下⾯的函数及参数就更加容易理解了。
尚德实验学校
spectrogram
参数列表
先来看spectrogram函数，在更早期的版本中，这个函数的名字是specgram，⼏种常⽤的⽤法如下：
humanaspectrogram(x)
s = spectrogram(x)
s = spectrogram(x, window)
s = spectrogram(x, window, noverlap)
争做文明学生s = spectrogram(x, window, noverlap, nfft)
s = spectrogram(x, window, noverlap, nfft, fs)
[s, f, t] = spectrogram(x, window, noverlap, nfft, fs)
[s, f, t] = spectrogram(x, window, noverlap, f, fs)
[s, f, t, p] = spectrogram(x, window, noverlap, f, fs)
其中，
x表⽰输⼊信号；
window表⽰窗函数，如果window的值是⼀个整数，那么被分段的x的每⼀段的长度都等于window，并采⽤默认的Hamming窗；如果window是⼀个向量，那么被分段后每⼀段的长度都等于length(window)，且输⼊的向量即为所要加的窗函数；
overlap表⽰两段之间的重合点数，overlap的值必须要⼩于窗长，如果没有指定overlap，默认是窗长的⼀半，即50%的overlap；
nfft表⽰fft的点数，fft的点数跟窗长可以是不同的，当没有指定该参数时，Matlab会取max(256,
2^(ceil(log2(length(window)))))，即当窗长⼩于256时，fft的点数是256；当窗长⼤于256时，fft的点数取⼤于窗长的最⼩的2的整数次幂；
fs表⽰采样率，⽤来归⼀化显⽰使⽤；
f表⽰显⽰的频谱范围，f是⼀个向量，长度跟s的⾏数相同；
当x是实信号且nfft为偶数时，s的⾏数为(nfft/2+1)
当x是实信号且nfft为奇数时，s的⾏数为(nfft+1)/2
当x是复信号时，s的⾏数为nfft
当在输⼊的参数列表中指定f后，函数会在f指定的频率处计算频谱图，返回的f跟输⼊的f是相同的；
t表⽰显⽰的时间范围，是⼀个向量，长度跟s的列数相同；
p表⽰功率谱密度，对于实信号，p是各段PSD的单边周期估计；对于复信号，当指定F频率向量时，P为双边PSD;
Examples
⾸先，⽣成信号如下，4个点频信号拼接起来：
clc;clear all;clo all;
fs = 10e6;
n = 10000;
independentlyf1 = 10e3; f2 = 50e3; f3 = 80e3; f4 = 100e3;
t = (0:n-1)'/fs;
sig1 = cos(2*pi*f1*t);
sig2 = cos(2*pi*f2*t);
sig3 = cos(2*pi*f3*t);
sig4 = cos(2*pi*f4*t);
sig = [sig1; sig2; sig3; sig4];
信号的时域波形如下：
直接调⽤spectrogram(sig)，可得如下结果，图中默认横轴是频率，纵轴是时间truelove
为了绘图更灵活，我们不直接⽤spectrogram绘图，⽽且求出s后，再对s单独绘图，这次我们指定window的⼤⼩为256
s = spectrogram(sig, 256);
t = linspace(0, 4*n/fs, size(s,1));
f = linspace(0, fs/2, size(s,2));
figure;
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
freeze是什么意思
colorbar;
body language
noverlap默认是50%，现在我们把它设为window的长度减1，即每次的步进为1
s = spectrogram(sig, 256, 255);
t = linspace(0, 4*n/fs, size(s,1));
f = linspace(0, fs/2, size(s,2));
figure;
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
colorbar;暴发户 英语
再加上nfft和fs参数，我们指定fft点数就是窗长
s = spectrogram(sig, 256, 128, 256, fs);
这个的图形跟之前⼀样，不再画了
如果在返回值中，增加f和t，这样我们下⾯就不⽤再重新定义f和t了
[s, f, t] = spectrogram(sig, 256, 128, 256, fs);
figure;
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
colorbar;
从上⾯的图中我们可以看出，我们的4个信号的频率都⽐较⼩，⽽画出来的图显⽰的频谱范围⽐较⼤，导致下⾯很⼤⼀部分信息我们其实都不需要。这时，我们就可以通过指定f的区间来计算频谱。为了显⽰效果更好，我们把其他参数也调⼀下
window = 2048;
noverlap = window/2;
f_len = window/2 + 1;
f = linspace(0, 150e3, f_len);
[s, f, t] = spectrogram(sig, window, noverlap , f, fs);
figure;
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
colorbar;
最后再把功率谱密度的返回值加上
[s, f, t, p] = spectrogram(sig, window, nfft, f, fs);
figure;
imagesc(t, f, p);xlabel('Samples'); ylabel('Freqency');
colorbar;
stft
这个函数在Matlab的解释并不是很多，example也只写了两个，但⽤法⽐较简单：
window = 2048;
noverlap = window/2;
nfft = window;
[s, f, t, p] = spectrogram(sig, window, noverlap, nfft, fs);
figure;
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
title('使⽤spectrogram画出的短时傅⾥叶变换图形');
colorbar;
ss = stft(sig,fs,'Window',hamming(window),'OverlapLength',window/2,'FFTLength',nfft);
figure;
imagesc(t, f, 20*log10((abs(ss(1024:end,:)))));xlabel('Samples'); ylabel('Freqency');
title('使⽤stft画出的短时傅⾥叶变换图形');
colorbar;
bottom的反义词
欢迎关注微信公众号：Qunat_Times
inbound欢迎⼤家学习我的课程：
