MATLAB 实现同步压缩小波变换

同步压缩小波变换（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

---

六、参考文献与工具

1. 核心文献
   • Daubechies et al., Synchrosqueezing Transform for Instantaneous Frequency Analysis
   • Thakur et al., Multisynchrosqueezing Transform for Multi-Component Signals
2. 参考代码 ：同步压缩小波变换程序 www.youwenfan.com/contentcsf/112998.html
3. MATLAB工具箱
   • Wavelet Toolbox（需安装）
   • Signal Processing Toolbox