企业官网建设流程全解析

前几天一个做故障诊断的同行问我，能不能把同步压缩变换的时频图做得再聚焦一点。他的信号里有两个频率靠得非常近的分量，短时傅里叶变换糊成一片，标准同步压缩变换勉强能分开，但边缘和幅值一直不满意。我给他推荐了基于群延迟估计做重分配的路线，顺手在 MATLAB R2018A 上把改进群延迟估计同步压缩变换，也就是我下面统称的 IGD-SST，整个跑通了一遍，数据、代码、参考一起整理出来。这篇文章就是这个过程的完整记录。

这套方法不是什么高不可攀的花架子，它就是针对“时频分析能量不集中”这个老问题做的一次重分配规则升级。适合正在做机械故障诊断、生物医学信号处理、地球物理或者振动分析的工程师和研究生，尤其是信号里有多分量、强调频、低信噪比情况的同学。如果你只想快速跑出来一张漂亮的时频图，也可以直接把我下面的测试信号换成你自己的数据，核心代码不用动。

1. 先理解同步压缩：它到底在压缩什么

1.1 从短时傅里叶变换说起

所有时频分析问题的起点，几乎都是短时傅里叶变换。打个比方，一段信号就像一卷长长的电影胶片，瞬时频率是胶片上每一秒的真实剧情。短时傅里叶变换的干法是，拿一个开了一段小窗口的放大镜，沿着时间轴一段一段地看胶片，把每一段做傅里叶变换，得到“哪个时刻大概有哪些频率成分”。问题在于，放大镜一旦选定了宽度，频率分辨率和时间分辨率就打架了：窗口拉长，频率看得清但时间模糊；窗口缩短，时间定位准了但频率含糊。

这个矛盾不是参数调一调就能完全绕开的，它是短时傅里叶变换本身的数学约束。实际信号里最常见的麻烦，是两个频率差十几赫兹的分量，窗口稍微短一点就叠在一起，肉眼无法分辨；窗口长一点，它们中间的瞬态变化又会被抹平。同步压缩变换的初衷，就是在这种分辨率限制下，再往前推一把，把已经“糊”在时频图上的能量，按某种准则重新塞回它该去的位置。

1.2 同步压缩的核心动作：重分配

同步压缩变换第一次被大家广泛认识，靠的是 Daubechies 等人 2011 年那篇把小波变换结果做重分配的经典工作。它做的事情用一句话说：原来时频图上每个点的能量，并不是自己老实待在原地的，很多能量是从相邻区域“漏”过来的，那么只要我能估计出每个点对应的真实瞬时频率，就能把这些漏掉能量全收回到瞬时频率曲线上。

这个“估计真实瞬时频率”的动作，在同步压缩里叫频率重分配算子。标准的做法是取短时傅里叶变换的相位信息，对时间求偏导，算出每一个时频点的瞬时频率偏移量。整个算法最后输出一张经过同步压缩的时频图，看起来比原始 STFT 图锐利得多，线状分量变成了一条条细亮的脊线。很多用过同步压缩工具包的人，第一感受就是“图像像是被锐化过了”，其实背后是实实在在的重分配，不是后处理滤镜。

1.3 为什么只沿频率轴压不够

我刚开始用标准同步压缩变换时，也觉得挺惊艳，直到碰上一个线性调频信号才意识到问题。标准同步压缩只沿频率方向重分配，它假设在窗函数覆盖的很短时间内，信号的瞬时频率几乎不变，可以当成恒定频率来处理。可是遇到调频斜率很陡的信号，比如频率从 100 Hz 一路扫到 400 Hz，这个假设就不太成立了。

同样的道理，群延迟信息也被忽略了。群延迟可以粗浅理解成“某个频率成分的真实到达时刻”。标准同步压缩只在频率方向做文章，相当于只把同一时刻附近漏出去的能量按频率收拢，却没有对时间方向的偏移做纠正。很多瞬时频率变化剧烈的信号，压缩之后虽然频率方向变细了，时间方向的拖尾依然存在，这就是为什么有人觉得同步压缩后的时频图仍然“不够干净”。

1.4 改进路线到底改在哪

我做这套 IGD-SST 的时候，核心思路是在原来频率重分配的基础上，把群延迟估计也拉进来，形成一个联合重分配：每个时频点不再只沿频率轴移动，而是同时沿频率轴和时间轴移动到自己最该待的位置。听起来简单，真正搞的时候还是踩了几个坑，比如群延迟估计很容易被相位卷绕带偏，阈值设得太低时噪声点会被无差别重分配，整个时频图全是毛刺。

所以我在实现里做了两个收紧动作。第一，用复谱比值来估计偏导数，避开直接对相位做差分的那些麻烦；第二，加了一个能量门槛，幅值低于门槛的点不参与重分配。这两点改进让整张时频图干净了很多，也正是我把这个方法叫做“改进”的主要原因。

2. 算法原理：群延迟怎么参与重分配

2.1 两个偏导谱的计算

既然要同时做频率方向和时间方向的重分配，就得先拿到两个关键量：短时傅里叶变换对时间的偏导谱，以及对频率的偏导谱。很多新手在这里容易蒙，觉得偏导谱是不是还要去摆弄相位解卷绕。其实不用，在 MATLAB 里可以直接绕开相位，用窗函数和一阶导数窗各做一次傅里叶变换就行。

假设我们有一帧加窗后的信号片段seg，窗函数是win，那么当前这一帧的短时傅里叶谱就是：

S = fft(seg .* win, nfft);

要估计对时间的偏导谱，就把窗换成它的一阶导数winD，再做一次 FFT：

St = fft(seg .* winD, nfft);

要估计对频率的偏导谱，就用时间加权的窗win .* tp再做一次 FFT，其中tp是以窗中心为零点的时间坐标：

Sf = fft(seg .* win .* tp, nfft);

这三条代码几乎就是整个改进算法的地基。用这种方式计算偏导谱，好处是不需要关心相位是怎样从 -pi 跳到 pi 的，复数的比值天然把相位卷绕问题吞掉了，跑起来特别稳。

2.2 瞬时频率与群延迟联合重分配准则

有了S、St、Sf三组复数谱，接下来就可以算重分配量。标准同步压缩里估计瞬时频率，我直接写成：

[ \hat{\omega}(t,\eta) = \eta + \frac{1}{2\pi} \operatorname{Im}\left( \frac{S_t(t,\eta)}{S(t,\eta)} \right) ]

这个式子很好理解：当前频率格点是eta，加上相位随时间变化带来的修正量，就得到该点能量真正的瞬时频率。群延迟估计则是另一个方向：

[ \hat{\tau}(t,\eta) = t - \frac{1}{2\pi} \operatorname{Im}\left( \frac{S_f(t,\eta)}{S(t,\eta)} \right) ]

它告诉我们这个频率成分相对当前窗中心的到达时间偏移。两个式子在代码里就是两行：

omega = eta + imag(St(k) ./ S(k)) / (2*pi); tau = t_n - imag(Sf(k) ./ S(k)) / (2*pi);

算出了omega和tau，这个时频点的能量就可以从原来的(t_n, eta)位置，整体挪到(tau, omega)位置。这里我把瞬时频率修正和时间方向修正合在一起，也就是联合重分配名称的由来。

2.3 完整流程串一遍

整个 IGD-SST 算法的执行顺序，我平时喜欢浓缩成五个步骤。第一步，对输入信号逐帧加窗做短时傅里叶变换，得到复数谱 S。第二步，用窗导数谱 St 和时间加权谱 Sf，分别得到两个偏导谱。第三步，对每个时频点按上面的公式计算瞬时频率和群延迟。第四步，判断幅值是否超过能量门槛，超过才执行重分配。第五步，把能量叠加到目标时频格点上，输出压缩后的时频谱。

之所以把能量门槛放在第四步而不是一开始就过滤，是因为计算门槛需要先看整个 S 的动态范围。我习惯取gamma = 0.01 * max(abs(S(:)))，信噪比特别低或者噪声比较重的时候，会把门槛调到0.03左右。这个值不是越大约好，门太高会把弱分量也丢光，后面会专门讲参数怎么调。

2.4 和其他重分配方法的对照

把 IGD-SST 和几个常见方法放在一起看，区别会更清楚。标准同步压缩变换 SST 只压缩频率方向，速度快但时间方向拖尾没救；重分配方法二分法 reassignment 最早由 Auger 和 Flandrin 提出，两个方向都压，但它没有后面同步压缩那套滤波重构机制；群延迟同步压缩变换 GD-SST 则把重心放在时间方向，对窄带瞬态信号效果好，但在强调频信号上频率方向依然不够集中。

我做 IGD-SST 的思路是两边都管，同时利用偏导谱做联合估计，再加门槛抑制噪声。和经典重分配方法相比，它保留了同步压缩变换的可解释性，也保留了后续信号重构的可能性。下面这张表是我自己常用的对比思路：

3. 在 MATLAB R2018A 上动手实现

3.1 环境准备

我在这个项目里用的是 MATLAB R2018A，运行环境是 Windows 10 的机器。需要确认 Signal Processing Toolbox 已经装好，因为后面要用的gausswin、kaiser这类窗函数都在这个工具箱里。如果只有 MATLAB 基础环境，可以考虑手动建一个高斯窗向量，但没必要和工具箱过不去，还是装上最省事。

另外要注意一个细节：R2018A 里还没有后来新版 MATLAB 那种封装得比较完整的stft命令，至少我自己在这个版本里更喜欢直接用fft自己写短时傅里叶变换。这样看起来麻烦，实际上有一个很大的好处，就是我可以随时拿到复数谱的中间结果，后面的偏导谱、重分配全都基于这些复数谱来做，不需要额外去猜工具箱内部的实现。

脚本存放路径上，我建议不要放到中文路径或者带空格的目录里，MATLAB 对这类路径偶尔会出一些莫名其妙的问题。我自己一开始放在“桌面\新文件夹”下面，结果脚本运行半路报错找不到函数，把整个文件夹复制到D:\IGDSST之后就一切正常了。

3.2 测试信号：数据怎么造

这个项目的“含数据”部分，我准备了一个带两个分量的测试信号。第一个分量是线性调频信号，频率从 150 Hz 扫到 250 Hz；第二个分量是正弦调频信号，中心频率 350 Hz，调制幅度 30 Hz，调制频率 3 Hz。采样率 1024 Hz，时长 1 秒。为了让测试更接近真实场景，我加了很小的白噪声，信噪比大约 20 dB。

这样的数据构造能同时考验两件事：线性调频分量检验群延迟联合重分配对强调频信号的处理能力，正弦调频分量检验算法能不能追踪快速变化的瞬时频率曲线。直接把下面的代码跑一遍，就能生成x这个行向量，也方便后面替换自己的数据。

fs = 1024; t = (0:1023) / fs; f0 = 150; f1 = 250; x_chirp = chirp(t, f0, t(end), f1); x_fm = sin(2*pi*350*t + 2*pi*30/(2*pi*3) * sin(2*pi*3*t)); x = x_chirp + x_fm + 0.03 * randn(size(t));

如果你用的是自己的工程数据，直接把x定义成你自己的信号就行，但要保证它是行向量，采样率fs也要改对。我之前有一次信号是列向量，代码里索引全乱套，输出时频图直接歪了，所以要注意这一点。

3.3 主函数完整代码

我把核心函数命名为igdsst，输入信号、采样率、窗长、FFT 点数和能量门槛都作为参数传进去。首先对信号加高斯窗，高斯窗的参数可以根据主瓣宽度需求调。然后逐帧做 FFT，在循环里同时算出S、St、Sf，再进行联合重分配。为了控制篇幅，我没有做矩阵向量化的极致优化，但逻辑非常直观。

function [TFR, freq, tvec] = igdsst(x, fs, winLen, nfft, gamma) % IGD-SST: Improved Group Delay Synchrosqueezing Transform % x : 一维信号行向量 % fs : 采样率 % winLen: 窗长度，建议奇数 % nfft : FFT 点数 % gamma: 能量门槛 x = x(:).'; % 强制转为行向量 N = length(x); Lh = floor(winLen / 2); win = gausswin(winLen, 8).'; winD = gradient(win); % 窗函数一阶导数 tp = ((0:winLen-1) - Lh) / fs; % 以窗中心为零点的时间坐标 idx = (Lh+1) : (N-Lh); % 避开边界 nframe = length(idx); freq = (0:nfft-1) * fs / nfft; halfK = floor(nfft/2); TFR = zeros(halfK, nframe); tvec = (idx - 1) / fs; for m = 1:nframe n = idx(m); seg = x(n-Lh : n+Lh); S = fft(seg .* win, nfft); St = fft(seg .* winD, nfft); Sf = fft(seg .* win .* tp, nfft); t_n = (n - 1) / fs; for k = 1:halfK a = abs(S(k)); if a < gamma continue; end omega = freq(k) + imag(St(k) / S(k)) / (2*pi); tau = t_n - imag(Sf(k) / S(k)) / (2*pi); [~, kk] = min(abs(freq - omega)); [~, nn] = min(abs(tvec - tau)); TFR(kk, nn) = TFR(kk, nn) + a^2; end end end

调用主函数时，我用窗长winLen = 127，FFT 点数nfft = 1024，门槛gamma = 0.01 * max(abs(S(:)))。由于重分配输出时我用的是TFR(kk, nn)累加能量，最后用imagesc画图时通常还要取一下对数，否则动态范围太大，亮点会压掉弱分量。

3.4 关键参数怎么调

这套算法最关键的参数有三个：窗长、FFT 点数和能量门槛。窗长直接决定主瓣带宽和时间分辨率。我给出的 127 点高斯窗，在 1024 Hz 采样率下大约是 0.124 秒，这个宽度适合观察 150 Hz 到 450 Hz 频段的分量。如果两个分量频率差得特别近，我会把窗长加到 255 点甚至 511 点，让频率分辨率更高。

FFT 点数主要影响频率格点的细化程度。我习惯选不小于窗长的 2 的幂，比如 1024 或 2048。频率轴分辨率是fs/nfft，nfft 越大，重分配的目标格点越密，谱线定位越细腻，但耗时会增加。实测下来，对 1 秒信号、1024 点 FFT 的单核循环，跑一次大概需要十几秒，在调试阶段完全能接受。

能量门槛gamma是我建议优先调的参数。它的意义是避免把噪声能量随机分配到某个格点上，形成一堆虚假亮点。设太低，图面毛糙；设太高，真实小幅值分量会被削掉。我的经验是先不设门槛跑一次，看最大幅值量级，再取它的百分之一到千分之一作为门槛。噪声重的时候取百分之三左右比较稳。

3.5 结果验证与绘图

算法跑完，一定要用可视化验证，不要只看数值。我习惯把原始信号、标准 STFT 时频图和 IGD-SST 时频图画在同一个 figure 里对比。标准 STFT 直接用spectrogram画，IGD-SST 用imagesc画，坐标轴加上axis xy防止图像上下颠倒。

figure; subplot(2,2,1); plot(t, x); title('原始信号'); subplot(2,2,2); spectrogram(x, gausswin(winLen, 8), winLen-1, nfft, fs, 'yaxis'); title('短时傅里叶变换'); [TFR, freq, tvec] = igdsst(x, fs, 127, 1024, gamma); subplot(2,2,3); imagesc(tvec, freq(1:floor(nfft/2)), 10*log10(TFR + eps)); axis xy; title('改进群延迟同步压缩变换');

正常情况下，IGD-SST 图里的两个分量会呈现为两条细亮曲线，比 STFT 图清楚得多。如果那条线性调频分量在频率斜率大的位置仍然发虚，不要先去怀疑代码，先检查窗长是不是太短、门槛是不是太高。调试这类算法，最忌讳的就是不看中间结果直接改参数。

4. 常见问题与排查技巧

4.1 时频图上的能量不集中

这是一个出现频率非常高的问题。时频图上明明有两条分量，但线条很粗，或者局部发虚，通常不是算法错了，而是窗长和信号局部调频特性不匹配。我之前遇到一次，测试信号的调频斜率特别大，窗长却还是按低频信号的经验设成 127 点，结果压缩出来的线性调频分量在中段明显变粗，边缘还有轻微分叉。

解决办法是把窗长调短，或者改用自适应窗。调短窗会让时间分辨率变好，频率主瓣变宽，但只要后续的同步压缩能把频率方向能量收拢，整体分辨率还是能保住。如果两个分量频率差很近，不适合盲目调短窗，那就提高 FFT 点数和窗长，优先保住频率分辨率，牺牲一点时间定位。

4.2 两端出现明显的暗区和畸变

边界效应是老生常谈，但实际处理时很多人还是容易忽略。我的代码里用idx = (Lh+1) : (N-Lh)主动丢掉了边界帧，所以边缘畸变基本不会出现。不过丢边界帧也有代价：你丢失了信号开头和结尾的一小段时频信息。如果这段信息很重要，可以通过对信号做镜像延拓或对称延拓，把边缘补出来。

我自己在工程项目里更倾向保留边界帧但做加权补偿，计算同一帧在不同窗偏移下的重叠贡献，这样时间边缘信息会保留更多。这个方法实现起来要小心归一化，一不小心能量就失真。对大多数分析场景来说，直接丢边界帧是最省心最稳的做法。

4.3 多分量信号出现交叉干扰

当信号里有两个非平稳分量，瞬时频率曲线在时频平面上交叉时，交叉点附近的重分配容易把能量分错家，出现虚假的亮斑。这不是同步压缩变换独有的问题，任何重分配方法在交叉区域都会遇到所谓的“交叉项干扰”。IGD-SST 能缓解一部分，因为它在交叉附近同时用了频率修正和时间修正，两个曲线交叉瞬间的瞬时频率估计会相对准确一些，但做不到完全消除。

实际工程中我通常先做分量分离，比如用带通滤波或者脊提取，把两个分量大致分开，再分别做同步压缩分析。脊提取本身又会涉及瞬时频率追踪，可以简单粗暴地在时频图上选峰值点，也可以用动态规划找最优路径。等两条脊线都确定下来，把脊线附近的时频系数保留，别的地方清零，再做逆变换重构单分量，这样得到的时频图会非常干净。

4.4 运行速度慢和内存占用过高

IGD-SST 在 MATLAB 里最容易出现的性能瓶颈就是那层逐帧、逐频率点的双重循环。信号一长，循环次数急剧增加，跑半天出不来结果。我在调试时先用 1 秒信号把算法逻辑跑通，然后才放大到 10 秒甚至 60 秒信号。真到长信号阶段，我会把内层循环用向量化改写，一次处理一整帧的所有正频率点。

如果你安装了 Parallel Computing Toolbox，R2018A 支持parfor，把外层循环改成parfor就能明显提速。但要注意，parfor对变量切片有要求，TFR的累积更新需要改成局部累加再合并的方式，不然 MATLAB 会报错。另外nfft不要盲目拉到 8192，很多场景 1024 到 2048 已经足够，频率轴过于细密只会让计算量飞涨，视觉提升却很有限。

4.5 从 R2018A 迁移到新版 MATLAB 的坑

代码在新版本上运行，最常见的报错是函数名冲突或者参数语法变化。比如新版 MATLAB 的信号处理工具箱加入了stft函数，如果你之前的脚本里自己写过同名函数，新旧定义会打架。我自己的做法是给核心函数加上独立命名空间，或者直接放在单独的文件夹里，避免和工具箱函数重名。

还有一个容易踩的坑是chirp函数的参数行为在不同版本上有细微差别，如果测试信号是用chirp生成的，迁移后最好先对比一下信号波形是否一致。包括gausswin的参数单位、spectrogram的返回结果形式，这些细节在不同版本间都有演变。总体而言，IGD-SST 的核心代码只用fft、gradient这些非常稳定的基础函数，迁移起来比我想象中顺利得多。

5. 参考资料与扩展使用建议

5.1 值得去查的经典文献方向

做同步压缩变换，绕不开 Daubechies、Lu 和 Wu 在 2011 年提出的同步压缩小波变换那篇经典工作，它把经验模态分解的思路和重分配方法结合得非常漂亮。再往上追溯，Auger 和 Flandrin 在 1995 年提出的时频重分配方法也值得精读，群延迟估计和频率重分配准则的基本框架在那里已经有了完整表达。

群延迟估计方向，我建议去看重分配方法里相位导数的推导，以及它在非平稳信号处理中的应用文献。市面上很多关于时频重分配、同步压缩变换的综述都包含这一块。你不需要把每篇论文的公式都手推一遍，但至少要把重分配算子的推导逻辑看完，这样遇到参数异常时，你知道是该调窗长还是该调门槛，而不是盲目试。

5.2 还可以扩展哪些方向

这套方法做完之后，很自然的一个扩展方向是把群延迟估计和脊提取结合起来，做成自动瞬时频率提取工具。时频图只是视觉产物，真正用到故障诊断或者参数识别时，还需要把每条时频脊线提取出来，再计算脊线频率随时间的变化趋势。IGD-SST 由于能量集中度更高，提取出来的脊线会比标准 SST 稳定得多。

另一个值得尝试的方向是逆变换和信号重构。标准同步压缩变换因为保留了频率方向的压缩算子，可以做近似信号重构；IGD-SST 由于同时压缩了时间方向，重构的时候需要额外处理时间方向的补偿。我目前用的是一种简化策略：只保留脊线附近的压缩系数，把它们映射回原始的时频网格再做 ISTFT，效果在多数测试信号上都不错，但严格的理论重构分析还需要再推一轮公式。

最后分享一个我实际踩过的坑：第一次调通 IGD-SST 时，我把能量门槛设得特别低，整张图全是毛刺，差点以为算法本身有 bug。后来我把门槛设成最大幅值的百分之一，画面立刻干净很多。之后每换一批新数据，我都习惯先看一眼标准 STFT 的动态范围，再决定门槛取多少。这个习惯帮我避开了很多不必要的弯路，你也可以试试。