前阵子处理一组变转速轴承振动数据的时候我用短时傅里叶变换STFT看了时频图转频附近的能量糊成一片边频带几乎看不见。当时我以为是窗长没调好折腾了好几组窗函数结果时间和频率分辨率还是跷跷板顾得了这头就顾不了那头。后来换成连续小波变换带宽是自适应的但高频段依然发虚单个冲击成分更是被拉成一条亮带。直到尝试了同步挤压变换Synchrosqueezing Transform, SST和重新分配方法Reassignment时频平面上的能量才真正“归位”到瞬时频率所在的位置整个图一下就清爽了。这个项目说到底就是围绕“海森堡不确定性原理HUP给时频分析带来的限制”展开的。传统线性时频变换比如短时傅里叶变换、小波变换都逃不开窗函数时宽与带宽约束窗选短了频率看不准窗选长了时间定位又模糊。同步挤压变换和重新分配方法属于一类“后处理”思路它们不做无损改良而是利用信号的局部相位信息在变换完成后把分散的能量重新分配到真实的瞬时频率位置让时频图更锐利。适合看非平稳信号、做故障诊断、分析生物医学信号或者处理地震勘探数据的人参考同时也非常适合想理解时频分析内部机制的研究生和工程师。1. 项目整体设计与思路拆解1.1 海森堡不确定性原理对时频分析的约束到底是什么很多人一听到海森堡不确定性原理第一反应是量子力学。实际上信号处理里的时频不确定性是更早被意识到的一件事任意一个有限能量信号它的时间范围和时间域的局部化程度没办法同时做到任意小。用均方根定义的时宽 (\Delta_t) 和频宽 (\Delta_f) 满足[ \Delta_t \cdot \Delta_f \ge \frac{1}{4\pi} ]这个不等式意味着窗口函数在时域越窄频率域就越宽反过来也一样。短时傅里叶变换就是在整个时间轴上滑动一个固定窗函数相当于每次用这个窗选出一段数据再计算这一段里的频率成分。窗函数选短一点时间定位会好一些但频率分辨率差窗函数选长一点频率分辨率上去了时间上又只能得到一段平均的结果。这个约束不是Matlab代码能绕过去的也不是加个更好窗函数就消失的它是由线性变换本身的结构决定的。时频图上那些看起来“模糊”的区域本质上就是窗函数和信号卷积后留下的涂抹痕迹。一个理想冲击信号经过STFT后会被窗函数拉成一条有一定宽度的竖直亮带一个理想单频连续波则会被拉成一条水平亮带。带内的能量分布就是窗函数的主瓣和旁瓣形状。我们看到的时频图其实是信号的真实时频分布和窗函数相互卷积的结果。既然卷积已经发生能不能在事后想办法把能量放回原来的位置这就是重新分配思想要解决的事情也是同步挤压变换的核心动机。1.2 重新分配方法把能量搬到重心位置重新分配方法最早可以追溯到Kodera等人在上世纪70年代的工作后来Auger和Flandrin把它系统化应用到了短时傅里叶变换和Wigner-Ville分布上。它的核心逻辑很朴素时频平面上某个点 ((t, f)) 处的频谱值并不是从这一点本身产生的而是来自附近某个“能量重心” ((t_0, f_0))。问题在于我们不知道这个真实重心在哪但可以通过局部相位对时间和频率的偏导数估计出重心偏移量然后把 ((t, f)) 处的能量移动到估计出来的 ((t_0, f_0)) 位置。这个概念有点像城市夜景的灯箱广告你在某个位置看见一片光晕它可能并不是光源本身而是广告灯被雾散射后的结果。重新分配要做的事情就是根据散射光的方向反推出灯真实在哪再把亮度搬回去。对STFT来说重新分配之后的时频图会有惊人的锐化效果正弦波的谱线几乎变成一条细线脉冲也会被压缩成一个点。我常用Matlab自带的spectrogram函数中的reassigned选项来做这一步只需要一行代码非常方便。但重新分配方法也有代价。第一它是一个非线性后处理不保留原始变换的相位结构所以很难直接重建信号第二它把很多点的能量搬走后时频平面会出现大量空点直接显示出稀疏的散点图像如果不做平滑会显得很乱第三它对噪声敏感因为相位导数在低能量区域被噪声主导容易出现虚假重定位。因此实际使用时要先适当降噪或者对估计的偏移量做平滑。1.3 同步挤压变换只朝频率方向“压缩”同步挤压变换是Daubechies等人在2011年前后系统化的一类方法它和重新分配思想同源但有一个关键的不同重新分配是在时间和频率两个方向同时移动能量而同步挤压只把能量沿频率方向压缩时间轴保持不变。这样做的好处是保留了一定的反变换能力可以从挤压后的时频表示中重建出原始信号或某个分量信号这在特征提取和信号分离场景里非常有用。具体思路是先用连续小波变换得到不同尺度下的系数 (W_x(a, b))其中 (a) 是尺度(b) 是时间平移。对于一个理想的单频分量CWT系数会在对应尺度附近形成一条具有一定宽度的“脊”。同步挤压通过计算系数沿时间方向的相位导数给出每个时间-尺度点上的“候选瞬时频率”[ \omega_x(a,b) -\frac{i}{2\pi} \frac{\partial_b W_x(a,b)}{W_x(a,b)} ]然后对每个时间点把所有尺度上候选频率落在某个小频带内的CWT系数叠加到该频带的中心频率处。这样原本宽的尺度脊就会被压缩成一条几乎看不见宽度的频率线。效果上同步挤压变换非常擅长把多分量信号里的密集频率成分分开比如轴承故障诊断中常常需要从变转速信号中分清转频、倍频和边频带。它和重新分配方法并不是互相替代的关系更像是两种不同取舍。如果只关心“看得清”重新分配方法简单粗暴能量聚焦效果明显如果关心“看清单个分量还能把它提取出来”同步挤压变换更有优势因为挤压过程保证了从时频图到信号的逆映射存在。1.4 为什么用Matlab做验证平台这个项目我选择Matlab来实现主要有三个原因。第一Matlab的信号处理工具箱里已经内置了连续小波变换以及不少时频分析函数比如cwt、spectrogram、wsst连同步挤压小波变换都有现成实现wsst新手可以直接调用而自己手写核心代码时也容易对照验证。第二Matlab的可视化能力很强画时频图、对比分辨率提升效果甚至做动画演示瞬时频率估计都比在其它语言里方便很多。第三学术和工程领域用Matlab做算法验证的比例很高贴出来的代码容易被接手复现。不过我也建议不要只停留在调用内置函数。自己手写一遍简化版同步挤压过程和重新分配过程能帮你真正理解相位导数和能量重定位的细节。这篇博文后面给出的示例就是偏教学性质的简化实现实际工程中可以直接用内置函数提升性能但理解原理一定得自己推导一遍。2. 核心细节解析与实操要点2.1 连续小波变换与同步挤压的数学关联连续小波变换的定义是[ W_x(a,b) \frac{1}{\sqrt{a}} \int x(t), \psi^*\left(\frac{t-b}{a}\right) dt ]其中 (\psi(t)) 是小波母函数。对纯谐波信号 (x(t) A\cos(2\pi f_0 t)) 来说CWT系数会在尺度轴接近 (a_0 f_c / f_0) 附近形成峰值其中 (f_c) 是小波中心频率。但小波本身有一定的带宽所以这个峰不是一条线而是有一定宽度的分布。同步挤压要做的就是根据 (W_x(a,b)) 的相位随时间的变化估计出该时间点的实际频率 (f_0)然后再把属于这个频率的所有小波系数在频率方向上累加。这里有一个容易忽略的细节相位导数估计的前提是 (W_x(a,b)) 不能太小否则相位无意义。所以实际代码里都会设一个阈值比如只处理 (\lvert W_x(a,b) \rvert \gamma) 的系数。否则噪声区域的随机相位会被计算成毫无规律的候选频率挤压之后得到的是杂乱无章的伪能量点时频图反而更难看。在实现中还要注意尺度向量的分布。同步挤压最终是把尺度轴映射到频率轴如果尺度只有稀疏的几个映射到频率网格后会出现空洞和齿状效应。我通常按“每倍频程多少个尺度”来控制常见取16、32、64。这个参数叫做voices per octave是决定频率方向精细程度的核心参数。2.2 同步挤压变换Matlab实现3个关键参数我把自己写同步挤压代码时的三个关键参数单独拿出来说因为调参往往比背公式更影响结果。第一个是母小波的中心频率。Matlab里如果使用cwtfilterbank并选用幅度调制复小波amor即Morlet小波默认的中心频率相当于6 rad/sample左右。中心频率越高小波在频率域的带宽越窄频率选择更精细但时间分辨率会下降。反之中心频率降低会让时间定位更准频率定位变差。遇到密集的谐波成分我一般会把中心频率调高一点比如用WaveletParameters, [8 1]代价是瞬态成分的起止时刻变得模糊。第二个是scale个数或者voices per octave。这个参数直接影响频率轴分辨率。如果设置为32每个倍频程内会有32个尺度频率轴可以画得比较平滑。取16时计算快但峰顶略粗糙取64时更精细但计算量也更大。对于轴承振动、语音信号这类频带较宽的数据32一般是一个不错的起点。第三个是挤压阈值gamma。它的作用是把无效小波系数过滤掉。我通常取CWT系数模的最大值的一个小比例比如0.01到0.05。阈值太高会把弱的分量丢掉阈值太低又会让噪声产生的随机频率混进结果。具体取值可以在一定范围内扫一下观察时频图的信噪比变化。2.3 重新分配方法的核心公式与代码逻辑重新分配方法在STFT框架下有比较直观的实现。假设STFT为[ S(t,f) \int x(\tau) w(\tau - t) e^{-i2\pi f \tau} d\tau ]窗函数 (w(t)) 的导数窗和 (t w(t)) 窗分别参与计算可以估计出局部群延迟和局部瞬时频率[ \hat{\tau}(t,f) t \mathrm{Re}\left(\frac{S_{tw}(t,f)}{S(t,f)}\right) ][ \hat{f}(t,f) f - \mathrm{Im}\left(\frac{S_{dw}(t,f)}{2\pi S(t,f)}\right) ]其中 (S_{dw}) 是用 (\frac{dw}{dt}) 作为窗函数得到的STFT(S_{tw}) 是用 (t w(t)) 作为窗函数得到的STFT。把 ((t,f)) 处的能量搬移到 ((\hat{\tau}, \hat{f}))就完成了重新分配。在Matlab里如果只是使用直接调用[ps, f, t, p, fc, tc] spectrogram(x, win, noverlap, nfft, fs, reassigned);返回值里fc是重定位后的频率tc是重定位后的时间。如果要做散点图可以用scatter(tc(:), fc(:), 1, ps(:))。如果要把重定位结果画成时频图而不是散点则可以用histogram2或者对tc/fc做二维累积。但需要注意spectrogram内置实现里窗函数导数的计算已经封装好了你不需要自己写出 (S_{dw}) 和 (S_{tw})。自己实现时最常出的问题是导数窗与原始窗的采样对齐以及用FFT时对信号补零的方式。我在早期实现重分配时没有对窗函数导数做同样的归一化结果重定位频率整体偏移了一个常数排查了很久才发现是导数窗忘了除以采样周期。3. 实操过程与核心环节实现3.1 构造一个能说明问题的测试信号为了直观展示同步挤压和重新分配带来的效果我构造了一个多分量测试信号采样率取1000 Hz时长2秒。第一个分量为正弦调频信号频率从80 Hz扫到180 Hz第二个分量为固定频率的谐波分量250 Hz第三个分量为一个短时冲击信号出现在0.6秒附近。这个组合能同时考验算法对调频轨迹、密集谐波和瞬态成分的分辨能力。fs 1000; t (0:1/fs:2-1/fs); x cos(2*pi*(80*t 30*t.*t/2)) 0.8*cos(2*pi*250*t); x(600) x(600) 5; % 冲击这里调频项如果用瞬时频率表示是 (80 30t) Hz所以1秒处频率约110 Hz2秒处约140 Hz。你可以先画一下STFT时频图会看到调频轨迹较粗250 Hz那条线始终是一条粗带子冲击在时频图上则显示为一片垂直亮带。3.2 从零开始写同步挤压变换代码我下面给一个教学用简化版本的同步挤压变换实现它基于Matlab的cwtfilterbank计算CWT系数然后做频率方向的挤压。注意这不是工业级完整实现但逻辑和真正的SST一致适合手写理解。function [Tfr, fgrid, tgrid] my_sst(x, fs, nvoice, gamma) % my_sst 简化版同步挤压小波变换 % 输入x信号fs采样率nvoice每倍频程尺度数gamma阈值系数 % 输出Tfr时频矩阵fgrid频率轴tgrid时间轴 x x(:).; N length(x); dt 1/fs; tgrid (0:N-1)*dt; % 构造小波滤波器组使用复Morlet小波 fb cwtfilterbank(SignalLength, N, ... SamplingFrequency, fs, ... VoicesPerOctave, nvoice, ... Wavelet, amor); [W, fa] fb.wt(x); % fa 是每个尺度对应的频率 nscales length(fa); % 时间方向中心差分估计相位导数 dW zeros(size(W)); dW(:, 2:end-1) (W(:, 3:end) - W(:, 1:end-1)) / (2*dt); % 候选瞬时频率 eps_val 1e-10; omega -imag(dW ./ (W eps_val)); % 定义频率网格 fmin min(fa); fmax max(fa); fgrid linspace(fmin, fmax, 512); df fgrid(2) - fgrid(1); Tfr zeros(length(fgrid), N); % 能量阈值 thr gamma * max(abs(W(:))); for b 1:N for a 1:nscales if abs(W(a,b)) thr fest omega(a,b); if fest fmin fest fmax idx round((fest - fmin)/df) 1; if idx 1 idx length(fgrid) Tfr(idx,b) Tfr(idx,b) W(a,b); end end end end end end这段代码里有几个细节要重点解释。第一dW中心差分会减小数据长度所以我保留了边界点的零值这也意味着边界附近的挤压效果会变差。第二候选瞬时频率计算里的imag结果跟“频率”单位有关这里得到的是Hz还是rad/s取决于CWT系数内部归一化方式我通常会在输出前做一个单位修正或用内置wsst对照检查。第三把尺度轴上的系数叠加到频率网格时我用的是幅度加权没有加 (a^{-3/2}) 这类尺度归一化因子严格来说会带来幅值偏差但对看时频图锐化效果影响不大。如果想直接使用Matlab内置的高质量实现下面的代码就够用了[sst, fwsst] wsst(x, fs, Waveletamor);wsst返回的是同步挤压小波变换结果和对应的频率轴。内置版本性能好边界处理也更好适合正式实验。3.3 重新分配方法的Matlab实现重新分配方法我一般直接用内置函数来跑它比手写稳定得多。核心代码很简单win hamming(256); noverlap 250; nfft 1024; [ps, f, t, p, fc, tc] spectrogram(x, win, noverlap, nfft, fs, reassigned); figure; scatter(tc(:), fc(:), 1.5, ps(:), filled); ylim([0 500]); xlabel(Time (s)); ylabel(Frequency (Hz));画出来的散点图能明显看到原来的连续亮带变成一系列非常集中的亮斑调频轨迹只剩一条很细的线。这个效果在能量重心散点图上非常震撼。需要提醒的是scatter里点的大小会同时影响视觉效果点太大又会糊成一团我一般取1到3之间的值。如果你像我一样希望把重分配结果转成类似时频矩阵的格式可以用二维累积edgesT linspace(0, 2, 512); edgesF linspace(0, 500, 512); T_reassign hist3([tc(:), fc(:)], {edgesF(2:end), edgesT(2:end)});这样得到的矩阵可以和STFT/CWT结果用同一套绘图代码展示。但要注意这样只保留了能量计数没保留幅度权重所以只能用于观察结构不能用于幅值对比。3.4 结果对比与分辨率提升的量化分析做完三种时频表示后我来对比一下效果。STFT用256点Hamming窗重叠250点FFT点数1024CWT使用复Morlet小波voices per octave取32SST用内置wsst重新分配用spectrogram的reassigned选项。从图形上能直观看到STFT的调频轨迹是一条有宽度的斜线边缘模糊CWT的轨迹在低频段较窄在高频段也变宽但总体上比STFT精细SST的轨迹几乎是一条极细的线250 Hz分量也锐化成了清晰的直线重新分配散点图则把冲击周围的能量集中成一个极小的区域调频轨迹同样非常清晰。如果要用数值指标来比较“分辨率集中度”我常用一个简单指标将时频矩阵的每个元素平方后求和再除以时频矩阵总能量的平方。这个值越大说明能量越集中function conc concentration(Tfr) P abs(Tfr).^2; conc sum(P(:).^2) / (sum(P(:)).^2); end在我构造的测试信号上STFT的集中度通常在0.02到0.05之间SST可以到0.3以上重分配散点图如果按矩阵累积算也能到0.2左右。这个差距非常明显说明能量确实被重新聚集了。4. 常见问题与排查技巧实录4.1 为什么挤压后的时频图出现竖直条纹我第一次用自写SST时时频图里总是出现很多竖直的亮条纹尤其在冲击附近。后来发现是相位导数在低幅度区域不稳定导致的。冲击信号会在某个时刻激发大量频率成分但每个成分的CWT系数幅度在边缘处迅速下降一旦低于阈值相位估计就被噪声主导候选项随机散布挤压后表现为竖直条纹。解决办法有几个方向。一是提高gamma阈值把不稳定的低幅度系数过滤掉二是对omega做中值滤波或平滑减少瞬时频率估计的抖动三是增加voices per octave让尺度网格更细避免映射时舍入误差放大。如果问题依然严重可以检查小波母函数中心频率是否过高太窄的频率带会让时间域的局部化变差冲击附近的干扰更明显。4.2 噪声环境下SST参数怎么调我拿一组信噪比大约6 dB的仿真信号试过直接跑默认SST时频图上会出现大量细碎噪声点弱分量几乎被淹没。此时最关键的不是继续调SST参数而是先做信号预处理。比如用带通滤波把感兴趣频带截出来或者先做一次短时傅里叶变换并降噪再对干净信号做SST。如果不想做预处理那可以尝试以下几组调整增加voices per octave到64让噪声点被分散到更多频率通道适当提高gamma让弱系数不参与挤压减小频率网格范围只保留目标频带这样高频噪声不会扩散到全图。还需要注意噪声环境下瞬时频率估计偏差会变大SST的脊可能出现断裂这时可以对挤压后的时频图再做一个低通平滑比如用imgaussfilt处理。4.3 边界效应和频率混叠的处理边界效应几乎是所有时频方法逃不掉的坑。小波变换在信号两端因为需要延拓边界处的系数不可靠之后挤压出来的边界频率也会失真。我一般有两种处理方式。第一在做SST之前对信号做对称延拓处理完后再裁剪到原始长度第二只关注时频图中间80%的区域边界部分不用于分析。在实际做故障诊断时我不会把边界处的结论当作判断依据。频率混叠问题主要出现在频率网格上限设置过高时。由于SST是在尺度域挤压如果最高频率超过奈奎斯特率的一半或者小波滤波器组的有效频带不含该频率就会出现虚假峰。解决方法是把fgrid上限限制在0.45倍的采样率以下同时用wsst默认返回的频率轴来检查。如果用的是自写代码最好根据小波滤波器组的中心频率范围来定义频率网格不要凭空设上限。4.4 计算速度太慢怎么办同步挤压的双循环在Matlab里效率很低尤其是我的my_sst函数里针对每个时间点循环所有尺度实际数据量一大就跑不动。我的经验有三条优化路线。第一用内置函数。wsst是经过优化的能比自写代码快一个数量级。第二矩阵化重写。挤压过程本质上是一个累加散点映射可以用accumarray代替双重循环代码虽然绕一点但速度提升非常明显。第三分频带处理。如果只关心某个频段可以先把数据降采样或者带通滤波再在小频率范围内做SST能省去大量冗余尺度计算。如果你需要处理长时间序列我建议改用分块策略把信号切成若干段每段做SST后再拼接但要注意每段的边界重叠区域要适当裁剪否则拼接处会出现接缝。4.5 内置wsst和自写SST结果不一致正常吗这个问题我被问过很多次。Matlab内置的wsst使用了更严格的尺度归一化和更精细的挤压逻辑同时它默认对边界做了延拓因此结果也会更平滑。自写版本如果只做教学演示幅度上与内置版本不一致是完全正常的关键要看瞬时频率轨迹是否对齐。如果你发现轨迹有偏移先检查中心频率的换算尤其要注意CWT滤波器组返回的频率轴和小波中心频率之间的关系。wsst返回的是Hz而自写代码中从相位导数得到的可能也是Hz但任何单位错误都会导致倍频偏差这个是最常见的坑。我个人在实际项目中的习惯是先用自己的简化代码理解算法行为等把参数逻辑理顺之后再切到内置版本做正式计算。这样既有对原理的掌控又能保证效率和稳定性。这个方向还能往多地方扩展。比如把同步挤压变换和变分模态分解结合先提取瞬时频率脊再重构各个分量用在旋转机械故障特征提取上效果很好也可以把重新分配方法的能量重心信息用于语音基频检测比单纯看STFT峰值稳定得多。时频分析的核心就是“在不确定性原理的约束下尽量恢复真实时频结构”而同步挤压变换和重新分配方法给了我们两个很实用的工具。真要说有什么心得那就是不要被“解决海森堡不确定性原理”这句话吓住它们并没有违背不确定性原理它们只是在信号已知的前提下用相位信息把被窗函数涂抹掉的能量重新放回原来的位置。理解这一点你后面调节参数、诊断问题时思路就会清晰很多。 同步挤压变换与重新分配：突破海森堡不确定性原理的时频分析
前阵子处理一组变转速轴承振动数据的时候我用短时傅里叶变换STFT看了时频图转频附近的能量糊成一片边频带几乎看不见。当时我以为是窗长没调好折腾了好几组窗函数结果时间和频率分辨率还是跷跷板顾得了这头就顾不了那头。后来换成连续小波变换带宽是自适应的但高频段依然发虚单个冲击成分更是被拉成一条亮带。直到尝试了同步挤压变换Synchrosqueezing Transform, SST和重新分配方法Reassignment时频平面上的能量才真正“归位”到瞬时频率所在的位置整个图一下就清爽了。这个项目说到底就是围绕“海森堡不确定性原理HUP给时频分析带来的限制”展开的。传统线性时频变换比如短时傅里叶变换、小波变换都逃不开窗函数时宽与带宽约束窗选短了频率看不准窗选长了时间定位又模糊。同步挤压变换和重新分配方法属于一类“后处理”思路它们不做无损改良而是利用信号的局部相位信息在变换完成后把分散的能量重新分配到真实的瞬时频率位置让时频图更锐利。适合看非平稳信号、做故障诊断、分析生物医学信号或者处理地震勘探数据的人参考同时也非常适合想理解时频分析内部机制的研究生和工程师。1. 项目整体设计与思路拆解1.1 海森堡不确定性原理对时频分析的约束到底是什么很多人一听到海森堡不确定性原理第一反应是量子力学。实际上信号处理里的时频不确定性是更早被意识到的一件事任意一个有限能量信号它的时间范围和时间域的局部化程度没办法同时做到任意小。用均方根定义的时宽 (\Delta_t) 和频宽 (\Delta_f) 满足[ \Delta_t \cdot \Delta_f \ge \frac{1}{4\pi} ]这个不等式意味着窗口函数在时域越窄频率域就越宽反过来也一样。短时傅里叶变换就是在整个时间轴上滑动一个固定窗函数相当于每次用这个窗选出一段数据再计算这一段里的频率成分。窗函数选短一点时间定位会好一些但频率分辨率差窗函数选长一点频率分辨率上去了时间上又只能得到一段平均的结果。这个约束不是Matlab代码能绕过去的也不是加个更好窗函数就消失的它是由线性变换本身的结构决定的。时频图上那些看起来“模糊”的区域本质上就是窗函数和信号卷积后留下的涂抹痕迹。一个理想冲击信号经过STFT后会被窗函数拉成一条有一定宽度的竖直亮带一个理想单频连续波则会被拉成一条水平亮带。带内的能量分布就是窗函数的主瓣和旁瓣形状。我们看到的时频图其实是信号的真实时频分布和窗函数相互卷积的结果。既然卷积已经发生能不能在事后想办法把能量放回原来的位置这就是重新分配思想要解决的事情也是同步挤压变换的核心动机。1.2 重新分配方法把能量搬到重心位置重新分配方法最早可以追溯到Kodera等人在上世纪70年代的工作后来Auger和Flandrin把它系统化应用到了短时傅里叶变换和Wigner-Ville分布上。它的核心逻辑很朴素时频平面上某个点 ((t, f)) 处的频谱值并不是从这一点本身产生的而是来自附近某个“能量重心” ((t_0, f_0))。问题在于我们不知道这个真实重心在哪但可以通过局部相位对时间和频率的偏导数估计出重心偏移量然后把 ((t, f)) 处的能量移动到估计出来的 ((t_0, f_0)) 位置。这个概念有点像城市夜景的灯箱广告你在某个位置看见一片光晕它可能并不是光源本身而是广告灯被雾散射后的结果。重新分配要做的事情就是根据散射光的方向反推出灯真实在哪再把亮度搬回去。对STFT来说重新分配之后的时频图会有惊人的锐化效果正弦波的谱线几乎变成一条细线脉冲也会被压缩成一个点。我常用Matlab自带的spectrogram函数中的reassigned选项来做这一步只需要一行代码非常方便。但重新分配方法也有代价。第一它是一个非线性后处理不保留原始变换的相位结构所以很难直接重建信号第二它把很多点的能量搬走后时频平面会出现大量空点直接显示出稀疏的散点图像如果不做平滑会显得很乱第三它对噪声敏感因为相位导数在低能量区域被噪声主导容易出现虚假重定位。因此实际使用时要先适当降噪或者对估计的偏移量做平滑。1.3 同步挤压变换只朝频率方向“压缩”同步挤压变换是Daubechies等人在2011年前后系统化的一类方法它和重新分配思想同源但有一个关键的不同重新分配是在时间和频率两个方向同时移动能量而同步挤压只把能量沿频率方向压缩时间轴保持不变。这样做的好处是保留了一定的反变换能力可以从挤压后的时频表示中重建出原始信号或某个分量信号这在特征提取和信号分离场景里非常有用。具体思路是先用连续小波变换得到不同尺度下的系数 (W_x(a, b))其中 (a) 是尺度(b) 是时间平移。对于一个理想的单频分量CWT系数会在对应尺度附近形成一条具有一定宽度的“脊”。同步挤压通过计算系数沿时间方向的相位导数给出每个时间-尺度点上的“候选瞬时频率”[ \omega_x(a,b) -\frac{i}{2\pi} \frac{\partial_b W_x(a,b)}{W_x(a,b)} ]然后对每个时间点把所有尺度上候选频率落在某个小频带内的CWT系数叠加到该频带的中心频率处。这样原本宽的尺度脊就会被压缩成一条几乎看不见宽度的频率线。效果上同步挤压变换非常擅长把多分量信号里的密集频率成分分开比如轴承故障诊断中常常需要从变转速信号中分清转频、倍频和边频带。它和重新分配方法并不是互相替代的关系更像是两种不同取舍。如果只关心“看得清”重新分配方法简单粗暴能量聚焦效果明显如果关心“看清单个分量还能把它提取出来”同步挤压变换更有优势因为挤压过程保证了从时频图到信号的逆映射存在。1.4 为什么用Matlab做验证平台这个项目我选择Matlab来实现主要有三个原因。第一Matlab的信号处理工具箱里已经内置了连续小波变换以及不少时频分析函数比如cwt、spectrogram、wsst连同步挤压小波变换都有现成实现wsst新手可以直接调用而自己手写核心代码时也容易对照验证。第二Matlab的可视化能力很强画时频图、对比分辨率提升效果甚至做动画演示瞬时频率估计都比在其它语言里方便很多。第三学术和工程领域用Matlab做算法验证的比例很高贴出来的代码容易被接手复现。不过我也建议不要只停留在调用内置函数。自己手写一遍简化版同步挤压过程和重新分配过程能帮你真正理解相位导数和能量重定位的细节。这篇博文后面给出的示例就是偏教学性质的简化实现实际工程中可以直接用内置函数提升性能但理解原理一定得自己推导一遍。2. 核心细节解析与实操要点2.1 连续小波变换与同步挤压的数学关联连续小波变换的定义是[ W_x(a,b) \frac{1}{\sqrt{a}} \int x(t), \psi^*\left(\frac{t-b}{a}\right) dt ]其中 (\psi(t)) 是小波母函数。对纯谐波信号 (x(t) A\cos(2\pi f_0 t)) 来说CWT系数会在尺度轴接近 (a_0 f_c / f_0) 附近形成峰值其中 (f_c) 是小波中心频率。但小波本身有一定的带宽所以这个峰不是一条线而是有一定宽度的分布。同步挤压要做的就是根据 (W_x(a,b)) 的相位随时间的变化估计出该时间点的实际频率 (f_0)然后再把属于这个频率的所有小波系数在频率方向上累加。这里有一个容易忽略的细节相位导数估计的前提是 (W_x(a,b)) 不能太小否则相位无意义。所以实际代码里都会设一个阈值比如只处理 (\lvert W_x(a,b) \rvert \gamma) 的系数。否则噪声区域的随机相位会被计算成毫无规律的候选频率挤压之后得到的是杂乱无章的伪能量点时频图反而更难看。在实现中还要注意尺度向量的分布。同步挤压最终是把尺度轴映射到频率轴如果尺度只有稀疏的几个映射到频率网格后会出现空洞和齿状效应。我通常按“每倍频程多少个尺度”来控制常见取16、32、64。这个参数叫做voices per octave是决定频率方向精细程度的核心参数。2.2 同步挤压变换Matlab实现3个关键参数我把自己写同步挤压代码时的三个关键参数单独拿出来说因为调参往往比背公式更影响结果。第一个是母小波的中心频率。Matlab里如果使用cwtfilterbank并选用幅度调制复小波amor即Morlet小波默认的中心频率相当于6 rad/sample左右。中心频率越高小波在频率域的带宽越窄频率选择更精细但时间分辨率会下降。反之中心频率降低会让时间定位更准频率定位变差。遇到密集的谐波成分我一般会把中心频率调高一点比如用WaveletParameters, [8 1]代价是瞬态成分的起止时刻变得模糊。第二个是scale个数或者voices per octave。这个参数直接影响频率轴分辨率。如果设置为32每个倍频程内会有32个尺度频率轴可以画得比较平滑。取16时计算快但峰顶略粗糙取64时更精细但计算量也更大。对于轴承振动、语音信号这类频带较宽的数据32一般是一个不错的起点。第三个是挤压阈值gamma。它的作用是把无效小波系数过滤掉。我通常取CWT系数模的最大值的一个小比例比如0.01到0.05。阈值太高会把弱的分量丢掉阈值太低又会让噪声产生的随机频率混进结果。具体取值可以在一定范围内扫一下观察时频图的信噪比变化。2.3 重新分配方法的核心公式与代码逻辑重新分配方法在STFT框架下有比较直观的实现。假设STFT为[ S(t,f) \int x(\tau) w(\tau - t) e^{-i2\pi f \tau} d\tau ]窗函数 (w(t)) 的导数窗和 (t w(t)) 窗分别参与计算可以估计出局部群延迟和局部瞬时频率[ \hat{\tau}(t,f) t \mathrm{Re}\left(\frac{S_{tw}(t,f)}{S(t,f)}\right) ][ \hat{f}(t,f) f - \mathrm{Im}\left(\frac{S_{dw}(t,f)}{2\pi S(t,f)}\right) ]其中 (S_{dw}) 是用 (\frac{dw}{dt}) 作为窗函数得到的STFT(S_{tw}) 是用 (t w(t)) 作为窗函数得到的STFT。把 ((t,f)) 处的能量搬移到 ((\hat{\tau}, \hat{f}))就完成了重新分配。在Matlab里如果只是使用直接调用[ps, f, t, p, fc, tc] spectrogram(x, win, noverlap, nfft, fs, reassigned);返回值里fc是重定位后的频率tc是重定位后的时间。如果要做散点图可以用scatter(tc(:), fc(:), 1, ps(:))。如果要把重定位结果画成时频图而不是散点则可以用histogram2或者对tc/fc做二维累积。但需要注意spectrogram内置实现里窗函数导数的计算已经封装好了你不需要自己写出 (S_{dw}) 和 (S_{tw})。自己实现时最常出的问题是导数窗与原始窗的采样对齐以及用FFT时对信号补零的方式。我在早期实现重分配时没有对窗函数导数做同样的归一化结果重定位频率整体偏移了一个常数排查了很久才发现是导数窗忘了除以采样周期。3. 实操过程与核心环节实现3.1 构造一个能说明问题的测试信号为了直观展示同步挤压和重新分配带来的效果我构造了一个多分量测试信号采样率取1000 Hz时长2秒。第一个分量为正弦调频信号频率从80 Hz扫到180 Hz第二个分量为固定频率的谐波分量250 Hz第三个分量为一个短时冲击信号出现在0.6秒附近。这个组合能同时考验算法对调频轨迹、密集谐波和瞬态成分的分辨能力。fs 1000; t (0:1/fs:2-1/fs); x cos(2*pi*(80*t 30*t.*t/2)) 0.8*cos(2*pi*250*t); x(600) x(600) 5; % 冲击这里调频项如果用瞬时频率表示是 (80 30t) Hz所以1秒处频率约110 Hz2秒处约140 Hz。你可以先画一下STFT时频图会看到调频轨迹较粗250 Hz那条线始终是一条粗带子冲击在时频图上则显示为一片垂直亮带。3.2 从零开始写同步挤压变换代码我下面给一个教学用简化版本的同步挤压变换实现它基于Matlab的cwtfilterbank计算CWT系数然后做频率方向的挤压。注意这不是工业级完整实现但逻辑和真正的SST一致适合手写理解。function [Tfr, fgrid, tgrid] my_sst(x, fs, nvoice, gamma) % my_sst 简化版同步挤压小波变换 % 输入x信号fs采样率nvoice每倍频程尺度数gamma阈值系数 % 输出Tfr时频矩阵fgrid频率轴tgrid时间轴 x x(:).; N length(x); dt 1/fs; tgrid (0:N-1)*dt; % 构造小波滤波器组使用复Morlet小波 fb cwtfilterbank(SignalLength, N, ... SamplingFrequency, fs, ... VoicesPerOctave, nvoice, ... Wavelet, amor); [W, fa] fb.wt(x); % fa 是每个尺度对应的频率 nscales length(fa); % 时间方向中心差分估计相位导数 dW zeros(size(W)); dW(:, 2:end-1) (W(:, 3:end) - W(:, 1:end-1)) / (2*dt); % 候选瞬时频率 eps_val 1e-10; omega -imag(dW ./ (W eps_val)); % 定义频率网格 fmin min(fa); fmax max(fa); fgrid linspace(fmin, fmax, 512); df fgrid(2) - fgrid(1); Tfr zeros(length(fgrid), N); % 能量阈值 thr gamma * max(abs(W(:))); for b 1:N for a 1:nscales if abs(W(a,b)) thr fest omega(a,b); if fest fmin fest fmax idx round((fest - fmin)/df) 1; if idx 1 idx length(fgrid) Tfr(idx,b) Tfr(idx,b) W(a,b); end end end end end end这段代码里有几个细节要重点解释。第一dW中心差分会减小数据长度所以我保留了边界点的零值这也意味着边界附近的挤压效果会变差。第二候选瞬时频率计算里的imag结果跟“频率”单位有关这里得到的是Hz还是rad/s取决于CWT系数内部归一化方式我通常会在输出前做一个单位修正或用内置wsst对照检查。第三把尺度轴上的系数叠加到频率网格时我用的是幅度加权没有加 (a^{-3/2}) 这类尺度归一化因子严格来说会带来幅值偏差但对看时频图锐化效果影响不大。如果想直接使用Matlab内置的高质量实现下面的代码就够用了[sst, fwsst] wsst(x, fs, Waveletamor);wsst返回的是同步挤压小波变换结果和对应的频率轴。内置版本性能好边界处理也更好适合正式实验。3.3 重新分配方法的Matlab实现重新分配方法我一般直接用内置函数来跑它比手写稳定得多。核心代码很简单win hamming(256); noverlap 250; nfft 1024; [ps, f, t, p, fc, tc] spectrogram(x, win, noverlap, nfft, fs, reassigned); figure; scatter(tc(:), fc(:), 1.5, ps(:), filled); ylim([0 500]); xlabel(Time (s)); ylabel(Frequency (Hz));画出来的散点图能明显看到原来的连续亮带变成一系列非常集中的亮斑调频轨迹只剩一条很细的线。这个效果在能量重心散点图上非常震撼。需要提醒的是scatter里点的大小会同时影响视觉效果点太大又会糊成一团我一般取1到3之间的值。如果你像我一样希望把重分配结果转成类似时频矩阵的格式可以用二维累积edgesT linspace(0, 2, 512); edgesF linspace(0, 500, 512); T_reassign hist3([tc(:), fc(:)], {edgesF(2:end), edgesT(2:end)});这样得到的矩阵可以和STFT/CWT结果用同一套绘图代码展示。但要注意这样只保留了能量计数没保留幅度权重所以只能用于观察结构不能用于幅值对比。3.4 结果对比与分辨率提升的量化分析做完三种时频表示后我来对比一下效果。STFT用256点Hamming窗重叠250点FFT点数1024CWT使用复Morlet小波voices per octave取32SST用内置wsst重新分配用spectrogram的reassigned选项。从图形上能直观看到STFT的调频轨迹是一条有宽度的斜线边缘模糊CWT的轨迹在低频段较窄在高频段也变宽但总体上比STFT精细SST的轨迹几乎是一条极细的线250 Hz分量也锐化成了清晰的直线重新分配散点图则把冲击周围的能量集中成一个极小的区域调频轨迹同样非常清晰。如果要用数值指标来比较“分辨率集中度”我常用一个简单指标将时频矩阵的每个元素平方后求和再除以时频矩阵总能量的平方。这个值越大说明能量越集中function conc concentration(Tfr) P abs(Tfr).^2; conc sum(P(:).^2) / (sum(P(:)).^2); end在我构造的测试信号上STFT的集中度通常在0.02到0.05之间SST可以到0.3以上重分配散点图如果按矩阵累积算也能到0.2左右。这个差距非常明显说明能量确实被重新聚集了。4. 常见问题与排查技巧实录4.1 为什么挤压后的时频图出现竖直条纹我第一次用自写SST时时频图里总是出现很多竖直的亮条纹尤其在冲击附近。后来发现是相位导数在低幅度区域不稳定导致的。冲击信号会在某个时刻激发大量频率成分但每个成分的CWT系数幅度在边缘处迅速下降一旦低于阈值相位估计就被噪声主导候选项随机散布挤压后表现为竖直条纹。解决办法有几个方向。一是提高gamma阈值把不稳定的低幅度系数过滤掉二是对omega做中值滤波或平滑减少瞬时频率估计的抖动三是增加voices per octave让尺度网格更细避免映射时舍入误差放大。如果问题依然严重可以检查小波母函数中心频率是否过高太窄的频率带会让时间域的局部化变差冲击附近的干扰更明显。4.2 噪声环境下SST参数怎么调我拿一组信噪比大约6 dB的仿真信号试过直接跑默认SST时频图上会出现大量细碎噪声点弱分量几乎被淹没。此时最关键的不是继续调SST参数而是先做信号预处理。比如用带通滤波把感兴趣频带截出来或者先做一次短时傅里叶变换并降噪再对干净信号做SST。如果不想做预处理那可以尝试以下几组调整增加voices per octave到64让噪声点被分散到更多频率通道适当提高gamma让弱系数不参与挤压减小频率网格范围只保留目标频带这样高频噪声不会扩散到全图。还需要注意噪声环境下瞬时频率估计偏差会变大SST的脊可能出现断裂这时可以对挤压后的时频图再做一个低通平滑比如用imgaussfilt处理。4.3 边界效应和频率混叠的处理边界效应几乎是所有时频方法逃不掉的坑。小波变换在信号两端因为需要延拓边界处的系数不可靠之后挤压出来的边界频率也会失真。我一般有两种处理方式。第一在做SST之前对信号做对称延拓处理完后再裁剪到原始长度第二只关注时频图中间80%的区域边界部分不用于分析。在实际做故障诊断时我不会把边界处的结论当作判断依据。频率混叠问题主要出现在频率网格上限设置过高时。由于SST是在尺度域挤压如果最高频率超过奈奎斯特率的一半或者小波滤波器组的有效频带不含该频率就会出现虚假峰。解决方法是把fgrid上限限制在0.45倍的采样率以下同时用wsst默认返回的频率轴来检查。如果用的是自写代码最好根据小波滤波器组的中心频率范围来定义频率网格不要凭空设上限。4.4 计算速度太慢怎么办同步挤压的双循环在Matlab里效率很低尤其是我的my_sst函数里针对每个时间点循环所有尺度实际数据量一大就跑不动。我的经验有三条优化路线。第一用内置函数。wsst是经过优化的能比自写代码快一个数量级。第二矩阵化重写。挤压过程本质上是一个累加散点映射可以用accumarray代替双重循环代码虽然绕一点但速度提升非常明显。第三分频带处理。如果只关心某个频段可以先把数据降采样或者带通滤波再在小频率范围内做SST能省去大量冗余尺度计算。如果你需要处理长时间序列我建议改用分块策略把信号切成若干段每段做SST后再拼接但要注意每段的边界重叠区域要适当裁剪否则拼接处会出现接缝。4.5 内置wsst和自写SST结果不一致正常吗这个问题我被问过很多次。Matlab内置的wsst使用了更严格的尺度归一化和更精细的挤压逻辑同时它默认对边界做了延拓因此结果也会更平滑。自写版本如果只做教学演示幅度上与内置版本不一致是完全正常的关键要看瞬时频率轨迹是否对齐。如果你发现轨迹有偏移先检查中心频率的换算尤其要注意CWT滤波器组返回的频率轴和小波中心频率之间的关系。wsst返回的是Hz而自写代码中从相位导数得到的可能也是Hz但任何单位错误都会导致倍频偏差这个是最常见的坑。我个人在实际项目中的习惯是先用自己的简化代码理解算法行为等把参数逻辑理顺之后再切到内置版本做正式计算。这样既有对原理的掌控又能保证效率和稳定性。这个方向还能往多地方扩展。比如把同步挤压变换和变分模态分解结合先提取瞬时频率脊再重构各个分量用在旋转机械故障特征提取上效果很好也可以把重新分配方法的能量重心信息用于语音基频检测比单纯看STFT峰值稳定得多。时频分析的核心就是“在不确定性原理的约束下尽量恢复真实时频结构”而同步挤压变换和重新分配方法给了我们两个很实用的工具。真要说有什么心得那就是不要被“解决海森堡不确定性原理”这句话吓住它们并没有违背不确定性原理它们只是在信号已知的前提下用相位信息把被窗函数涂抹掉的能量重新放回原来的位置。理解这一点你后面调节参数、诊断问题时思路就会清晰很多。