www.pudn.com > fft_WD.rar > fft_WD.m

%利用FFT计算WVD分布(非平稳随机信号分析与处理第45页) 
clc 
clear all 
 
NN=100; 
 
w1=0.001; 
a1=0.0005; 
w2=0.001; 
a2=-0.005; 
j=sqrt(-1); 
nn=-3*NN-1:1:3*NN-1;              %构造信号 
%n=nn-3*NN; 
x=exp(j*(2*pi*(w1*nn+1/2*a1*nn.^2)));   %s1=ss1(nn-3*NN) 
 
%ss1c=conj(ss1); 
L=64;      %加矩形窗 
N=2*L; 
w=zeros(N); 
w(1)=0; 
l=-(L-1):1:(L-1); 
w(2:length(l))=1;     %2L-1窗 
%s1=ss1(3*NN-(L-1)+1:3*NN+(L-1)+1); 
Wx=zeros(N,N); 
G=zeros(N,N); 
GG=zeros(N,N);  %conj(G) 
f=zeros(N,N); 
 
for n=0:N-1       %  n=nn-3*NN 
%     l=-(L-1):1:(L-1) 
%     G=w.*x(3*NN+n+l)  %G(n,l)=w(l).*x(n+1) 
 %for 
 l=0:1:L-1; 
% ll=L+l; 
% nn=3*NN+n; 
G(n+1,l+1)=w(L+1+l).*x(3*NN+n+L+1+l); 
GG(n+1,l+1)=conj(w(L+1-l).*x(3*NN+n+L+1-l)); 
f(n+1,l+1)=G(n+1,l+1).*GG(n+1,l+1); 
  
l=L:1:2*L-1; 
G(n+1,l+1)=w(L+1+l-2*L).*x(3*NN+n+L+1+l-2*L); 
GG(n+1,l+1)=conj(w(L+1-l+2*L).*x(3*NN+n+L+1-l+2*L)); 
f(n+1,l+1)=G(n+1,l+1).*GG(n+1,l+1); 
 
end 
 
F=f';           %FFT按列计算矩阵的傅立叶变换 
Wx=fft(F); 
Wx=Wx'; 
n=0:N-1; 
k=0:N-1; 
 
figure; 
grid; 
mesh(n,k,Wx); 
 
figure; 
contour(Wx);