没有直接的函数画,需要自己编成。短时傅立叶变换,gabor变换和WD变换都可以
%%短时傅立叶变换
%解析信号x=exp(j*pi*k*t.^2)
clear,clc,close all
figure(1)
k=4;T=5;
fc=k*T;
fs=3*fc;
Ts=1/fs;
N=T/Ts;
x=zeros(1,N);
t=0:N-1;
x=exp(j*k*pi*(t*Ts).^2);
% x=awgn(x,-3,'measured');
subplot(221)
plot(t*Ts,real(x))
X=fft(x);
X=fftshift(X);
subplot(222)
plot((t-N/2)*fs/N,abs(X))
Nw=20;
L=Nw/2;
Tn=(N-Nw)/L+1;
nfft=32;
TF=zeros(Tn,nfft);
for i=1:Tn
xw=x((i-1)*10+1:i*10+10);
temp=fft(xw,nfft);
temp=fftshift(temp);
TF(i,:)=temp;
end
subplot(223)
fnew=((1:nfft)-nfft/2)*fs/nfft;
tnew=(1:Tn)*L*Ts;
[F,T]=meshgrid(fnew,tnew);
mesh(F,T,abs(TF))
subplot(224)
contour(F,T,abs(TF))
%用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))
[ 本帖最后由 yejet 于 2006-7-9 17:40 编辑 ] |