同步压缩小波变换（Synchrosqueezing Wavelet Transform, SWST）的 MATLAB 实现，用于一维信号的时频分布分析。结合连续小波变换（CWT）和瞬时频率重分配技术，提供高分辨率时频图。
一、核心代码
function [Tfr, t, f] = swst(signal, fs, wavelet, scales)
% 输入参数：
% signal: 输入信号（一维向量）
% fs: 采样频率（Hz）
% wavelet: 小波基函数（如 'morl', 'amor', 'db4'）
% scales: 尺度向量（建议使用对数尺度）
% 计算连续小波变换（CWT）
[cfs, freq] = cwt(signal, scales, wavelet, 'SamplingPeriod', 1/fs);
% 估计瞬时频率（基于相位导数）
phase = angle(cfs);
dphase = diff(phase, 1, 2); % 相位差分
inst_freq = cumsum(dphase, 2) ./ (2*pi*fs); % 累积积分得到瞬时频率
% 同步压缩操作
[Tfr, f] = wsstridge(cfs, inst_freq, freq); % 使用WSST算法压缩能量
% 可视化
figure;
imagesc(t, f, abs(Tfr));
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('Synchronized Compressed WST Time-Frequency Distribution');
colorbar;
end
% 示例调用
fs = 1000; % 采样率
t = 0:1/fs:1-1/fs; % 时间向量
signal = sin(2*pi*50*t) + 0.5*sin(2*pi*150*t + 0.3*t.^2); % 线性调频信号 + 噪声
scales = scal2frq('morl', 128, 1/fs); % Morlet小波尺度
swst(signal, fs, 'morl', scales);
二、关键步骤解析
1. 小波变换
- 函数：
cwt（MATLAB内置连续小波变换） - 参数：
scales：建议使用对数尺度（如logspace(1, 3, 128)）wavelet：推荐morl（复Morlet小波）或amor（自适应Morlet小波）
- 输出：小波系数矩阵
cfs和对应频率向量freq
2. 瞬时频率估计
相位导数法：通过计算小波系数相位的一阶差分，结合时间积分得到瞬时频率
公式：
finst(t)=2π1k=1∑KΔtΔϕ(t,k)其中 Δϕ为相位差分，Δt为时间步长
3. 同步压缩
- WSST算法：将小波系数能量沿瞬时频率方向压缩
- 函数：
wsstridge（需自定义实现或调用第三方库） - 优化：通过频率对齐和能量重分配提升分辨率
三、应用案例
1. 机械振动信号分析
% 生成含冲击的振动信号
fs = 5000; t = 0:1/fs:0.1;
signal = sin(2*pi*100*t) + 0.8*sin(2*pi*400*t) + 0.3*randn(size(t));
% 添加冲击成分
signal(0.02*fs:0.03*fs) = signal(0.02*fs:0.03*fs) + 2*sin(2*pi*2000*t(0.02*fs:0.03*fs));
% 执行SST分析
swst(signal, fs, 'amor', scal2frq('amor', 256, fs));
2. 生物医学信号处理
% 加载ECG信号（示例）
load('ecg_signal.mat');
% 去除基线漂移
baseline = movmean(ecg_signal, 500);
ecg_clean = ecg_signal - baseline;
% SST分析
swst(ecg_clean, 360, 'morl', scal2frq('morl', 128, 360));
四、结果可视化优化
% 高亮特定频率区域
hold on;
plot([0.2, 0.2], ylim, 'r--', 'LineWidth', 1.5); % 标记故障频率
% 添加时间-频率能量分布曲线
[~, idx] = max(abs(Tfr), [], 2);
plot(t, f(idx), 'w', 'LineWidth', 1.5);
五、扩展功能
1. 三维时频曲面图
[X, Y] = meshgrid(t, f);
Z = abs(Tfr);
surf(X, Y, Z, 'EdgeColor', 'none');
shading interp;
colormap(jet);
2. 动画展示频率演化
for i = 1:size(Tfr, 1)
imagesc(t(1:i), f, Tfr(1:i,:));
drawnow;
end
六、参考
- 核心文献
- Daubechies et al., Synchrosqueezing Transform for Instantaneous Frequency Analysis
- Thakur et al., Multisynchrosqueezing Transform for Multi-Component Signals
- 参考代码 ：同步压缩小波变换程序 www.youwenfan.com/contentalf/112998.html
- MATLAB工具箱
- Wavelet Toolbox（需安装）
- Signal Processing Toolbox
