短时傅里叶变换原理及其MATLAB实现(Short Time Fourier Transform,STFT)
1.短时Fourier变换原理(STFT原理)
信号x(t)短时Fourier变换定义为:
其中w(τ)为窗函数。
X(ω,t)中的时间t表示窗函数w(τ−t)的位置,随着窗函数在整个区间上的滑动,可获得信号x(τ)在 t 附近区域上对应的频谱。
信号短时Fourier变换是一种常用的信号时频分析方法。
2.DFT中的STFT原理
信号x(t)的STFT是一个积分运算,在实际计算中也可通过DFT来实现,即
其中:窗函数w[k]的宽度为N,x[k]为连续信号x(t)的抽样。
若抽样频率为fs,则存在 t=kT, T=1/fs。
3.时间分辨率和频谱分辨率
(2)时间分辨率
时间分辨率由时窗宽度Tp决定,公式为
Tp越小,时间分辨率越高。
(2)频谱分辨率
频谱分辨率是指分辩信号中相邻谱峰的能力,公式为
∆fc越小,频谱分辨率越高。
4.不足
信号的STFT虽然能在一定程度上改善Fourier变换的不足,实现信号的时频分析,但其时间分辨率固定不变,因而不能有效地反映信号的突变程度,其应用受到局限。
小波分析拓展了信号STFT,实现了一种新的时频分析方法,其时窗可以随着信号的频率增高而缩小,频率降低而增大,有效地解决信号短时Fourier变换的缺陷,因而得到广泛应用。
5.Matlab代码实现
clear;clc;Fs=1e4;t=0:1/Fs:4;t1=0:1/Fs:1;%10Hz的信号x1=[sin(2*pi*10*t1)';zeros(Fs*3,1)]; %1000Hz的信号x2=[zeros(Fs,1);sin(2*pi*1000*t1)';zeros(Fs*2,1)];%2000Hz的信号x3=[zeros(Fs*2,1);sin(2*pi*2000*t1)';zeros(Fs,1)];%3000Hz的信号x4=[zeros(Fs*3,1);sin(2*pi*3000*t1)'];figure(1);subplot(411);plot(t,x1);subplot(412);plot(t,x2);subplot(413);plot(t,x3);subplot(414);plot(t,x4);x=x1+x2+x3+x4;figure(2);subplot(211);plot(t,x);title('时域图'); axis([0 4 -2 2]);%FFT代码X=fft(x);N=length(x);f=(1:N/2)*Fs/N;subplot(212);plot(f,abs(X(1:N/2)));xlabel('f (Hz)');ylabel('|X(f)|');title('频谱图'); %STFT关键代码figure(3);spectrogram(x,128,120,128,Fs,'yaxis');title('时频图');