基于MATLAB的二阶同步挤压小波变换（WSST2）实现

function [Tfr, t, f] = wsst2(signal, fs, wavelet, scales)
% 二阶同步挤压小波变换（WSST2）MATLAB实现
% 输入参数：
% signal: 输入信号（一维向量）
% fs: 采样频率（Hz）
% wavelet: 小波基函数（'morl','amor','db4'等）
% scales: 尺度向量（建议使用对数尺度）
% 输出参数：
% Tfr: 时频能量矩阵
% t: 时间向量
% f: 频率向量

%% 步骤1：计算连续小波变换（CWT）
[cfs, freq] = cwt(signal, scales, wavelet, 'SamplingPeriod', 1/fs);

%% 步骤2：计算相位导数（一阶和二阶）
phase = angle(cfs);
dphase = diff(phase, 1, 2);      % 一阶导数
d2phase = diff(dphase, 1, 2);    % 二阶导数

%% 步骤3：瞬时频率估计（二阶修正）
inst_freq = cumsum(dphase, 2) ./ (2*pi*fs);
inst_freq2 = inst_freq + ...     % 二阶修正项
    (d2phase .* cumsum(dphase.^2, 2)) ./ (2*(2*pi*fs).^3);

%% 步骤4：频率轴准备
[~, idx] = sort(freq, 2);
f = freq(idx(:,1), :);

%% 步骤5：二阶同步挤压操作
Tfr = zeros(size(cfs));
for i = 1:size(cfs,1)
    for j = 1:size(cfs,2)
        % 二阶频率对齐
        delta_b = 1/(2*pi*fs) * ( ...
            (d2phase(i,j) * (inst_freq(i,j) - f(j)) ) / ...
            ( (inst_freq(i,j) - f(j)).^2 + eps ) );
        idx2 = find(f >= inst_freq2(i,j), 1);
        Tfr(idx2,j) = Tfr(idx2,j) + abs(cfs(i,j))^2 * exp(-delta_b^2);
    end
end

%% 步骤6：时频图可视化
figure;
imagesc(t, f, Tfr);
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('WSST2时频谱');
colorbar;

end

%% 测试案例
if __name__ == "__main__"
    % 生成测试信号
    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);
    
    % 执行WSST2分析
    [Tfr, t, f] = wsst2(signal, fs, 'morl', scales);
    
    % 叠加原始信号时频图
    figure;
    subplot(2,1,1);
    imagesc(t, f, Tfr);
    title('WSST2时频谱');
    subplot(2,1,2);
    spectrogram(signal, 256, 250, 256, fs);
    title('STFT对比');
end

关键改进说明（相比传统SST）：

1. 二阶导数修正

   通过计算相位二阶导数（d2phase），改进瞬时频率估计精度：

   d2phase = diff(dphase, 1, 2);  % 二阶导数计算
   inst_freq2 = inst_freq + ...   % 二阶修正公式
       (d2phase .* cumsum(dphase.^2, 2)) ./ (2*(2*pi*fs).^3);
   

2. 频率对齐优化

   采用指数衰减函数进行二次频率对齐，提升能量集中效果：

   delta_b = 1/(2*pi*fs) * ( ... 
       (d2phase * (inst_freq - f) ) / ... 
       ( (inst_freq - f).^2 + eps ) );
   

3. 抗噪增强处理

   在能量重分配时引入高斯衰减因子：

   Tfr(idx2,j) = Tfr(idx2,j) + abs(cfs(i,j))^2 * exp(-delta_b^2);
   

应用场景：

1. 机械故障诊断

   分析轴承振动信号中的故障特征频率（如内圈损伤的3倍频）：

   % 加载振动信号
   load('bearing_fault.mat');
   [Tfr, t, f] = wsst2(vibration_signal, 10000, 'amor', scal2frq('amor', 512, 10000));
   

2. 海洋波浪分析

   提取畸形波的瞬时频率演化特征：

   % 加载海浪数据
   load('draupner_wave.mat');
   [Tfr, t, f] = wsst2(wave_height, 2, 'morl', scal2frq('morl', 128, 2));
   

参考代码 同步挤压小波变换的扩展 www.youwenfan.com/contentcni/65037.html

优化建议：

1. GPU加速

   对大规模数据使用gpuArray加速计算：

   cfs_gpu = gpuArray(cfs);
   Tfr_gpu = wsst2_gpu(cfs_gpu, ...);  % 自定义GPU版本函数
   

2. 多尺度分解

   结合EMD分解实现多尺度WSST2：

   imf = emd(signal);
   for i = 1:size(imf,1)
       [Tfr(:,:,i), t, f] = wsst2(imf(i,:), fs, 'morl', scales);
   end