%%%%%%%%%%%
%用WD變換解析信號x=exp(j*pi*k*t.^2)
clear,close all
k=4;T=4;
fc=k*T;fs=4*fc;%采樣頻率大于4倍的信號寬度
Ts=1/fs;N=T/Ts;
x=zeros(1,N);
t=0:N-1;
x=exp(j*k*pi*(t*Ts).^2);
subplot(221),plot(t*Ts,real(x));
X=fftshift(fft(x));
subplot(222),plot((t-N/2)*fs/N,abs(X))
R=zeros(N);
for n=0:N-1
    M=min(n,N-1-n);
    for k=0:M
        R(n+1,k+1)=x(n+k+1)*conj(x(n-k+1));
    end
    for k=N-1:-1:N-M
        R(n+1,k+1)=conj(R(n+1,N-k+1));
    end
end
TF=zeros(N);
for n=0:N-1
   temp=fftshift(fft(R(n+1,:)));
   TF(n+1,:)=temp;
end
fnew=(t-N/2)*fs/2/N;
tnew=(0:N-1)*Ts;
[F,T]=meshgrid(fnew,tnew);
subplot(223),mesh(F,T,abs(TF))
subplot(224),contour(F,T,abs(TF))

%實信號數字信號N=200，在(50,100)和(120,180)處分別有兩個頻率不同的正弦信號
clear
figure(2)
N=200;t=0:N-1;
x=zeros(1,N);
x(50:100)=cos(pi*(t(50:100)-50)/10);
x(120:180)=cos(pi*(t(120:180)-300)/5);
subplot(221),plot(x)
X=fft(x);
X=fftshift(X);
subplot(222)
plot((t-N/2)/N,abs(X))
R=zeros(N);
for n=0:N-1
    M=min(n,N-1-n);
    for k=0:M
        R(n+1,k+1)=x(n+k+1)*conj(x(n-k+1));
    end
    for k=N-1:-1:N-M
        R(n+1,k+1)=conj(R(n+1,N-k+1));
    end
end
TF=zeros(N);
for n=0:N-1
   temp=fftshift(fft(R(n+1,:)));
   TF(n+1,:)=temp;
end
fnew=(t-N/2)/N;%????
tnew=0:N-1;
[F,T]=meshgrid(fnew,tnew);
subplot(223),mesh(F,T,abs(TF))
subplot(224),contour(F,T,abs(TF))

%解析信號
clear
figure(3)
N=200;t=0:N-1;
x=zeros(1,N);
x(50:100)=exp(j*pi*(t(50:100)-50)/10);
x(120:180)=exp(j*pi*(t(120:180)-300)/5);
subplot(221),plot(real(x))
X=fft(x);
X=fftshift(X);
subplot(222)
plot((t-N/2)/N,abs(X))
R=zeros(N);
for n=0:N-1
    M=min(n,N-1-n);
    for k=0:M
        R(n+1,k+1)=x(n+k+1)*conj(x(n-k+1));
    end
    for k=N-1:-1:N-M
        R(n+1,k+1)=conj(R(n+1,N-k+1));
    end
end
TF=zeros(N);
for n=0:N-1
   temp=fftshift(fft(R(n+1,:)));
   TF(n+1,:)=temp;
end
fnew=(t-N/2)/N;%????
tnew=0:N-1;
[F,T]=meshgrid(fnew,tnew);
subplot(223),mesh(F,T,abs(TF))
subplot(224),contour(F,T,abs(TF))
%%WD存在著交叉干擾項，出現了一些干擾頻率，尤其對實信號更加嚴重