时频分析matlab

模拟一段信号，分别进行小波变换、短时傅里叶变换、S变换、WV变换，获得其时频图。

% 不同分析方法对比
clear,clc;close all

% 产生chirp信号
fs = 4000;
t = 0:1/fs:1;
s1 = chirp(t,0,1,fs/2,'quadratic');
s2 = chirp(t,0,1,fs/2,'linear');
s3 = sin(2*pi*100*t);
s4 = sin(2*pi*200*t);
s = zeros(size(s1));  % 合成的非线性非平稳信号
s(400:450) = s3(400:450);      % 0.1s-0.1125s
s(650:800) = s3(650:800) + s4(650:800);  % 0.1625s-0.2s
s(800:1000) = s2(800:1000);    % 0.2s-0.25s
s(2000:2800) = s1(2000:2800);  % 0.5s-0.7s

figure(1)
subplot 511
plot(t,s1,'k');title('信号1')
subplot 512
plot(t,s2,'k');title('信号2')
subplot 513
plot(t,s3,'k');title('信号3')
subplot 514
plot(t,s4,'k');title('信号4')
set(0,'defaultfigurecolor','w')
subplot 515
plot(t,s,'k');title('合成信号')
set(0,'defaultfigurecolor','w')

%%信号加噪
SNR = 5;
r2 = randn(length(s),1);                      % 选择白噪声
r2 = r2/sqrt(sum(r2.^2))*sqrt(sum(s.^2)/(exp(SNR/10)));
r1 = r2+s;

% 连续小波变换时频图
wavename = 'cmor3-3';  % cmor是复Morlet小波，其中3－3表示Fb－Fc，Fb是带宽参数，Fc是小波中心频率。
totalscal = 256;
Fc = centfrq(wavename);  % 小波的中心频率
c = 2*Fc*totalscal;
scals = c./(1:totalscal);
f = scal2frq(scals,wavename,1/fs); % 将尺度转换为频率
coefs = cwt(s,scals,wavename); % 求连续小波系数
figure(2)
subplot 221
imagesc(t,f,abs(coefs)/max(max(abs(coefs))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('小波时频图');

% 短时傅里叶变换时频图
% figure
% spectrogram(s,256,250,256,fs);
% 时频分析工具箱里的短时傅里叶变换
f = 0:fs/2;
tfr = tfrstft(s');
tfr = tfr(1:floor(length(s)/2), :);
subplot 222
imagesc(t, f, abs(tfr)/max(max(abs(tfr))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('短时傅里叶变换时频图');

% st
[st_fre,st_times,st_frequencies] = st(s,4);
subplot 223
imagesc(t,f,abs(st_fre)/max(max(abs(st_fre))));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('广义S变换时频图');

% W-V
sig = s;
N=length(sig);
% 归一化
sig=(sig-mean(sig))/std(sig,1);
% 计算伪Wigner-Ville分布
sig = (sig-mean(sig))/std(sig,1);
sig = hilbert(sig);
[tfr0,t0,f0] = tfrpwv(sig');
subplot 224
imagesc(t',fliplr(f),abs(tfr0)/max(max(abs(tfr0))));
colorbar;
axis xy;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('W-V分布');

得到时频图：