%例程6－4  計算非平穩(wěn)隨機信號的WVD時頻譜
% e.g.6-4.m for example6-4;
 
clear;
N=500;      %length of the signal x.
t=0:N-1;
x=zeros(size(t));
%Generate two sine signal with different frequencies.
x(50:150)=cos(pi*(t(50:150)-50)/10);   
x(250:350)=cos(pi*(t(250:350)-250)/20);
subplot(221);
plot(t,x);
axis([ 0 500 -1 1]);
grid;
X=fft(x);
 %Shift zero-frequency component to center of spectrum
X=fftshift(X);        
subplot(222);
plot(abs(X));
axis([0 500 0 60 ]);
grid;
R=zeros(N,N);
for i=1:N-1
   M=min(i,N-1-i);
   for j=0:M
       R(i+1,j+1)=x(i+j+1)*conj(x(i-j+1));
   end
   for j=N-1:-1:N-M
       R(i+1,j+1)=conj(R(i+1,N-j+1));
   end
end
TF=zeros(N,N);
for i=0:N-1
    temp=fft(R(i+1,:));
    temp=fftshift(temp);
    TF(i+1,:)=temp;
end
f_new=(t-N/2)/N;
t_new=0:N-1;
[F,T]=meshgrid(f_new,t_new);  %generate F and T arrays for 3-D plots
subplot(223);
mesh(F,T,abs(TF));
subplot(224);
contour(abs(TF));
grid;