快速同步压缩变换：时频降采样与选择性重分配加速振动信号分析

做旋转机械故障诊断和结构振动监测的同行应该都有过这种体验：从齿轮箱或者轴承座上采回来一段多分量振动信号，时域波形看是密密麻麻的调制花样，频谱图上几十根谱线挤在一起，到底哪个是故障特征、哪个是转频谐波，光靠FFT根本分辨不清。于是我们习惯把信号搬到时频平面上看，同步压缩变换（SST）就是这类工具里效果最狠的一个，它能把模糊的时频带重新压回真实的瞬时频率脊线，多分量振动信号的分离和辨识一下子清楚很多。但SST有一个长期让人头疼的问题：计算代价高。尤其采样率一提上来、数据量一长，跑一次全网格SST要等到怀疑人生。我最近在MATLAB里实践了一套加速方案：基于时频降采样与选择性重分配做快速同步压缩变换，核心就两句话——先用降采样后的粗时频网格把有能量的区域圈出来，再只对圈内的时频点做精细的同步压缩重分配。这篇文章就把这套方案从头到尾拆开讲清楚。

1. 多分量振动信号时频分析，为什么需要“快速版”同步压缩变换

1.1 多分量振动信号带来的时频分析难题

多分量振动信号最常见的来源是齿轮箱、滚动轴承、往复机械和风力发电机传动链。齿轮啮合会产生啮合频率及其高次谐波，每一阶谐波周围还挂着一堆边频带，边频带的间距对应着故障轴的转频；轴承故障则在共振频带里激起周期性冲击，表现为若干条调频调幅明显的窄带分量。把这些成分合在一起，时域信号是非平稳、多模态、强耦合的，传统FFT只能给出整个时窗内的平均频率分布，完全无法回答“某个频率成分在哪个时刻出现、又沿时间怎么变化”这类问题。

时频分析方法里，短时傅里叶变换（STFT）实现简单、鲁棒性好，但受窗函数约束，时间分辨率与频率分辨率相互牵制，两个相邻分量靠得近了，时频图上的能量带就会糊成一团。Wigner-Ville分布的分辨率很高，却对多分量信号产生严重的交叉项干扰，图上会凭空多出许多“幽灵频率”。小波变换可以变尺度分析，但能量仍然沿着尺度轴扩散，脊线不如人意。同步压缩变换之所以被工程界追捧，本质上是把STFT扩散掉的能量通过相位信息“重新收拢”到瞬时频率曲线上，属于后处理重分配技术里理论和实际效果都比较成熟的一支。

不过成熟归成熟，SST在长数据上并不友好。你拿一段几分钟的连续监测信号去做全网格SST，先要算完整张STFT复数矩阵，接着对每一个时频格点做瞬时频率估计，最后还要逐点把能量搬到新的频率坐标上。这套流程把“分析高分辨率时频图”和“高计算量”绑定在一起，很多在线的故障诊断项目根本没耐心等它跑完，最后只能退回到普通STFT或谱峭度这类快速方法，这是很可惜的事。

1.2 SST计算量到底大在哪

先给一组直观数字。假设采样率fs=4096 Hz，连续监测10分钟，信号长度约245.8万点。用512点窗、32点帧移做STFT，单边频率点数是257，时间帧数约7.67万，总时频格点接近1972万个。标准SST对每个格点都要做这几件事：计算相位导数、做一次复数除法、取虚部换算瞬时频率、在频率轴上搜索最近目标格点、把系数累加到新的位置。哪怕MATLAB里用矩阵运算一次性算出全部瞬时频率矩阵，重分配阶段那1900万个点的“读-算-写”操作仍然会让内存带宽吃紧，实际耗时极其夸张。

如果为了提高频率分辨率把窗长加到1024，频率点变成513，时间帧数减半但总数反而更大。再加上很多研究者在论文里用带噪声的窄带信号测试，幅值很低的噪声区也要参与全网格重分配，而噪声区的相位是随机的，算出来的瞬时频率几乎毫无意义，不仅浪费计算还严重污染输出图。这就是问题的核心：标准SST把80%以上的计算浪费在了没有有效分量的时频区域上，而我们需要做的只是让它在真正有脊线的地方保持精度够用。

1.3 整体设计思路：先定位，再聚焦

快速SST的方案设计可以概括为“先粗后细、先选后搬”。第一阶段对信号做一个降采样的粗网格STFT，频率轴和时间轴都按倍数抽稀，在这个低分辨率时频图上用峰值搜索和能量检测识别出各个分量的候选频带；第二阶段回到原始细网格，只对候选频带内的时频点计算瞬时频率，并通过三条件判据“选择性”地执行能量重分配，其余格点直接忽略。这样做有三个好处：计算量按降采样倍数和时间抽稀倍数成倍下降，重分配操作只覆盖少数活跃区域；噪声区的无效瞬时频率估计被彻底排除，时频图反而比标准SST更干净；分量的瞬时频率脊线位置和幅值保持得不错，下游脊线提取、瞬时频率轨迹重构的精度几乎不受影响。

这套思路不改变STFT和同步压缩的核心数学过程，只是用“工程裁剪”的方式砍掉了大量冗余计算。后面我会按照原理、降采样实现、选择性重分配、完整代码和调试经验分别展开，代码全部基于MATLAB信号处理工具箱，大家拿去改参数就能跑。

2. 同步压缩变换原理与重分配核心逻辑

2.1 从STFT到同步压缩：能量为什么要重分配

同步压缩变换的原理解释起来并不复杂。对信号x(t)做加窗STFT得到时频系数G(t,eta)，其中t是时间中心，eta是频率索引。理想情况下，如果信号里只有一个慢慢变化的单频分量，STFT的系数在频率方向上不会只出现在那个瞬时频率点上，而是被窗函数主瓣“抹开”成一条有一定宽度的能量带。这正是时频不确定性的直接后果：你想看清频率在哪，就必须付出频率方向展宽的代价。

SST的精妙之处在于利用相位信息把这个展宽“收回去”。对STFT系数沿时间方向求偏导，可以得到每个时频格点处的相位变化率，它和真实瞬时频率之间存在一个可解析的换算关系。换算出来之后，原本分散在瞬时频率附近的能量带就可以被逐列重分配，全部集中到瞬时频率对应的单一频率点上。打个不太严谨的比方：STFT像是把一束激光打在毛玻璃上，光斑被散射成一片；同步压缩则是根据相位信息反推出每束光的原始方向，再把它们重新聚焦到各自的源位置。

对多分量振动信号来说，这个“聚焦”效果极其直观。两个调频分量在STFT图中可能互相重叠，能量带连成一片分不清边界；经过同步压缩后，每个分量被压缩成一条锐利的脊线，重叠区域也被拆分成两条独立轨迹。这就是为什么SST在机械故障诊断、模态参数识别和语音信号分析里都特别受用的原因——它极大地缓解了多分量时频图的可读性问题。

2.2 相位法瞬时频率估计的数值实现

实现SST绕不开的关键步骤是瞬时频率估计。最常用的方法是利用STFT系数的相位对时间的偏导数，公式可以写成：omega(t, eta) = eta - (1/(2*pi)) * imag(dG(t,eta)/dt / G(t,eta))。这里eta是当前栅格频率，dG/dt是STFT系数沿时间方向的导数，imag表示取虚部。物理含义是：在窗函数没有额外调频的情况下，STFT系数的时间演化速率直接携带了真实瞬时频率的信息。

在MATLAB里实现这段逻辑时，最容易踩的坑是时间方向差分的归一化。很多半路出家的代码省略了时间步长dt，导致估出来的瞬时频率整体偏大或偏小。正确做法是用中心差分除以后帧时间间隔：dGdt = (G(:,3:end) - G(:,1:end-2)) / (2dt)，其中dt = hop/fs，hop是帧移点数，fs是采样率。除以dt之后，dGdt的量纲才是1/s，取虚部除以2pi后得到的是Hz单位的频率修正量。第二个坑是除零保护，幅值接近零的格点在做复数除法时会爆炸，所以分母要加上eps。第三个坑是噪声区相位完全随机，在这些点上算出的瞬时频率毫无意义。这三个坑就是后面“选择性重分配”判据存在的根本原因——不是每个点都值得被重分配。

2.3 标准SST的MATLAB参考实现

在讲加速之前，先给一份教学级标准SST实现，方便大家对后续的优化有对比基准。这个版本不追求性能，只求逻辑清晰、可读性强：

function T_sst = standard_sst(x, fs, win_len, hop, amp_thresh) % x: 输入信号; fs: 采样率 % win_len: 分析窗长度; hop: 帧移; amp_thresh: 参与重分配的最低幅值阈值 n = length(x); nf = win_len/2 + 1; nt = floor((n - win_len)/hop) + 1; w = hann(win_len, 'periodic')'; dt = hop / fs; G = zeros(nf, nt); for k = 1:nt seg = x((k-1)*hop + (1:win_len)) .* w; spec = fft(seg, win_len); G(:, k) = spec(1:nf); end f_axis = (0:nf-1) * fs / win_len; % 时间方向中心差分估计瞬时频率 dGdt = zeros(nf, nt); dGdt(:, 2:end-1) = (G(:, 3:end) - G(:, 1:end-2)) / (2*dt); dGdt(:, 1) = dGdt(:, 2); dGdt(:, nt) = dGdt(:, nt-1); omega = f_axis(:) - imag(dGdt ./ (G + eps)) / (2*pi); T_sst = zeros(nf, nt); df = f_axis(2) - f_axis(1); for k = 1:nt for j = 1:nf if abs(G(j,k)) > amp_thresh % 幅值过低不参与搬移 [~, jj] = min(abs(f_axis - omega(j,k))); T_sst(jj, k) = T_sst(jj, k) + G(j,k); end end end end

注意这份代码的重分配部分用的是最原始的双层循环，纯粹为了演示原理。实际使用中，标准的SST重分配在MATLAB里往往要靠accumarray或者离散化bin索引来向量化，否则速度不堪入目。但仔细观察会发现，即便向量化了，只要时频矩阵规模一大，内存中间变量同样会爆炸。这也是我转向“时频降采样+选择性重分配”路线的主要原因。

3. 时频降采样：如何在粗网格上圈出有效能量区域

3.1 降采样的两个维度：频率轴与时间轴

时频降采样不是新概念，但多数文章只是把它作为预处理技巧一笔带过，几乎没有人把它的工程收益讲清楚。我把它拆成两个维度看：频率方向降采样和时间方向降采样。

频率方向降采样的本质是“降低FFT点数”。标准SST用win_len=512的FFT，单边频率点数是257；如果把参与定位的FFT长度直接减半到128，单边频率点数降到65，频率分辨率从原来的8 Hz变到32 Hz（以fs=4096为例）。有人一听降分辨率就摇头，但请注意：这一步的目标不是看清脊线细节，而是“圈出有能量的频带”，32 Hz的粗分辨率足以判断一个分量大致在哪个频段，根本不需要精细到零点几赫兹。时间方向降采样则更简单——相邻时间帧之间的STFT谱高度相关，调频分量在几十毫秒内不会有明显位移，所以可以先每隔一帧取一帧，比如帧移从64点放大到128点，时间帧数直接减半。

组合起来，粗网格STFT的规模大约是细网格的1/8（频率抽稀4倍、时间抽稀2倍）。这个规模下做峰值搜索和能量检测，MATLAB基本是瞬间出结果。粗网格的计算成本和后续节省的细网格重分配成本完全不成比例，两三百毫秒的粗定位开销能换走后续几秒甚至几十秒的重分配计算。

3.2 粗网格STFT实现与参数选择

粗网格的实现可以直接复用标准STFT的代码，只需要调整三个参数。我常用的参数组合如下表：

粗网格的FFT长度选128、窗长也是128，意味着窗函数只覆盖约31 ms信号（fs=4096时），频率分辨率比较粗糙但时间分辨率好，这对于捕捉冲击和调频的快速变化反而有利。帧移从细网格的64放大到128，时间帧数减少一半，定位计算量进一步压缩。需要注意，粗网格的窗长和FFT长度要一致，否则会出现频谱泄漏的怪异现象，定位结果就不准了。

从理论角度解释一下为什么粗定位能容忍这样的分辨率损失：STFT的时间-频率单元本质上是一个“可分辨面积”，粗细网格覆盖同一段信号，粗网格单元更大，但每个单元内如果确实存在一个调频分量，它的能量仍然会集中出现在对应频率附近的几个单元里。峰值搜索只要在这几个单元里找到一个局部极大值，就能把分量的中心频率带圈出来。后续细网格重计算会在原始分辨率上重新确定脊线的精确位置，所以粗网格导致的频率误差是可以通过第二阶段修正的。

3.3 活跃区域检测：能量掩码的构造方法

粗网格STFT算完之后，要把它转成一张“哪些时频位置有活跃能量”的掩码图。直接用全局阈值是新手最容易犯的错，因为振动信号中各分量能量差异可能非常大，强分量压过弱分量，全局阈值会直接把弱分量整条脊线吞掉。我的做法是逐时间帧做频率方向的局部峰值检测，以每个局部峰值中心向两侧扩展若干频率单元，形成候选频带掩码。

% 粗网格定位: 返回活跃频带掩码 function mask_band = coarse_localization(x, fs, nfft_c, hop_c) n = length(x); nf = nfft_c/2 + 1; nt = floor((n - nfft_c)/hop_c) + 1; w = hann(nfft_c, 'periodic')'; G = zeros(nf, nt); for k = 1:nt seg = x((k-1)*hop_c + (1:nfft_c)) .* w; G(:, k) = abs(fft(seg, nfft_c)); end mask_band = false(nf, nt); base_thresh = max(median(G(:)) * 3, 1e-6); % 基础阈值 for k = 1:nt [pks, locs] = findpeaks(G(:,k), ... 'MinPeakHeight', base_thresh, ... 'MinPeakDistance', 3); for i = 1:length(locs) lo = max(1, locs(i) - 3); hi = min(nf, locs(i) + 3); mask_band(lo:hi, k) = true; end end end

这里用findpeaks找局部峰，MinPeakDistance设为3，是为了避免把同一个展宽峰分裂成多个候选，MinPeakHeight用全图幅值中位数的3倍作为底线，既适应噪声水平变化，又不会定得太高以至于漏掉弱分量。每个峰左右各扩展3个粗频率单元，对应实际频率宽度约96 Hz（fs=4096、粗分辨率32 Hz），对绝大多数机械信号的脊线宽度来说是够的。如果你处理的信号分量特别多、频带密集，可以把扩展宽度和MinPeakDistance同时调小，搜索会细一些但漏检风险也随之上升。

4. 选择性重分配：三条件判据与能量搬移细节

4.1 为什么不能对所有点做重分配

标准SST之所以慢，根源在于它对时频平面上的“每一个”能量点都做了重分配。但实际操作起来你会发现，绝大多数点根本不应该参与搬移。以一段带噪声的轴承振动信号为例：信号里真正有物理意义的分量可能只占据时频平面的10%到20%，剩余部分要么是宽带噪声，要么是干扰冲击，要么是分量边缘的平滑过渡区域。对噪声点做瞬时频率估计，得到的是一个随机数，把它搬到随机的目标频率上去，除了在时频图上制造雪花一样的噪点，没有任何价值；对分量边缘的低幅值点做重分配，它们对应的相位估计同样不稳定，搬移结果会让脊线毛糙、边缘发虚。

所以我主张把“全网格重分配”改成“选择性重分配”，用一套显式的判据决定哪些点值得搬。这样做一方面减少了85%左右的重分配运算量，大幅提升速度；另一方面通过排除不可信点，让输出时频图的信噪比反而更高。这个思路特别适合工程场景：我们可以接受在细节上少一点理论优雅，但换来的是一个又快又干净的结果。

4.2 选择性判据的设计：幅值、掩码与频率偏差

我实际用的选择性判据一共三条，缺一不可。

第一条是幅度判据。参与重分配的时频点幅值必须显著高于噪声底数倍，推荐是5到20倍噪声标准差。幅值太低的点相位不可信，搬移只会制造伪脊线。

第二条是活跃掩码判据。该点必须落在粗定位生成的候选频带内。这条判据把重分配限制在粗网格检测到的脊线邻域内，直接从空间上切掉绝大多数非活跃区域。由于粗定位用的是降采样网格，掩码存在一定的位置模糊，但好处是绝不会因为粗网格分辨率不够而漏掉真实分量。

第三条是频率偏差判据。瞬时频率估计值omega与当前栅格频率eta的偏差必须在合理范围内，比如|omega - eta|小于3倍细网格频率间隔。这个判据的本质是判断相位估计是否“可信”——如果估计出来的瞬时频率离当前栅格太远，说明该点处于模态混叠区或者相位噪声区，强行搬运非但不能聚集能量，还会造成跨频带的伪峰。

三条判据组合成完整的选择逻辑后，重分配执行范围被压缩得非常干净。我在实际测试中曾用一段包含3个分量、信噪比约10 dB的仿真信号对比，全网格SST的时频图噪声斑点密集，选择性SST的图几乎只保留脊线本身，两侧背景干净到可以直接提取轨迹。

4.3 重分配坐标映射与能量保持细节

重分配操作的最终形态是把筛选后的STFT系数从原频率格点搬到瞬时频率对应的目标格点。实现时要注意保持系数本身是复数，而不是先取幅值再搬。原因很简单：如果后续要基于时频表示做信号重构或相位分析（比如提取某个分量的瞬时相位做阶比跟踪），必须保留复数系数里的相位信息。很多论文里只画时频幅值图，就容易写成搬|G|，等到要做重构时发现相位信息已经丢了，不得不从头再跑一遍。

边界处理上还需要留意瞬时频率估计值可能落在频率轴范围之外。接近直流分量或者超过奈奎斯特频率的ω，如果不加处理直接找最近格点，会让能量堆积在频谱两端形成虚假的亮线。稳妥的做法是遍历时遇到omega < f_axis(1)或omega > f_axis(end)就直接跳过，不参与搬运。别小看这个细节，很多程序跑出来的时频图在低频端有一条贯穿全图的亮带，十有八九就是这个原因。

5. 完整MATLAB实现与效果实测对比

5.1 总体流程与函数清单

快速SST的整体流程我拆成三步：粗网格定位、细网格STFT、选择性重分配。主脚本负责生成仿真信号和调用这三个环节，粗定位函数输出掩码，选择性重分配函数读取掩码和细网格STFT系数，输出稀疏的同步压缩时频矩阵。整个流程不依赖额外的第三方库，只要MATLAB装了信号处理工具箱就能跑。

流程设计的顺序是有讲究的。细网格STFT要等粗定位确定掩码之后再做，但细网格STFT本身可以一次性全频带算完，因为FFT用矩阵运算做并不慢。真正慢的是逐点相位估计和重分配搬运，所以细网格STFT全算完并不吃亏，反而可以用矩阵运算最大化吞吐。选择性重分配只在候选点上做循环，MATLAB的循环开销被控制在很小的规模内，整体性能趋于最优。

5.2 核心代码：主脚本与关键函数

下面给出一套完整可运行的实现。仿真信号是两个调频分量加高斯白噪声，参数设置方便大家复现：

% fast_sst_demo.m % 快速同步压缩变换: 时频降采样 + 选择性重分配 % 仿真信号: 两个调频分量 + 高斯噪声 clear; close all; clc; rng(42); fs = 2048; t = (0:8191)/fs; % 4秒信号 x = 1.0*sin(2*pi*(80*t + 10*sin(2*pi*0.4*t))) ... % 低频调频分量 + 0.7*sin(2*pi*(220*t + 15*sin(2*pi*0.25*t))) ... % 高频调频分量 + 0.15*randn(size(t)); % 噪声 noise_sigma = 0.15; win_len = 512; % 细网格FFT窗长 hop = 64; % 帧移 R_f = 4; % 频率降采样倍数 R_t = 2; % 时间降采样倍数 %% 1. 粗网格定位 nfft_c = win_len / R_f; hop_c = hop * R_t; mask_band = coarse_localization(x, fs, nfft_c, hop_c); %% 2. 细网格STFT [G_full, f_axis, t_axis] = my_stft(x, fs, win_len, hop); %% 3. 选择性重分配 T_sst = selective_sst(G_full, f_axis, t_axis, mask_band, ... noise_sigma, fs, R_f, R_t); %% 4. 画图对比 figure; subplot(2,1,1); imagesc(t_axis, f_axis, abs(G_full)); axis xy; ylabel('Frequency (Hz)'); xlabel('Time (s)'); title('STFT'); colorbar; subplot(2,1,2); imagesc(t_axis, f_axis, abs(T_sst)); axis xy; ylabel('Frequency (Hz)'); xlabel('Time (s)'); title('Fast SST (Downsampling + Selective Reassignment)'); colorbar;

% my_stft.m - 单边STFT实现 function [G, f_axis, t_axis] = my_stft(x, fs, win_len, hop) n = length(x); nf = win_len/2 + 1; nt = floor((n - win_len)/hop) + 1; w = hann(win_len, 'periodic')'; G = zeros(nf, nt); for k = 1:nt seg = x((k-1)*hop + (1:win_len)) .* w; G(:, k) = fft(seg, win_len); end G = G(1:nf, :); f_axis = (0:nf-1) * fs / win_len; t_axis = (0:nt-1) * hop / fs; end

% selective_sst.m - 选择性重分配 function T = selective_sst(G, f_axis, t_axis, mask_band, ... noise_sigma, fs, R_f, R_t) [nf, nt] = size(G); df = f_axis(2) - f_axis(1); dt = (t_axis(2) - t_axis(1)); % 将粗掩码块状放大到细网格尺寸 mask_full = kron(double(mask_band), ones(R_f, R_t)) > 0.5; mask_full = mask_full(1:nf, 1:nt); % 时间方向中心差分 dGdt = zeros(nf, nt); dGdt(:, 2:end-1) = (G(:, 3:end) - G(:, 1:end-2)) / (2*dt); dGdt(:, 1) = dGdt(:, 2); dGdt(:, nt) = dGdt(:, nt-1); % 瞬时频率估计 omega = f_axis(:) - imag(dGdt ./ (G + eps)) / (2*pi); T = zeros(nf, nt); amp = abs(G); for k = 1:nt for j = 1:nf if mask_full(j,k) && amp(j,k) > 5*noise_sigma ... && abs(omega(j,k) - f_axis(j)) < 3*df if omega(j,k) < f_axis(1) || omega(j,k) > f_axis(end) continue; end [~, jj] = min(abs(f_axis - omega(j,k))); T(jj, k) = T(jj, k) + G(j,k); end end end end

粗定位函数沿用前面3.3节里的coarse_localization。整套代码跑完后，上半张STFT图里两条分量是两条宽窄不一的能量带，下半张快速SST图里两条分量被压成两条锐利的细线，背景噪声明显变淡。如果频率偏差阈值设得合适，脊线的锐度跟标准SST几乎一致。

5.3 效果评价：速度与精度怎么权衡

我在自己的台式机上（i5-12400处理器、32 GB内存、MATLAB R2023b）用上面这份仿真信号做了对比。标准SST（前面2.3节的参考代码）耗时约4.2秒，其中一小半花在矩阵运算上，一大半花在双层循环的重分配阶段；快速SST总耗时约1.3秒，粗定位只占0.15秒，细网格STFT占0.4秒，选择性重分配占0.7秒，整体加速约3.2倍。如果把降采样倍数继续加大到R_f=8、R_t=4，理论上还能更快，但需要接受弱分量漏检风险的上升。

精度方面我没有发现快速SST有不可接受的退化。对两条调频分量做脊线峰值频率提取，快速SST与标准SST的最大偏差小于0.5 Hz，瞬时频率轨迹的重合度很高。尤其在高信噪比场景下，选择性重分配因为排除了噪声点干扰，时频图的能量集中度（用Rényi熵衡量）反而优于标准SST。低信噪比场景下，快速SST需要你把幅度判据的倍数调低一点，或者把粗定位的阈值放宽，这样虽然多引入一些计算点，但仍比全网格方案快得多。

6. 参数调试实战：常见问题与避坑心得

6.1 粗定位阶段漏掉弱分量

我最初测试三分量信号时，其中一个幅值很低的分量在结果图里完全消失了。排查半天发现是粗定位的MinPeakHeight定得太高，弱分量的谱峰没达到阈值。这个问题在振动信号里尤其常见：齿轮箱啮合频率的高次谐波能量一次比一次弱，最后一个可分辨的谐波分量可能比基频低了20 dB以上。解决思路是别用全局阈值一刀切，改用逐帧自适应阈值，比如取当前帧幅值中位数的若干倍作为该帧的峰值底线；同时把MinPeakDistance从3适当降到2，让密集频带里的弱峰也有出头机会。

如果是分量真实频率间隔本来就小于粗网格分辨率，漏检属于定位精度的硬限制，这时候要把R_f从4降回2，代价是粗定位计算量翻倍。我个人经验是，先跑一次R_f=4，如果发现弱分量漏检或相邻分量在粗定位图里糊成一团，再退到R_f=2最省事。

6.2 重分配后时频图出现条纹噪声

条纹噪声的典型表现是脊线附近出现细密的斜向亮纹，或者背景里有随机分布的亮点带。原因多半是频率偏差判据放行过宽：瞬时频率估计值跑偏了好几倍频率间隔，仍然被搬到了远处。解决办法是把第三个判据从3df收紧到1.5df或2df，跑完再观察脊线是不是变干净。如果脊线边缘出现断裂，说明阈值收得太紧，部分本应搬回脊线的点被丢弃了，适当回调到2df即可。这个参数几乎不需要数学推演，直接在测试信号上二分搜索几轮就能找到适合你信号的最优值。

另一种条纹来自掩码块状放大的边界阶梯效应。粗掩码的每个块在细网格上放大成矩形区域，矩形边界处瞬时频率估计被强行截断，容易产生细碎的搬移假象。可以用conv2对掩码做一个简单的平滑膨胀：

mask_full = conv2(double(mask_full), ones(7,7)/49, 'same') > 0.5;

这个操作把边界磨平，重分配时的连续性明显改善。代价是计算量小幅上升，但相比全网格重分配仍然是极小的开销。

6.3 参数敏感性速查表

我把整套方案涉及的核心参数整理成一张速查表，方便大家迁移到自己的数据上：

这套参数组合并非某种最优解，但它给出的是一个可以快速启动的起点。拿到新信号后，我会先跑一遍默认参数，看时频图哪里脏、哪里断，再针对性调两个参数，通常两三轮就能调到满意状态。

6.4 噪声标准差估计：一个容易被忽视的细节

选择性重分配的幅度判据写的是5*noise_sigma，但实际工程信号里的噪声标准差不是已知的。仿真里可以直接用噪声的生成标准差，真实数据却需要估计。我常用的估计方法是在粗定位阶段顺手统计：取粗网格STFT幅值矩阵中低于全局中位数的那部分幅值，计算其绝对中位差，再乘以1.4826换算成高斯噪声标准差。这个估计在纯信号区会被主瓣能量污染，所以最好拿低幅值部分算，既简单又稳，不用额外跑噪声估计工具。

还有一个相关的小细节：如果信号本身几乎没有噪声，比如来自仿真器或者经过强滤波的测试数据，幅度判据反而可能把真实分量也过滤掉，导致输出全是零。解决办法是给幅度判据加一个下限保护：当噪声估计值低于全局最大幅值的1%时，直接把倍数约束关掉，仅依赖掩码和频率偏差两个判据。这样纯信号下也能正常出结果，不会刹不住车。

6.5 从单段数据扩展到连续监测的工程建议

这套快速SST方法最终的价值要落到连续监测数据上。十分钟甚至几小时的振动数据如果一次性载入内存，细网格STFT矩阵本身就会吃光内存，顺序处理几乎不可行。我的建议是数据分块，每块处理2到5秒，块与块之间保留少量重叠（比如半个窗长），然后直接拼接时频矩阵。粗定位在每一块内独立进行，块间脊线的连续性靠重叠区自然衔接。如果机器有多核，用parfor把每一块的数据分派给不同工作线程，速度还能再翻倍。我在连续文件测试中把整段两小时的轴承数据切成了2400块，并行处理后总耗时从原先估计的数十分钟压缩到五分钟上下，已经接近在线监测的可接受范围。

这个扩展在实际项目里非常实用：粗定位天然具备自适应能力，每一块信号的候选频带按块内能量动态生成，即使整段信号里某个分量的频率漂移很大，也不会出现全局参数失配。我目前的项目已经把快速SST封装成函数，直接替换原来的标准SST调用，下游的脊线提取和瞬时频率轨迹重构接口完全不需要改动。对长期维护的诊断系统来说，这是最让人放心的改进方式。