上个月在复现一台减速箱的故障诊断算法时我又因为时频图不够锐利被折腾了一整天。标准同步压缩变换Synchrosqueezing Transform确实比短时傅里叶变换好看太多能量线细、分量清晰但它离真正的“实用”还差着一口气——噪声一大瞬时频率曲线就像喝多了酒相邻的调频分量还会互相“抢能量”。后来我把群延迟估计这个不算新的概念重新捡起来融合进同步压缩的频率重排环节才算把这个问题压下去。这篇文章就把这套基于MATLAB R2018A的改进方法完整拆开讲透包括算法原理、核心代码骨架、参数调优和我在复现过程中踩过的坑。如果你也在做机械故障诊断、语音分析、生物医学信号或任何非平稳信号的时频分析这篇内容应该能让你少走不少弯路。1. 标准同步压缩变换在实测信号里的三个“软肋”1.1 频率重排的原理回顾为什么SST比STFT锐利先把标准SST的底子过一遍。它的思路其实很简单信号经过短时傅里叶变换后时频图上的能量并不都落在真实瞬时频率上而是因为窗函数的频域泄漏散在真实频率周围的一小片区域里。SST做的事情就是把这个分散的能量“拿起来”按每个点估计出来的瞬时频率重新放回去。用数学点的话说对STFT结果S(t,ω)计算它相位关于时间的导数就得到局部瞬时频率估计f̂(t,ω)。然后对所有满足f̂某个频率f的系数做累加把能量集中到f上。这本质上是一种“频率方向的重新分配”所以它得到的时频图比STFT锐利得多也比小波尺度图更容易直接读数。我在R2018A里写这个功能时最直观的感受是它就像你先在一张纸上盖了很多模糊的印章然后拿一把尺子把每个印章上那些“跑偏”的墨迹按尺子刻度推回它本来应该在的格子里。尺子准不准决定了最后图面干不干净。1.2 软肋一噪声一上来相位导数就开始胡说标准SST的频率重排核心依赖相位对时间的导数。但是相位这个东西在低幅值区域非常不稳定我再打个比方一个信号的幅值一旦接近噪声底它的相位就基本被噪声主导了算出来的瞬时频率变成一堆随机数。这些随机数把能量搬到了乱七八糟的频率位置时频图上就会出现一堆“伪目标”看起来像下雪一样。我在仿真里给双分量信号加了5dB的高斯白噪声标准FSST做出来瞬时频率脊线在低频段勉强能看但到了高频段抖得厉害脊线周围布满散斑。对比真值一看频率估计误差比想象中大很多。这个问题的本质不是SST的框架错了而是它在重排前只用了一个“单方向”的相位信息信息量不够。1.3 软肋二交叉项与弱分量被强分量“吃掉”多分量信号是另一个重灾区。两个分量在时频平面上离得比较近时STFT系数在交叉区域是两个分量的叠加结果相位也是两个分量的“混合相位”。用混合相位去算瞬时频率重排位置大概率落在两个真实频率之间的某个空白地带形成伪脊线。更麻烦的是弱分量。当一个弱调频分量和一个强调频分量频率接近时弱分量的STFT系数幅值远小于强分量重排时能量很容易被判到强分量那边去。我一开始用标准SST处理一个“强线性调频弱正弦调频”的双分量信号时弱分量的时频脊在不少时段里直接消失了输出图里只能看到一条脊。这让我意识到单纯靠频率方向的一阶相位导数解决不了幅值对比悬殊的情况。2. 群延迟估计给同步压缩装上“第二只眼”2.1 群延迟的物理意义和很多人印象里不太一样群延迟这个概念在滤波器设计里常被说成“系统对不同频率成分的时间延迟”公式是τg(ω)−dφ(ω)/dω其中φ是系统相频特性。但在时频分析里如果把STFT的复系数看成同时包含幅度和相位的二维场我们其实可以问一个更基本的问题这个时频点的相位随着频率改变时到底在怎么变这个变化率就对应时频域里的群延迟。直观理解频率方向的相位变化率携带的是“信号能量沿时间轴应该从哪里聚集”的信息。瞬时频率告诉你“这根脊线的高度在哪里”群延迟告诉你“这根脊线的时间偏移有多大”。所以它们是互补的两条线索。标准SST只盯着瞬时频率相当于只看纵坐标移动从不理会横坐标的调整信息。2.2 为什么群延迟能补上频率重排的盲区我在测试中发现群延迟信息对噪声的敏感度和瞬时频率估计的敏感度不太一样。瞬时频率估计在高幅值时很准但低幅值区域乱跳群延迟估计在强噪声下虽然也会波动但它对“时间方向能量重心”的指示相对稳定特别是在瞬态成分、冲击脉冲这类信号上群延迟的定位能力比瞬时频率强很多。这让我想到一个改进方向把群延迟不当作最终参数而是当作频率重排的“校验信号”——用群延迟估计构造一个置信度权重把那些瞬时频率估计不可靠、而群延迟又明显异常的时频点给抑制掉。说得直白点相当于多装了一只眼睛两只眼睛看同一个点意见一致才采信。这个思路实现简单效果却很实在。2.3 改进算法的工作流程总览我这套改进方法的整体流程如下对输入信号做带窗STFT得到复系数矩阵S(t,ω)。在复系数矩阵上同时估计时间方向的相位导数得到瞬时频率场和频率方向的相位导数得到群延迟场。用瞬时频率场完成标准频率重排。用群延迟场构造逐点置信度掩码对重排结果做加权或门限处理。输出改进后的时频表示。流程不复杂但每一步的实现细节都能让人翻车。尤其是第二步直接决定整条链路的稳定性。我在最初版本里试图先用angle()取相位再做二维unwrap最后差分结果在含噪信号上完全不行。后来换成了复数域的一阶导数线性估计才把这个问题解决。具体实现见下一章。3. MATLAB R2018A下的核心实现与关键代码3.1 代码整体结构与数据组织先说项目里的代码组织。我习惯把这类工程拆成四个功能模块加一个演示入口结构如下improved_sst_ged/ ├── main_demo.m # 主脚本生成测试信号并跑对比 ├── functions/ │ ├── my_stft.m # STFT实现 │ ├── phase_derivative_est.m # 复数域相位导数估计 │ ├── improved_sst.m # 改进同步压缩主函数 │ └── plot_tfr_map.m # 时频图绘制 ├── data/ │ ├── demo_signal.mat # 双分量噪声仿真信号 │ └── case_bearing.mat # 实测轴承振动片段 └── references/ └── references.txt # 参考文献整理main_demo.m里我先把仿真信号载入跑三个版本的时频分析标准STFT、标准FSST、改进SST然后画在同一张对比图上。这样参数有没有调对一眼就能看出来。3.2 复数域相位导数估计核心中的核心相位导数的数值估计关键是不要直接对angle()结果做差分。正确做法是利用复数微分关系对于复函数S它的一阶导数与相位导数之间有d(arg S)/dx Im(dS/dx / S)等价地实际计算时可以写成d(arg S)/dx imag(conj(S).*dS_dx) ./ (abs(S).^2 eps)这里conj(S).*dS_dx的虚部除以|S|²就得到相位导数。这个形式避开了相位跳变问题数值上要稳得多。下面是我写的函数骨架function [dS_dt, dS_dw] phase_derivative_est(S, dt, dw) % PHASE_DERIVATIVE_EST 复数域相位导数估计 % S : STFT复数矩阵尺寸 [n_freq, n_time] % dt: 时间采样间隔秒 % dw: 频率采样间隔rad/s如果没有也可传1只在后续定标用 [n_freq, n_time] size(S); dS_dt zeros(n_freq, n_time); dS_dw zeros(n_freq, n_time); % 时间方向中心差分 if n_time 3 dS_dt(:, 2:end-1) (S(:, 3:end) - S(:, 1:end-2)) / (2 * dt); % 边界用单边差分 dS_dt(:, 1) (S(:, 2) - S(:, 1)) / dt; dS_dt(:, end) (S(:, end) - S(:, end-1)) / dt; end % 频率方向中心差分 if n_freq 3 dS_dw(2:end-1, :) (S(3:end, :) - S(1:end-2, :)) / 2; dS_dw(1, :) (S(2, :) - S(1, :)); dS_dw(end, :) (S(end, :) - S(end-1, :)); end % 如果dw不是1还需要统一量纲 dS_dw dS_dw * (1 / dw); end注意这个函数返回的是S的偏导数场不是最终的瞬时频率和群延迟。后续在improved_sst里再逐点计算IF imag(conj(S).*dS_dt) ./ (abs(S).^2 epsv); % 瞬时频率相关量 GD imag(conj(S).*dS_dw) ./ (abs(S).^2 epsv); % 群延迟相关量这里的epsv建议取1e-10甚至1e-12太小会除出尖峰太大又会把低幅值点的真实信息磨掉。我调试的时候用的是1e-10效果比较平衡。3.3 群延迟辅助重排的实现得到IF和GD之后改进的核心就在重排阶段。标准FSST只按IF把第(t,ω)点的能量搬到(t, ωIF·某个标定系数)的新频率槽里。我这边在两个地方做了改动。第一对每个时频点计算一个置信度权重weight 1 ./ (1 alpha * (GD - GD_median).^2)其中GD_median是同一时刻所有频率点上GD的中位数。这个权重的想法是当某点的群延迟异常偏离整体水平时说明该点很可能处于交叉项或噪声主导区频率重排的可信度低就把它压下去。第二频率重排时不再简单累加原始能量而是累加abs(S)的平方乘以weight。这样弱分量在正常工作区不会被压死强噪声区的伪目标会被明显抑制。核心代码大致是这样function TFA improved_sst(S, fvec, tvec, params) dt tvec(2) - tvec(1); dw 2 * pi * (fvec(2) - fvec(1)); % 根据频率轴定标 [dS_dt, dS_dw] phase_derivative_est(S, dt, 1); epsv params.epsv; IF imag(conj(S).*dS_dt) ./ (abs(S).^2 epsv); GD imag(conj(S).*dS_dw) ./ (abs(S).^2 epsv); % 群延迟置信度权重 alpha params.alpha; GD_med median(GD, 1); weight 1 ./ (1 alpha * (GD - GD_med).^2); % 频率重排 n_freq size(S, 1); TFA zeros(n_freq, size(S, 2)); for tIdx 1:size(S, 2) for fIdx 1:n_freq intIF round(IF(fIdx, tIdx) / (2*pi*(fvec(2)-fvec(1)))); newIdx fIdx intIF; if newIdx 1 newIdx n_freq TFA(newIdx, tIdx) TFA(newIdx, tIdx) ... abs(S(fIdx, tIdx)).^2 .* weight(fIdx, tIdx); end end end end当然两层for循环在信号较长时很慢。我在实际工程里做了矩阵化先把IF矩阵整体映射到目标索引再用accumarray做累加。下面这个写法更贴近我的实现intIF round(IF / (2*pi*(fvec(2)-fvec(1)))); idx_f repmat((1:n_freq), 1, size(S,2)) intIF; valid idx_f 1 idx_f n_freq; idx_f(~valid) 1; % 临时占位 lin_idx idx_f (0:size(S,2)-1) * n_freq; vals abs(S).^2 .* weight; TFA accumarray(lin_idx(:), vals(:), [n_freq*size(S,2), 1]); TFA reshape(TFA, n_freq, []);这样跑起来比for循环快十倍以上在R2018A下也能直接用。3.4 参数配置与运行Demo的路径参数配置我写在main_demo.m开头的结构体里params struct(); params.epsv 1e-10; params.alpha 1e-6; params.window_len 256; params.hop 64; params.nfft 1024; params.sigma 0.05; % 高斯窗标准差比例你第一次跑demo时建议先固定窗口长度逐步加大alpha从1e-8到1e-3观察噪声底的变化和弱分量是否丢失。alpha越大对群延迟异常点的压制越狠但压过头会把真实边缘也给抹掉。我的经验是alpha在1e-6附近对大多数双分量加噪信号都比较稳。4. 实验验证从仿真到噪声鲁棒性4.1 测试信号设计不要只做“干净”信号为了验证改进效果我构造了两个测试信号。第一个是双分量信号一个线性调频分量从100Hz扫到400Hz另一个正弦调频分量在250Hz附近往复波动幅值只有前者的三分之一。这专门用来考察交叉项和弱分量。第二个是模拟轴承早期故障的脉冲串信号脉冲重复频率约30Hz同时叠加两个谐波分量再加入高斯白噪声整体信噪比从10dB一直调到0dB。这里有个经验测试信号里一定要包含“弱分量强噪声时变频率”三者同时出现的工况。很多算法在漂亮的双分量无噪信号上表现接近完美一上实测就崩就是因为没覆盖到这些极端条件。4.2 评价指标Rényi熵和频率脊线误差算法好坏不能只靠“看着顺眼”判断。我用了两个量化指标。第一个是Rényi熵它衡量时频图的能量聚集程度。能量越集中在少数时频格点上熵越小散成一片熵就大。计算时用三阶Rényi熵R3 -1/2 * log2(sum(TFA.^3) / sum(TFA))第二个是频率脊线提取误差。从输出时频图里用峰值搜索提取两条脊线和真实瞬时频率比算均方根误差。这个指标直接说明算法在弱分量上的定位能力。4.3 结果解读改进SST到底赢在哪下面是这次测试里比较有代表性的结果信号具体参数不同数值会浮动但趋势稳定测试条件标准FSST Rényi熵改进SST Rényi熵弱分量脊线RMSE改善幅度双分量SNR10dB6.96.3约30%双分量SNR5dB7.66.6约45%双分量SNR0dB9.27.3约50%脉冲串谐波SNR5dB8.47.1约38%从时频图上看改进SST最明显的改善有两个一是噪声底明显变浅时频图干净不少二是弱分量那条脊线在原本消失的时段里重新露了出来。代价则是计算量增加了约25%到30%因为多算了一个频率方向的导数场和逐点权重。4.4 什么场景下改进不划算我也得说实话不是所有信号都适合这套改进。对单分量、高信噪比、平稳频率的信号标准SST已经很锐利群延迟置信度权重带来的收益很有限反而增加计算时间。另外如果信号本身包含大量持续时间极短的脉冲GD中位数估计会被脉冲能量拉偏此时alpha要调得很小或者改成按局部邻域计算中位数否则会误杀真脉冲。这个边界我在第一次测试时就踩过。5. 复现这套方法最容易翻车的五个地方5.1 相位导数千万别用angle()unwrap()diff()我最早试验时图省事对STFT系数直接取angle然后沿频率方向做一维unwrap再用diff差分。结果在噪声稍大时unwrap本身就会出错差的导数出现巨大尖峰时频图上就是几条竖直的亮线。换成复数域估计后问题彻底消失。这个坑我建议你直接绕开不要重复。5.2 窗长选择的两难没有万能配置窗长太长频率分辨率高但时间分辨率差信号瞬时频率变化快时脊线会被抹宽窗长太短时间定位好但频率方向上相邻分量容易糊在一起。群延迟估计对窗的频域形状也很敏感高斯窗在这个方向上的表现明显比矩形窗好。我的习惯是先根据感兴趣的频率范围估算扫频速率再大致让窗长覆盖扫频速率的倒数比如扫频速率200Hz/s时窗长取256样点采样率2kHz下约128ms起步然后看Rényi熵再微调。5.3 矩阵化能救命但小心内存上面那版accumarray写法快是快但临时矩阵idx_f和lin_idx都是和STFT矩阵一样大的double矩阵。如果STFT矩阵是1024×3000一个矩阵就是约24MB四个临时变量加起来逼近100MB。R2018A在一般办公电脑上还能顶住但如果你处理长信号或批量处理建议分块跑或者用single精度或者只在必要时保留lin_idx。5.4 群延迟的符号和量纲随STFT定义不同会变STFT有各种约定窗函数有没有归一化、时频平面用的是频率Hz还是角频率rad/s、窗函数是共轭在左还是右都会让GD和IF的整体符号差一个负号或者缩放因子。这个问题很难从教科书公式直接判断最实用的办法是拿一个单音信号做标定输入100Hz正弦看IF场是否稳定在100附近GD场是否在正确的常数附近。不对就整体取负或乘以系数比较省事。5.5 R2018A环境的具体限制R2018A里虽然已经有不少信号处理工具但内置的同步压缩函数对我来说更像“黑盒”拿不到中间相位导数场所以我全部自己实现只用到了矩阵运算、convn、accumarray这类元老级函数。这样代码的可移植性反而很好。唯一要注意的是某些后续版本里增强的函数如按维度的rmoutliers行为在R2018A里可能表现不完全一致我在代码里一律避免依赖新函数只用最基础的能力。如果你装了多个版本建议统一在R2018A上跑。6. 配套数据、文献整理与后续扩展方向6.1 数据与MAT文件的组织逻辑配套数据里我不光放了原始信号还存了元信息。比如demo_signal.mat里有五个变量sig信号序列fs采样率t时间轴f1_true分量1的真实瞬时频率f2_true分量2的真实瞬时频率。存真实频率的目的是让你复现时能直接算脊线RMSE不需要自己再估。case_bearing.mat里有振动信号和转速曲线但没有真实频率真值标定起来麻烦一点适合做定性观察。6.2 参考文献的整理思路我把参考文献按三个主题归档同步压缩变换经典文献、群延迟与时频重排相关文献、轴承故障诊断应用类文献。整理时重点记录每篇的“方法前提”比如哪篇用的是CWT域、哪篇是STFT域、窗函数用了什么形式这比记结论更重要。后续你要换到小波域或者加二阶项时回查速度会快得多。6.3 可以继续扩展的方向这套方法还有不少可玩空间。一是把一阶相位导数换成二阶做成高阶同步压缩适合扫频速率更快的信号二是让窗长随频率自适应变化用低频长窗高频短窗的思路兼顾分辨率三是把单通道扩展到多通道用阵列信号做群延迟估计的协同滤波。我自己已经试过高阶版本效果不错但代码还要再打磨以后有空再单独写一篇。最后再分享一个我自己的操作习惯每次调完改进算法我都会先跑一遍标准FSST作为baseline再跑改进版本把两幅图并排对比。这个简单的动作能避免很多“自我感觉良好”——毕竟时频分析这个领域视觉效果容易骗人量化指标又容易挑对自己有利的报告只有同条件下反复对比才知道改动是真有效还是只是巧合。你做实验时也建议保留这步对比流程。