时频图程序(小波时频,短时傅里叶变换,s变换) MATLAB程序
先搞个模拟信号热热身:
fs = 1000; % 采样率别太小气 t = 0:1/fs:1; f1 = 20; f2 = 100; signal = sin(2*pi*f1*t).*(t<0.5) + sin(2*pi*f2*t).*(t>=0.5); % 分段信号 noise = 0.5*randn(size(t)); % 加点佐料 x = signal + noise;这信号前0.5秒20Hz,后0.5秒突变到100Hz,非常适合拿来折腾时频分析。
STFT暴力美学
短时傅里叶变换属于老牌劲旅,代码简单粗暴:
window = hamming(128); % 窗长别瞎搞 noverlap = 120; % 重叠量 nfft = 256; spectrogram(x,window,noverlap,nfft,fs,'yaxis') colorbar off; % 强迫症必改关键参数就仨:窗长决定时频分辨率平衡(鱼与熊掌问题),noverlap影响滑动流畅度。有个坑要注意——窗太长会导致时间分辨率扑街,特别是像我们这个突变信号,可能捕捉不到精确的跳变时刻。
小波变换的千层套路
连续小波变换(CWT)更擅长捕捉突变:
[cfs,frq] = cwt(x,'amor',fs); % 用墨西哥帽小波 contour(t,frq,abs(cfs)) set(gca,'YScale','log') % 对数坐标更直观 colorbar title('小波时频图')这里'amor'指的是解析Morlet小波。小波变换的尺度选择是门玄学,MATLAB的cwt函数已经帮我们自动计算了最佳尺度。注意颜色映射的绝对值代表能量强度,那些突然出现的亮斑就是信号突变的证据。
S变换的均衡术
ST算是STFT和小波的私生子:
[st_matrix, st_t, st_f] = st(x,fs); % 自建函数见下文 imagesc(st_t, st_f, abs(st_matrix)) axis xy; colormap(jet) % 反转Y轴 function [st_matrix, t, f] = st(x,fs) N = length(x); t = (0:N-1)/fs; f = (0:N/2-1)*(fs/N); st_matrix = zeros(length(f),N); for k = 1:length(f) sigma = 1/(2*pi*f(k)+eps); % 防止除零 window = exp(-(t - t(end)/2).^2 / (2*sigma^2)); st = fft(x .* window); st_matrix(k,:) = abs(st(1:N/2)); end end这个自编S变换实现用了高斯窗,并且窗宽随频率自适应变化。注意sigma的计算加了eps防止频率为零时报错。相比STFT固定窗长,S变换在低频区用宽窗(频率分辨率高),高频区用窄窗(时间分辨率高),算是个端水大师。
实战对比
把三张时频图放一起看:
- STFT在20Hz处拖尾严重(窗太长导致)
- 小波在100Hz跳变时刻有更尖锐的过渡
- S变换在低频区频率定位更准确
个人经验:处理振动信号首选STFT,检测信号突变用小波,搞地震信号分析试试S变换。但别死磕理论,多调参才是王道——有时候把窗长改个10个点,效果立竿见影。