mrbz2kbr.m

时频分析工具

%%短时傅立叶变换
%解析信号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))

%解析信号x=cos(pi*k*t.^2)
clear
figure(2)
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=cos(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))
%对于解析信号只有正频率，实信号还有负频率


%数字信号N=400，在(50,100)和(300,360)处分别有两个频率不同的正弦信号
clear
fs=1;Ts=1/fs;
figure(3)
N=400;t=0:N-1;
x=zeros(1,N);
x(50:100)=cos(pi*(t(50:100)-50)/10);
x(300:360)=cos(pi*(t(300:360)-300)/5);
subplot(221),plot(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))

%在信号的频谱中可以观测到有两个不同频率的信号，但不能得出他们的持续范围，而从信号的时频平面不仅
%可以看出这两种信号的频率，还可以看出他们的持续时间