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/contentcni/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