短时傅立叶变换试验

为了克服傅立叶变换的时频局部化方面的不足，也是为了对时域信号作局部分析，D.Gabor 于1946年提出了窗口傅立叶变换（简记为WFT ）。

WFT 的公式形式

()(,)()()j w t R G f w b f t w t b e d t -=-?()(,)()()j w t R G f w b f t w t b e d t

-=-? 其中，实函数w(t)为是时窗函数，窗函数w(t)具有较强的衰减性，所以要精心选择窗函数。

下面是一个短时傅立叶变换的代码程序

function timefreq(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);%时频矩阵

for i=1:Tn

if(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 处理 else

return;

end

temp=fft(Xw,nfft);%求fft

temp=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’);