短时傅立叶变换的代码程序

短时傅立叶变换试验为了克服傅立叶变换的时频局部化方面的不足，也是为了对时域信号作局部分析，D.Gabor于1946年提出了窗口傅立叶变换（简记为WFT）。WFT的公式形式()(,)()()jwtRGfwbftwtbedt()(,)()()jwtRGfwbftwtbedt其中，实函数w(t)为是时窗函数，窗函数w(t)具有较强的衰减性，所以要精心选择窗函数。下面是一个短时傅立叶变换的代码程序functiontimefreq(x,Nw,window)%待分析信号，行向量，Nw时窗宽度subplot(2,2,1);plot(real(x));%描绘待分析信号X=fft(x);%快速傅里叶变换X=fftshift(X);%调整0频位置subplot(2,2,2);plot(abs(X));%描绘幅度谱Lap=Nw/2;%重叠宽度Tn=(length(x)-Lap)/(Nw-Lap);%计算分段数目nfft=2^ceil(log2(Nw));%做fft的点数TF=zeros(Tn,nfft);%时频矩阵fori=1:Tnif(strcmp(window,'rec'))Xw=x((i-1)*10+1:i*10+10);%加窗矩形处理elseif(strcmp(window,'Hamming'))Xw=x((i-1)*10+1:i*10+10).*Hamming(Nw)';%加hamming处理elseif(strcmp(window,'Blackman'))Xw=x((i-1)*10+1:i*10+10).*Blackman(Nw)';%加black处理elseif(strcmp(window,'Gauss'))Xw=x((i-1)*10+1:i*10+10).*Gauss(Nw)';%加Gauss处理elsereturn;endtemp=fft(Xw,nfft);%求ffttemp=fftshift(temp);%调整0频位置TF(i,:)=temp;%保存分段fft结果end%绘制时频分析结果subplot(2,2,3);fnew=((1:nfft)-nfft/2)/nfft;tnew=(1:Tn)*Lap;[F,T]=meshgrid(fnew,tnew);mesh(T,F,abs(TF));xlabel('n');ylabel('w');zlabel('Gf');subplot(2,2,4);contour(T,F,abs(TF));xlabel('n');ylabel('w');例子：clc;clear;N=400;x=zeros(1,N);T=0:N-1;x=exp(j*4*pi*(T/80).^2);figure(1);timefreq(x,20,’rec’);figure(2);timefreq(x,20,’Blackman’);