%% 写一个多分量LFM信号,用多阶次STFRFT做时频分析,并与STFT、STFRFT比较
clc;
clear all;
close all;
fs=1000;
L=3;
ts=1/fs;
t1=1*fs:L*fs;
s1=exp(j*pi*100*(t1/fs).^2);
s2=exp(j*pi*(-50)*(t1/fs).^2);
% s=zeros(9,length(t1));
% sum=zeros(9,length(t1));
% for i=1:9; %调频率
% k(i)=-20+(i-1)*5;
% s(i,:)=exp(j*pi*k(i)*t1.^2+j*2*pi*t1);
% sum=sum+s;
% end
% %sum=s(i,:)+sum;
% sum=reshape(sum,1,9*length(t1));
s=s1+s2;
figure;
% plot(real(s1),'r');
hold on;plot(real(s2),'g');
T=length(s)/fs;
%% STFT
N=length(s);
windowlength=32;
shift=windowlength/32;
n_frame=(N-windowlength)/shift+1;
NFFT=32;
wham=gausswin(windowlength);
TF=zeros(NFFT,n_frame);
for i=1:n_frame
n1=shift*(i-1)+1;
n2=shift*(i-1)+windowlength;
sig=s(n1:n2);
sf=fft(sig,NFFT);
TF(:,i)=fftshift(abs(sf));
end
figure;
colormap(jet(256))
imagesc([0,T],[-fs/2,fs/2],abs(TF))
xlabel('time (s)')
ylabel('frequency (Hz)')
title('STFT')
axis xy
colorbar
clim=get(gca,'Clim');
set(gca,'Clim',clim(2)+[-20 0]);
drawnow
%
%
%% STFRFT
N=length(s);
windowlength=32;
shift=windowlength/32;
n_frame=(N-windowlength)/shift+1;
NFFT=32;
wham=hanning(windowlength); %rectwin gausswin
TF=zeros(NFFT,n_frame);
TF_sum=zeros(NFFT,n_frame);
p=0.5:0.5:1.5; %0.5:0.1:1.5;
for h=1:length(p)
for n=1:n_frame
n1=shift*(n-1)+1;
n2=shift*(n-1)+windowlength;
sig=s(n1:n2).*wham.';
sf=frft(sig,p(h));
TF(:,n)=fftshift(abs(sf));
end
TF_sum=TF_sum+TF;
figure;
colormap(jet(256))
imagesc([0,T],[-fs/2,fs/2],20*log10(fftshift(abs(TF),1)+eps))
xlabel('time (s)')
ylabel('frequency (Hz)')
title('STFRFT')
axis xy
clim=get(gca,'Clim');
set(gca,'Clim',clim(2)+[-20 0]);
drawnow
end
figure;
colormap(jet(256))
imagesc([0,T],[-fs/2,fs/2],20*log10(fftshift(abs(TF_sum),1)+eps))
xlabel('time (s)')
ylabel('frequency (Hz)')
title('STFRFT综合分析')
axis xy
clim=get(gca,'Clim');
set(gca,'Clim',clim(2)+[-20 0]);
drawnow
评论8