9y9oduf4.m

有关时频分析做三维图形的源程序

%%%%%%%%%%%
%用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存在着交叉干扰项，出现了一些干扰频率，尤其对实信号更加严重