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

原创于 2025-10-22 16:07:19 发布 · 349 阅读

CC 4.0 BY-SA版权

文章标签：

基于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）：

二阶导数修正

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

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

频率对齐优化

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

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

抗噪增强处理

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

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

应用场景：

机械故障诊断

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

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

海洋波浪分析

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

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

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

优化建议：

GPU加速

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

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

多尺度分解

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

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