%例程6－5  計算非平穩隨機信號的WVD時頻譜
% e.g.6-5.m for example6-5;
 
clear;
k=5;T=4;
fc=k*T;
%Case 1 for fs=2fc 
fs=2*fc;
Ts=1/fs;
N=T/Ts;
t=0:N-1;
x=zeros(size(t));
x=exp(j*k*pi*(Ts*t).^2);
subplot(221);
plot(Ts*t,real(x));grid;
X=fft(x);
 %Shift zero-frequency component to center of spectrum
X=fftshift(X);
subplot(222);
plot((t-N/2)*fs/N,abs(X));
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)*fs/(2*N);
t_new=(0:N-1)*Ts;
 %generate F and T arrays for 3-D plots
[F,T]=meshgrid(f_new,t_new); 
subplot(223);
mesh(F,T,abs(TF));
subplot(224);
contour(F,T,abs(TF));
grid; 

%Case 2 for fs=4fc 
fs=4*fc;
Ts=1/fs;
N=T/Ts;
t=0:N-1;
x=zeros(size(t));
x=exp(j*k*pi*(Ts*t).^2);
subplot(221);
plot(Ts*t,real(x));grid;
X=fft(x);
 %Shift zero-frequency component to center of spectrum
X=fftshift(X);
subplot(222);
plot((t-N/2)*fs/N,abs(X));
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)*fs/(2*N);
t_new=(0:N-1)*Ts;
 %generate F and T arrays for 3-D plots
[F,T]=meshgrid(f_new,t_new); 
subplot(223);
mesh(F,T,abs(TF));
subplot(224);
contour(F,T,abs(TF));
grid;