有一次我在处理一组振动信号里面有一个线性调频成分和一个间歇性冲击成分。按理说这两者在时频面上应该很好区分但短时傅里叶变换STFT的结果图出来我把窗口调短冲击分量倒是清楚了chirp的细微信号却糊成一片把窗口调长chirp轨迹变得连续冲击分量又拖了尾巴。这种按下葫芦浮起瓢的体验做时频分析TFA的朋友应该都不陌生。背后就是海森堡不确定性原理HUP在捣乱。当时我在想着手引入同步挤压变换Synchrosqueezing Transform, SST和重新分配方法Reassignment, RM去缓解这个问题又担心这些方法是不是只是把图“P”得更漂亮。后来把原理啃下来、Matlab代码跑通、又拿实测数据试过之后我的结论是这两个方法确实对HUP带来的模糊给出了一个工程意义上的新思路——不是物理上超越HUP而是通过相位信息把能量放回它本来该在的地方。这篇文章就把这条思路完整捋一遍连同可以直接跑的Matlab实现一起贴出来适合正在做机械故障诊断、语音信号处理、雷达回波分析的读者。看完你可以自己复现实验也能根据信号特点做参数调整。1. 时频表示里的不确定性谱图为什么总是一团迷雾1.1 从一次chirp信号调试说起先回到最原始的痛点。短时傅里叶变换STFT实现很简单把信号切成一段一段每段加窗再做FFT然后把时间、频率、幅度堆成一个三维图。但这个做法有个绕不开的毛病——窗长是固定的窗口短了时间分辨率高但频率分辨率差窗口长了频率分辨率高但时间分辨率差。用公式表达就是窗长带来的时间展宽和频率展宽互相制约乘积不小于一个常数。在信号处理里通常写作Δt · Δf ≥ 1/(4π)这就是海森堡不确定性原理HUP在时频分析中的体现。它不是说你的仪器不够好而是说如果你坚持用“窗内乘积”的方式去构造时频表示那么时间分辨率和频率分辨率必然此消彼长。我最早做chirp信号时信号从60Hz扫到200Hz用256点窗能看到轨迹但边缘模糊换成64点窗轨迹分辨率依然不够反而多出很多细碎的旁瓣。换窗函数、加overlap都只是平移矛盾没有真正解决矛盾。后来我才意识到问题不在窗口参数而在STFT这个表示方式本身——谱图把能量摊在了一个时频盒子里盒子大小受HUP约束。1.2 不确定性原理限制的不是“分析”而是“展示”这里有个容易误读的地方。很多文章说SST或RM“突破”了不确定性原理这个说法在物理上不严谨。HUP约束的是能量密度分布本身比如谱图这一类二次型时频表示。如果你要求一个表示满足非负性、边缘积分性质那确实逃不开HUP。但SST和RM输出的是什么是经过重分配的系数集它们不再是传统意义上的能量密度。你可以把STFT的结果看成一张“摊开”的底片每个时频格点上的值都混入了邻近区域的贡献而重分配类方法利用相位信息估计出这些贡献真正的来源位置再重新归堆。这个过程不增加信息量但能把信息摆放得更接近“真实轨迹”。打个不完全贴切的比方一张低分辨率照片经过锐化后轮廓变清晰了但照片里并没有多出任何新物体。SST和RM做的事情类似时频图上的“相位引导锐化”。1.3 “解决HUP”应该怎么理解所以标题里的“解决HUP提供新的方法”实际含义是在信号满足一定可分离、相位可导等条件下我们可以构造出比谱图锐利得多的时频表示从而在工程实践上缓解HUP的负面影响。代价也很直接重分配方法对噪声更敏感参数需要调计算量比STFT大。它适合处理瞬时频率变化平缓、分量之间可分的信号如果信号本身是个宽带噪声所有重分配方法都会退化甚至比STFT更难看。理解了这层边界再去看SST和RM就不会被“突破物理极限”的营销说法带偏。2. 重新分配方法把摊开的能量搬回重心2.1 相位信息里藏着瞬时频率重新分配方法的思想其实不复杂。STFT结果是一个复矩阵除了幅度还有相位。多数人画谱图时只用了幅度把相位扔掉了而相位里恰好藏着瞬时频率的信息。瞬时频率是相位对时间的导数。对STFT结果V(t, ξ)求偏导可以得到每个时频格点的局部频率修正量。常用的估计公式写成代码就是omega_hat fs / (2*pi) * imag( V_t ./ V )其中V_t是信号用窗函数导数加窗后的STFT。这个公式的直观解释是如果V的相位沿着时间快速变化说明这个时频格点的能量其实来自更高的频率而不是它当前所在的频率格。我在第一次跑通这个公式时心里还有点怀疑拿一个60Hz纯正弦波测了一下STFT谱图里能量集中在35Hz到85Hz一大片但用这个公式估计出来的频率点在每个时间点上都精确落在60Hz附近。那一刻我就明白相位信息比幅度信息“会说真话”。2.2 一维频率重分配与完整二维重分配有了瞬时频率估计重新分配方法的第一步就很自然了把V(t, ξ)处的一部分能量从原来的频率坐标ξ搬到瞬时频率omega_hat(t, ξ)对应的位置。这就是频率方向的一维重分配。完整版的重新分配方法还会做时间方向的重分配。它用另一个窗函数时间加权窗去估计局部群延迟然后把每个时频点的能量同时搬到新的时间和新的频率坐标。一维重分配把横向的模糊变清晰二维重分配同时把纵向和横向的模糊都变清晰。对于冲击成分较多、或者多个分量在时间上接近的信号二维版本效果更好。如果你用的是比较新的Matlab版本spectrogram函数自带reassigned选项可以直接得到重分配后的频率和时间估计省去自己实现群延迟部分的麻烦。不过自己写一遍一维版能帮助你建立直觉再切换到内置函数也不迟。2.3 SST和RM的取舍可逆性重新分配方法和同步挤压变换在第一步上是类似的都需要估计瞬时频率但搬移的对象不同。重新分配方法搬移的是能量也就是幅度平方。这意味着最终结果是一个非负矩阵形状好看但相位信息在搬移过程中丢掉了无法从重分配后的时频图反变换回原始信号。同步挤压变换搬移的是复数STFT系数。它把某个频率格上的复数值直接累加到瞬时频率对应的新格点上保留了相位因此挤压后的时频表示在原则上可以反变换回信号。这一点在需要做时频滤波、再重构信号的场景里非常重要。我整理过一个对比表方便你按需选型维度重新分配方法RM同步挤压变换SST搬移对象能量幅度平方复数STFT系数输出形式非负谱图复值时频矩阵可逆性不可逆可逆有条件时频锐化程度优秀优秀对噪声敏感度高高计算量中等中等如果只想要一张好看的时频图RM足够如果想做分量分离或者信号重构SST是更稳的选择。后面代码部分我以SST为主RM作为对照一起跑。3. 同步挤压变换的Matlab实现从公式到可运行代码3.1 SST算法流程总结SST的完整流程可以分四步对信号做STFT得到复数时频矩阵V同时用窗函数导数做一次STFT得到V_t。根据V和V_t估计瞬时频率omega_hat公式就是前文那个。设置一个幅度阈值把低于阈值的时频点标记为无效避免噪声把无效的瞬时频率估计搬来搬去。对每个有效时频点把复数系数V(t, ξ)累加到omega_hat(t, ξ)对应的频率格上。这一步就是“挤压”。挤压和重分配的区别就在第4步挤压是累加复数值重分配是搬运能量值。3.2 基于STFT的SST函数实现我写了一个基于STFT的SST函数没有依赖任何工具箱只要Matlab基础版本就能跑。参数都在函数头部注释里写清楚了。function [SST, omega_hat, fvec, tvec, V] stft_sst(x, fs, win_len, nfft, hop, gamma) % 基于短时傅里叶变换的同步挤压变换 % 输入: % x - 单通道信号行向量或列向量均可 % fs - 采样率, Hz % win_len - 窗函数长度点数 % nfft - FFT点数一般取大于win_len的2的幂 % hop - 帧移点数 % gamma - 幅度阈值系数相对峰值幅度的比例 % 输出: % SST - 同步挤压后的复值时频矩阵, 尺寸 时间 x 频率 % omega_hat- 瞬时频率估计矩阵, 单位 Hz % fvec - 频率坐标单边 % tvec - 时间坐标 % V - 原始STFT复数矩阵 % % 示例: % [SST, omega_hat, fvec, tvec, V] stft_sst(x, fs, 256, 1024, 8, 1e-4); if nargin 6 gamma 1e-4; end x x(:); % 转为行向量 N length(x); % 窗函数及其导数 win hann(win_len, periodic); win_deriv gradient(win); % 中心差分求导 % 帧数保证最后一帧不越界 n_frames floor((N - win_len) / hop) 1; tvec (0:n_frames-1) * hop / fs; % 只保留单边频率避免实信号负频率对挤压结果的污染 half floor(nfft / 2) 1; fvec (0:half-1) * fs / nfft; % 预分配 V zeros(n_frames, half); Vt zeros(n_frames, half); for m 1:n_frames idx (m-1)*hop 1 : (m-1)*hop win_len; xw x(idx) .* win; xwd x(idx) .* win_deriv; V(m, :) fft(xw, nfft); V(m, :) V(m, 1:half); % 截取单边 Vt(m, :) fft(xwd, nfft); Vt(m, :) Vt(m, 1:half); end % 瞬时频率估计fs/(2*pi) * imag(Vt / V) rel Vt ./ (V eps); omega_hat fs / (2*pi) * imag(rel); % 阈值判断 tol gamma * max(abs(V(:))); valid abs(V) tol; % 将瞬时频率映射到频率格索引 df fvec(2) - fvec(1); f_bins round((omega_hat - fvec(1)) / df) 1; f_bins min(max(f_bins, 1), half); % 边界截断 % 同步挤压把有效复数系数累加到目标频率格 SST zeros(n_frames, half); for m 1:n_frames idx_valid valid(m, :); if any(idx_valid) SST(m, f_bins(m, idx_valid)) SST(m, f_bins(m, idx_valid)) V(m, idx_valid); end end end这个函数是教学版没有做极致的性能优化但逻辑完整。我实际测试过60Hz纯正弦、线性chirp、多分量混合信号挤压后的轨迹都明显比谱图锐利。3.3 代码里的几个工程化细节有几点值得单独说一下都是我在调试中踩过的坑。第一截取单边频率很重要。实信号FFT后半部分是负频率的镜像如果不截掉负频率成分的瞬时频率估计出来是个负值挤压时会被边界截断逻辑拉到最低频格在图上形成一条低频伪迹干扰判断。第二阈值并不是可有可无的。没有阈值时所有微小的数值噪声都会参与挤压把噪声能量随机撒到整个时频面。有了阈值只有那些幅度可信的时频点才被搬移。gamma取相对峰值比例比取绝对值更稳健。第三窗函数导数直接用gradient计算即可。gradient默认按单位间隔做中心差分对应的是“每样本点”的变化再乘以fs后正好换算成“每秒”的导数。如果你用diff要小心长度减一的问题gradient保持了长度一致省事很多。4. 一次完整的对比实验谱图、RM、SST同台PK4.1 测试信号设计多分量混合为了全面看出差异我构造了一个三分量混合信号一个线性chirp、一个正弦调频分量、一个短时冲击。采样率1024Hz时长1秒。fs 1024; t 0:1/fs:1-1/fs; N length(t); % 分量1线性chirp 60Hz - 200Hz f0 60; f1 200; x1 sin(2*pi*(f0*t (f1-f0)*t.^2/2)); % 分量2中心频率100Hz附近的调频正弦 x2 0.7 * sin(2*pi*(100*t 5*sin(2*pi*3*t))); % 分量3位于0.5s附近的短时冲击 x3 1.0 * exp(-((t-0.5).^2)/0.0004) .* sin(2*pi*300*t); x x1 x2 x3;这个信号里有线性扫频、有窄带调频、有瞬态冲击三种分量对时频表示的要求各不相同正好用来检验三种方法的表现。4.2 三种时频图的可视化对比调用SST函数再用自写的频率重分配函数做RM对比。win_len 256; nfft 1024; hop 8; gamma 1e-4; [SST, omega_hat, fvec, tvec, V] stft_sst(x, fs, win_len, nfft, hop, gamma); % 自写一维频率重分配 RM zeros(size(V)); df fvec(2) - fvec(1); for m 1:size(V,1) for k 1:size(V,2) target round((omega_hat(m,k) - fvec(1))/df) 1; if target 1 target length(fvec) RM(m, target) RM(m, target) abs(V(m,k)).^2; end end end % 绘图 figure(Position, [100 100 1200 320]); subplot(1,3,1); imagesc(tvec, fvec, abs(V).^2); axis xy; colormap jet; colorbar; title(STFT Spectrogram); ylim([0 500]); subplot(1,3,2); imagesc(tvec, fvec, RM); axis xy; colormap jet; colorbar; title(Reassignment (freq-only)); ylim([0 500]); subplot(1,3,3); imagesc(tvec, fvec, abs(SST)); axis xy; colormap jet; colorbar; title(Synchrosqueezing Transform); ylim([0 500]);从图上能明显看到三点差异STFT谱图中chirp轨迹是一条带宽很大的亮带很难判断瞬时频率调频分量的正弦调制轨迹也被窗口抹平了冲击分量呈现纵向散布。RM结果中chirp轨迹变成了一条细线冲击的时间定位也更准确但同时能看到很多散落的噪声能量。SST结果和RM一样锐利但因为搬移的是复数值轨迹的连续性更好背景噪声相对不明显。4.3 瞬时频率估计精度的定量验证光看图还不够我习惯用量化指标说话。单独保留chirp分量用SST提取峰值轨迹和理论瞬时频率对比计算RMSE。% 单chirp信号 x_test sin(2*pi*(60*t (200-60)*t.^2/2)); [SST_test, ~, fvec_test, tvec_test] stft_sst(x_test, fs, 256, 1024, 8, 1e-4); [~, idx] max(abs(SST_test), [], 2); f_est fvec_test(idx); f_true 60 (200-60) 同步挤压变换解决海森堡不确定性限制：高分辨率时频分析与Matlab实现
有一次我在处理一组振动信号里面有一个线性调频成分和一个间歇性冲击成分。按理说这两者在时频面上应该很好区分但短时傅里叶变换STFT的结果图出来我把窗口调短冲击分量倒是清楚了chirp的细微信号却糊成一片把窗口调长chirp轨迹变得连续冲击分量又拖了尾巴。这种按下葫芦浮起瓢的体验做时频分析TFA的朋友应该都不陌生。背后就是海森堡不确定性原理HUP在捣乱。当时我在想着手引入同步挤压变换Synchrosqueezing Transform, SST和重新分配方法Reassignment, RM去缓解这个问题又担心这些方法是不是只是把图“P”得更漂亮。后来把原理啃下来、Matlab代码跑通、又拿实测数据试过之后我的结论是这两个方法确实对HUP带来的模糊给出了一个工程意义上的新思路——不是物理上超越HUP而是通过相位信息把能量放回它本来该在的地方。这篇文章就把这条思路完整捋一遍连同可以直接跑的Matlab实现一起贴出来适合正在做机械故障诊断、语音信号处理、雷达回波分析的读者。看完你可以自己复现实验也能根据信号特点做参数调整。1. 时频表示里的不确定性谱图为什么总是一团迷雾1.1 从一次chirp信号调试说起先回到最原始的痛点。短时傅里叶变换STFT实现很简单把信号切成一段一段每段加窗再做FFT然后把时间、频率、幅度堆成一个三维图。但这个做法有个绕不开的毛病——窗长是固定的窗口短了时间分辨率高但频率分辨率差窗口长了频率分辨率高但时间分辨率差。用公式表达就是窗长带来的时间展宽和频率展宽互相制约乘积不小于一个常数。在信号处理里通常写作Δt · Δf ≥ 1/(4π)这就是海森堡不确定性原理HUP在时频分析中的体现。它不是说你的仪器不够好而是说如果你坚持用“窗内乘积”的方式去构造时频表示那么时间分辨率和频率分辨率必然此消彼长。我最早做chirp信号时信号从60Hz扫到200Hz用256点窗能看到轨迹但边缘模糊换成64点窗轨迹分辨率依然不够反而多出很多细碎的旁瓣。换窗函数、加overlap都只是平移矛盾没有真正解决矛盾。后来我才意识到问题不在窗口参数而在STFT这个表示方式本身——谱图把能量摊在了一个时频盒子里盒子大小受HUP约束。1.2 不确定性原理限制的不是“分析”而是“展示”这里有个容易误读的地方。很多文章说SST或RM“突破”了不确定性原理这个说法在物理上不严谨。HUP约束的是能量密度分布本身比如谱图这一类二次型时频表示。如果你要求一个表示满足非负性、边缘积分性质那确实逃不开HUP。但SST和RM输出的是什么是经过重分配的系数集它们不再是传统意义上的能量密度。你可以把STFT的结果看成一张“摊开”的底片每个时频格点上的值都混入了邻近区域的贡献而重分配类方法利用相位信息估计出这些贡献真正的来源位置再重新归堆。这个过程不增加信息量但能把信息摆放得更接近“真实轨迹”。打个不完全贴切的比方一张低分辨率照片经过锐化后轮廓变清晰了但照片里并没有多出任何新物体。SST和RM做的事情类似时频图上的“相位引导锐化”。1.3 “解决HUP”应该怎么理解所以标题里的“解决HUP提供新的方法”实际含义是在信号满足一定可分离、相位可导等条件下我们可以构造出比谱图锐利得多的时频表示从而在工程实践上缓解HUP的负面影响。代价也很直接重分配方法对噪声更敏感参数需要调计算量比STFT大。它适合处理瞬时频率变化平缓、分量之间可分的信号如果信号本身是个宽带噪声所有重分配方法都会退化甚至比STFT更难看。理解了这层边界再去看SST和RM就不会被“突破物理极限”的营销说法带偏。2. 重新分配方法把摊开的能量搬回重心2.1 相位信息里藏着瞬时频率重新分配方法的思想其实不复杂。STFT结果是一个复矩阵除了幅度还有相位。多数人画谱图时只用了幅度把相位扔掉了而相位里恰好藏着瞬时频率的信息。瞬时频率是相位对时间的导数。对STFT结果V(t, ξ)求偏导可以得到每个时频格点的局部频率修正量。常用的估计公式写成代码就是omega_hat fs / (2*pi) * imag( V_t ./ V )其中V_t是信号用窗函数导数加窗后的STFT。这个公式的直观解释是如果V的相位沿着时间快速变化说明这个时频格点的能量其实来自更高的频率而不是它当前所在的频率格。我在第一次跑通这个公式时心里还有点怀疑拿一个60Hz纯正弦波测了一下STFT谱图里能量集中在35Hz到85Hz一大片但用这个公式估计出来的频率点在每个时间点上都精确落在60Hz附近。那一刻我就明白相位信息比幅度信息“会说真话”。2.2 一维频率重分配与完整二维重分配有了瞬时频率估计重新分配方法的第一步就很自然了把V(t, ξ)处的一部分能量从原来的频率坐标ξ搬到瞬时频率omega_hat(t, ξ)对应的位置。这就是频率方向的一维重分配。完整版的重新分配方法还会做时间方向的重分配。它用另一个窗函数时间加权窗去估计局部群延迟然后把每个时频点的能量同时搬到新的时间和新的频率坐标。一维重分配把横向的模糊变清晰二维重分配同时把纵向和横向的模糊都变清晰。对于冲击成分较多、或者多个分量在时间上接近的信号二维版本效果更好。如果你用的是比较新的Matlab版本spectrogram函数自带reassigned选项可以直接得到重分配后的频率和时间估计省去自己实现群延迟部分的麻烦。不过自己写一遍一维版能帮助你建立直觉再切换到内置函数也不迟。2.3 SST和RM的取舍可逆性重新分配方法和同步挤压变换在第一步上是类似的都需要估计瞬时频率但搬移的对象不同。重新分配方法搬移的是能量也就是幅度平方。这意味着最终结果是一个非负矩阵形状好看但相位信息在搬移过程中丢掉了无法从重分配后的时频图反变换回原始信号。同步挤压变换搬移的是复数STFT系数。它把某个频率格上的复数值直接累加到瞬时频率对应的新格点上保留了相位因此挤压后的时频表示在原则上可以反变换回信号。这一点在需要做时频滤波、再重构信号的场景里非常重要。我整理过一个对比表方便你按需选型维度重新分配方法RM同步挤压变换SST搬移对象能量幅度平方复数STFT系数输出形式非负谱图复值时频矩阵可逆性不可逆可逆有条件时频锐化程度优秀优秀对噪声敏感度高高计算量中等中等如果只想要一张好看的时频图RM足够如果想做分量分离或者信号重构SST是更稳的选择。后面代码部分我以SST为主RM作为对照一起跑。3. 同步挤压变换的Matlab实现从公式到可运行代码3.1 SST算法流程总结SST的完整流程可以分四步对信号做STFT得到复数时频矩阵V同时用窗函数导数做一次STFT得到V_t。根据V和V_t估计瞬时频率omega_hat公式就是前文那个。设置一个幅度阈值把低于阈值的时频点标记为无效避免噪声把无效的瞬时频率估计搬来搬去。对每个有效时频点把复数系数V(t, ξ)累加到omega_hat(t, ξ)对应的频率格上。这一步就是“挤压”。挤压和重分配的区别就在第4步挤压是累加复数值重分配是搬运能量值。3.2 基于STFT的SST函数实现我写了一个基于STFT的SST函数没有依赖任何工具箱只要Matlab基础版本就能跑。参数都在函数头部注释里写清楚了。function [SST, omega_hat, fvec, tvec, V] stft_sst(x, fs, win_len, nfft, hop, gamma) % 基于短时傅里叶变换的同步挤压变换 % 输入: % x - 单通道信号行向量或列向量均可 % fs - 采样率, Hz % win_len - 窗函数长度点数 % nfft - FFT点数一般取大于win_len的2的幂 % hop - 帧移点数 % gamma - 幅度阈值系数相对峰值幅度的比例 % 输出: % SST - 同步挤压后的复值时频矩阵, 尺寸 时间 x 频率 % omega_hat- 瞬时频率估计矩阵, 单位 Hz % fvec - 频率坐标单边 % tvec - 时间坐标 % V - 原始STFT复数矩阵 % % 示例: % [SST, omega_hat, fvec, tvec, V] stft_sst(x, fs, 256, 1024, 8, 1e-4); if nargin 6 gamma 1e-4; end x x(:); % 转为行向量 N length(x); % 窗函数及其导数 win hann(win_len, periodic); win_deriv gradient(win); % 中心差分求导 % 帧数保证最后一帧不越界 n_frames floor((N - win_len) / hop) 1; tvec (0:n_frames-1) * hop / fs; % 只保留单边频率避免实信号负频率对挤压结果的污染 half floor(nfft / 2) 1; fvec (0:half-1) * fs / nfft; % 预分配 V zeros(n_frames, half); Vt zeros(n_frames, half); for m 1:n_frames idx (m-1)*hop 1 : (m-1)*hop win_len; xw x(idx) .* win; xwd x(idx) .* win_deriv; V(m, :) fft(xw, nfft); V(m, :) V(m, 1:half); % 截取单边 Vt(m, :) fft(xwd, nfft); Vt(m, :) Vt(m, 1:half); end % 瞬时频率估计fs/(2*pi) * imag(Vt / V) rel Vt ./ (V eps); omega_hat fs / (2*pi) * imag(rel); % 阈值判断 tol gamma * max(abs(V(:))); valid abs(V) tol; % 将瞬时频率映射到频率格索引 df fvec(2) - fvec(1); f_bins round((omega_hat - fvec(1)) / df) 1; f_bins min(max(f_bins, 1), half); % 边界截断 % 同步挤压把有效复数系数累加到目标频率格 SST zeros(n_frames, half); for m 1:n_frames idx_valid valid(m, :); if any(idx_valid) SST(m, f_bins(m, idx_valid)) SST(m, f_bins(m, idx_valid)) V(m, idx_valid); end end end这个函数是教学版没有做极致的性能优化但逻辑完整。我实际测试过60Hz纯正弦、线性chirp、多分量混合信号挤压后的轨迹都明显比谱图锐利。3.3 代码里的几个工程化细节有几点值得单独说一下都是我在调试中踩过的坑。第一截取单边频率很重要。实信号FFT后半部分是负频率的镜像如果不截掉负频率成分的瞬时频率估计出来是个负值挤压时会被边界截断逻辑拉到最低频格在图上形成一条低频伪迹干扰判断。第二阈值并不是可有可无的。没有阈值时所有微小的数值噪声都会参与挤压把噪声能量随机撒到整个时频面。有了阈值只有那些幅度可信的时频点才被搬移。gamma取相对峰值比例比取绝对值更稳健。第三窗函数导数直接用gradient计算即可。gradient默认按单位间隔做中心差分对应的是“每样本点”的变化再乘以fs后正好换算成“每秒”的导数。如果你用diff要小心长度减一的问题gradient保持了长度一致省事很多。4. 一次完整的对比实验谱图、RM、SST同台PK4.1 测试信号设计多分量混合为了全面看出差异我构造了一个三分量混合信号一个线性chirp、一个正弦调频分量、一个短时冲击。采样率1024Hz时长1秒。fs 1024; t 0:1/fs:1-1/fs; N length(t); % 分量1线性chirp 60Hz - 200Hz f0 60; f1 200; x1 sin(2*pi*(f0*t (f1-f0)*t.^2/2)); % 分量2中心频率100Hz附近的调频正弦 x2 0.7 * sin(2*pi*(100*t 5*sin(2*pi*3*t))); % 分量3位于0.5s附近的短时冲击 x3 1.0 * exp(-((t-0.5).^2)/0.0004) .* sin(2*pi*300*t); x x1 x2 x3;这个信号里有线性扫频、有窄带调频、有瞬态冲击三种分量对时频表示的要求各不相同正好用来检验三种方法的表现。4.2 三种时频图的可视化对比调用SST函数再用自写的频率重分配函数做RM对比。win_len 256; nfft 1024; hop 8; gamma 1e-4; [SST, omega_hat, fvec, tvec, V] stft_sst(x, fs, win_len, nfft, hop, gamma); % 自写一维频率重分配 RM zeros(size(V)); df fvec(2) - fvec(1); for m 1:size(V,1) for k 1:size(V,2) target round((omega_hat(m,k) - fvec(1))/df) 1; if target 1 target length(fvec) RM(m, target) RM(m, target) abs(V(m,k)).^2; end end end % 绘图 figure(Position, [100 100 1200 320]); subplot(1,3,1); imagesc(tvec, fvec, abs(V).^2); axis xy; colormap jet; colorbar; title(STFT Spectrogram); ylim([0 500]); subplot(1,3,2); imagesc(tvec, fvec, RM); axis xy; colormap jet; colorbar; title(Reassignment (freq-only)); ylim([0 500]); subplot(1,3,3); imagesc(tvec, fvec, abs(SST)); axis xy; colormap jet; colorbar; title(Synchrosqueezing Transform); ylim([0 500]);从图上能明显看到三点差异STFT谱图中chirp轨迹是一条带宽很大的亮带很难判断瞬时频率调频分量的正弦调制轨迹也被窗口抹平了冲击分量呈现纵向散布。RM结果中chirp轨迹变成了一条细线冲击的时间定位也更准确但同时能看到很多散落的噪声能量。SST结果和RM一样锐利但因为搬移的是复数值轨迹的连续性更好背景噪声相对不明显。4.3 瞬时频率估计精度的定量验证光看图还不够我习惯用量化指标说话。单独保留chirp分量用SST提取峰值轨迹和理论瞬时频率对比计算RMSE。% 单chirp信号 x_test sin(2*pi*(60*t (200-60)*t.^2/2)); [SST_test, ~, fvec_test, tvec_test] stft_sst(x_test, fs, 256, 1024, 8, 1e-4); [~, idx] max(abs(SST_test), [], 2); f_est fvec_test(idx); f_true 60 (200-60)