Skip to content

短时傅里叶变换

先来看一下stft短时傅里叶变换

clc;clear;close all;
fs = 10e6;
n=10000;
f1=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];
t1 = (0:4*n-1)'/fs;
%time domain
plot(t1,sig)

figure_0.png

s = spectrogram(sig,256);
t = linspace(0, 4*n/fs, size(s,1));
f = linspace(0, fs/2, size(s,2));
imagesc(t,f,20*log10(abs(s)));
xlabel("samples");
%sample是时间，每个采样点对应的频率大小。二维函数。
ylabel("freqency");
colorbar;

figure_1.png

%可以修改一下得到更好的图片
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);
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('Freqency');
colorbar;

figure_2.png

%功率谱
[s, f, t, p] = spectrogram(sig, window, 256, f, fs);
figure;
imagesc(t, f, p);xlabel('Samples'); ylabel('Freqency');
colorbar;

figure_3.png

example10.9

alpha_0 = 15*pi*1e-6;
n=0:20000;
x = cos(alpha_0*n.^2);
plot(n,x);
xlim([0,20000]);
title("time domain");

figure_4.png

%整个频率域
X = fftshift(fft(x, 20000));
plot((abs(X)));
title("frenqency domain");

figure_5.png

%表示这401个点包含的频率（基本不变，不涉及其他部分）
window1 = hamming(401);
x_piece1 = x(5000:5000+400).*window1';
X_piece1 = fftshift(fft(x_piece1, 4096));
omega = -1:2*1/4096:1-2*1/4096;
plot(omega, abs(X_piece1));

figure_6.png

x_piece2 = x(15000:15000+400).*window1';
X_piece2 = fftshift(fft(x_piece2, 4096));
omega = -1:2*1/4096:1-2*1/4096;
plot(omega, abs(X_piece2));

figure_7.png

example10.10 谱图

下面手写一个stft做成的谱图

n1=0:20000;
n2=20001:25000;
n3=25001:30000;
alpha_0 = 15*pi*1e-6;
y1 = cos(alpha_0*n1.^2);
y2 = cos(0.2*pi*n2);
y3 = cos(0.2*pi*n3)+cos(0.23*pi*n3);
n=[n1,n2,n3];
y=[y1,y2,y3];

%401points window stft
L=401;
h_win = hamming(L);
nfft=1024;
res = zeros(nfft/2,length(n));
for i=0:n(end)-L
    y_piece = y(1+i:i+L).*h_win';
    Y_piece = fft(y_piece,nfft);
    Y_piece_half = Y_piece(1:nfft/2);
    res(:,i+1) = Y_piece_half';
end
f =  0:1/512:1-1/512;
imagesc(n,f,abs(res));

figure_8.png

%101points window stft
L=101;
h_win = hamming(L);
nfft=1024;
res = zeros(nfft/2,length(n));
for i=0:n(end)-L
    y_piece = y(1+i:i+L).*h_win';
    Y_piece = fft(y_piece,nfft);
    Y_piece_half = Y_piece(1:nfft/2);
    res(:,i+1) = Y_piece_half';
end
f =  0:1/512:1-1/512;
imagesc(n,f,abs(res));

figure_9.png

%matlab自带的函数
% s = spectrogram(y, window1);
% s = spectrogram(x,window,noverlap,nfft)
%noverlap是重叠的样本数，默认是floor(ns/4.5)
%nfft是fft点数。
% imagesc(abs(s))

看上面窗长100个点和400个点的结果，怎么分析？

每一个竖直方向上的切片，都对应着一次dtft。我这里1024个点fft只取0到pi的结果（反正对称）

竖直方向上亮线的粗细是有什么决定的呢？我们已经学过了就是由窗的主瓣长度决定的。而窗的主瓣长度，又是直接由窗长决定的：对于hamming窗，有 \(\Delta_{\textrm{ml}} =8\frac{\pi }{M}\) ，M就是窗长。

窗长为401个点的时候，主瓣很短，亮线很细，但是在20000个点的突变这里，由于窗比较长，所以突变会模糊一点。

窗长为101个点的时候，就反过来，对应更长的主瓣，更粗的亮线，但是对信号特性的变化描述更清晰。

窗长和频率分辨率成正比，和时间分辨率成反比。