S变换原理与MATLAB实现：从时频分析到信号重建的实战指南

简介：一份面向电力系统与信号处理研究者的S变换MATLAB实现代码，聚焦电压暂降检测、谐波分析及频率动态追踪等实际问题。该代码基于Stockwell变换原理，以可调的尺度参数取代传统短时傅里叶变换的固定窗，能够自适应调整时频分辨率，尤其适合捕捉非平稳信号中的瞬态变化。包内仅含1个m源文件，整体2KB，结构紧凑、参数清晰，方便用户直接运行或将其嵌入已有的数据预处理流程；脚本可提取基频幅值、相位跳变、突变点及频率幅值包络线等关键特征，配合MATLAB的可视化能力可生成时频图像，直观展示电压暂降发生的时刻、幅度与扰动方向。已有629人学习下载，适用于从事电能质量分析、故障诊断及谐波源定位的工程师和研究人员快速上手与二次开发。

1. S变换到底是个什么东西：从时频分析的老问题说起

做信号处理的人，十有八九都遇到过这种尴尬：拿到一段非平稳信号，用短时傅里叶变换（STFT）分析，窗宽调来调去，要么频率分辨率上不去，要么时间定位模糊；换用小波变换，尺度参数又和实际频率对不上号，画出来的时频图横轴还得换算半天。S变换（S-Transform，简称ST）就是冲着这个痛点来的——它把短时傅里叶变换的固定窗函数换成了随频率变化的高斯窗，窗口宽度随频率自适应调整：低频时窗宽、频率分辨率高；高频时窗窄、时间定位准。这个特性让它在电力暂态分析、地震信号处理、心电信号降噪、雷达回波分析这些场景里特别吃香。

S变换的核心公式写出来并不复杂，它本质上是短时傅里叶变换的推广，也是小波变换的一种相位修正形式。在MATLAB里实现S变换，不需要装任何额外的工具箱，几十行代码就能跑通；但如果只知道调现成函数而不理解参数背后的意义，换一个信号就翻车也是常有的事。这篇笔记就从公式拆解说起，把S变换的MATLAB实现、参数设置、逆变换重建，以及那些你大概率会踩的坑一次讲清楚。适合正在做毕业设计、课程设计，或者刚接手非平稳信号分析项目的工程师照着复现。

2. S变换公式拆解：从连续定义到离散实现

2.1 连续S变换的数学定义与每个符号的含义

S变换的定义式通常写作：

[ S(\tau, f) = \int_{-\infty}^{\infty} x(t) \cdot \frac{|f|}{\sqrt{2\pi}} e^{-\frac{f^2(\tau-t)^2}{2}} e^{-i2\pi ft} dt ]

这个式子猛一看有点吓人，但拆开看其实就三块：原始信号 (x(t))、一个随频率变化的高斯窗 (w(\tau-t, f) = \frac{|f|}{\sqrt{2\pi}} e^{-\frac{f^2(\tau-t)^2}{2}})、以及复指数核 (e^{-i2\pi ft})。这里的 (\tau) 是时间平移因子，控制窗函数在时间轴上的位置；(f) 是频率，直接决定窗宽——频率越高，高斯窗的标准差 (\frac{1}{|f|}) 越小，窗就越窄，时间分辨能力越强；频率越低，窗就越宽，频率分辨能力越强。这正好解决了STFT窗宽无法自适应的问题。

其中时间和频率是完全耦合的，不是简单的两个独立参数。理解这一点非常关键：S变换本质上是一组带宽随中心频率线性变化的小波滤波器组，但它比小波变换多了一个相位因子 (e^{-i2\pi ft})，这让S变换的结果和傅里叶变换之间保持了直接的联系——对S变换结果沿时间轴积分，就能得到原始信号的标准傅里叶谱 (X(f))。这个性质叫"无损可逆性"或者"边际性质"，是做逆变换和信号重建的理论基础。

2.2 离散S变换公式与MATLAB的向量化思路

实际在MATLAB里处理信号，信号是离散的采样点序列，所以必须用离散S变换（DST）。标准离散定义的推导思路是先用傅里叶变换改写S变换：令 (X(\alpha)) 为 (x(t)) 的傅里叶变换，则S变换可以写成频域形式：

[ S(\tau, f) = \int_{-\infty}^{\infty} X(\alpha + f) \cdot e^{-\frac{4\pi^2\alpha^2}{2f^2}} \cdot e^{i2\pi\alpha\tau} d\alpha ]

写成离散版本之后，就发现一个很关键的计算结构：对每个离散频率 (f)，把平移后的频谱 (X(\alpha+f)) 和高斯窗相乘，再做一次逆傅里叶变换，就得到该频率对应的时间序列。这意味着S变换的MATLAB实现可以完全基于fft和ifft完成，不需要任何循环卷积或者矩阵求逆操作，效率非常高。

我一般会直接按向量化写法实现，核心代码如下：

function ST = s_transform(x, fs, fmin, fmax, df) % S变换的向量化实现 % 输入：x - 输入信号（行向量） % fs - 采样率 % fmin, fmax - 分析的频率范围 % df - 频率分辨率 % 输出：ST - 复矩阵，行对应频率，列对应时间 N = length(x); X = fft(x); % 原始信号的FFT f = fmin:df:fmax; % 频率轴 t = (0:N-1) / fs; % 时间轴 nf = length(f); ST = zeros(nf, N); % 预分配复数矩阵 for k = 1:nf fk = f(k); % 频率平移：在频域做循环移位 X_shift = circshift(X, [0, round(fk / fs * N)]); % 高斯窗（频域形式） gauss = exp(-2 * pi^2 * (fk / fs .* (-N/2:N/2-1)).^2 / fk^2); % 乘以窗函数后再做逆FFT回到时域 ST(k, :) = ifft(X_shift .* gauss); end end

这段代码的向量化关键点在于：第一，平移频谱用circshift完成循环移位，避免手动写索引拼接的麻烦；第二，高斯窗直接用频域自变量生成，省掉了一重时域到频域转换；第三，所有运算都基于FFT，复杂度是 (O(N_f \cdot N \log N))，比直接按时域定义做卷积快一个数量级。需要注意的是，circshift的移位量必须换算成点数，即round(fk / fs * N)，如果fk / fs * N不是整数，四舍五入会带来轻微的频率误差，这在后面讲精度的坑时会详细说。

2.3 为什么用FFT实现：从计算效率到内存占用的实际考量

有人可能会问：既然定义式里直接就是时域积分，为什么非要用频域实现？直接双重循环不也能算吗？能算，但代价天差地别。按时域定义，每个频率点都要做一次全长的卷积运算，假设信号长度是 (N=4096)，频率点数是 (200)，那就是 (200 \times 4096^2) 约等于 33 亿次复数乘加，MATLAB跑下来要好几分钟；而频域实现利用FFT，同样的参数只需要 (200 \times (2 \times 4096 \log_2 4096)) 大约两千万次运算，秒级甚至毫秒级完成。这是把高斯窗先变换到频域再和信号频谱相乘的经典滤波思路，凡是做过快速卷积的人应该都很眼熟。

内存方面也要提前心里有数。S变换的输出是一个 (N_f \times N) 的复数矩阵，每个复数元素占16字节。比如采样率 (fs=1000\text{Hz})，信号时长10秒，(N=10000)，频率分析到 (f_{max}=500\text{Hz})，步长 (df=1\text{Hz})，那矩阵就是 (500 \times 10000)，占内存约 (500 \times 10000 \times 16 = 80\text{MB})。看起来还好，但如果把采样率提到10k，信号时长拉长到30秒，频率分辨率做到0.5Hz，矩阵规模就到了 (20000 \times 300000)，也就是96GB，直接爆内存。所以实际工程中，频率范围、频率步长和信号长度三者之间要做取舍，CPU能等但内存等不起。

3. 用MATLAB把S变换跑起来：最小实现与可视化

3.1 生成测试信号并调用S变换：一段完整可跑的脚本

首先生成一个包含频率突变和调频分量的复合测试信号，这是验证S变换效果的标准做法。信号包含三段：0~2秒是50Hz正弦波，2~4秒是120Hz正弦波并叠加频率从20Hz线性扫到80Hz的chirp信号，4~6秒是80Hz正弦波加微弱的高斯白噪声。用这样的信号能同时看出S变换的频率分辨能力和时间定位能力。

%% 生成非平稳测试信号 fs = 500; % 采样率 500Hz T = 6; % 总时长 6秒 N = fs * T; % 总采样点 t = (0:N-1) / fs; % 时间轴 % 三段信号拼起来 x1 = sin(2 * pi * 50 * t(1:fs*2)); x2 = sin(2 * pi * 120 * t(fs*2+1:fs*4)) + ... chirp(t(fs*2+1:fs*4), 20, 4, 80, 'linear'); x3 = sin(2 * pi * 80 * t(fs*4+1:fs*6)) + ... 0.1 * randn(1, fs*2); x = [x1, x2, x3]; %% 计算S变换 fmin = 0.5; fmax = 200; df = 0.5; ST = s_transform(x, fs, fmin, fmax, df); %% 可视化时频图 f_axis = fmin:df:fmax; figure; imagesc(t, f_axis, abs(ST)); axis xy; % y轴从下往上显示低频 xlabel('时间 (s)'); ylabel('频率 (Hz)'); title('S变换时频幅度谱'); colorbar; colormap(jet);

运行这段脚本的前提是 2.2 节的s_transform函数文件已经保存到当前工作路径或MATLAB搜索路径中。imagesc默认会缩放颜色映射，幅度小的分量在图上会显得很淡，如果细节看不清楚，可以手动调整caxis的范围，比如只显示幅度最大的前20%到80%的能量，突出主要分量。实际效果是：50Hz段显示一条均匀的亮线，2秒切换时刻的边界非常陡峭，2到4秒之间的chirp是一条从20Hz斜着上升的亮带，噪声背景呈细碎的颗粒状，这个"颗粒感"大小直观反映了S变换在高频段的时间分辨率和低频段的频率分辨率。

3.2 从时频矩阵提取瞬时频率曲线：幅度阈值法

S变换输出的复数矩阵不只是拿来画图的。最常用的一个下游操作是提取某个时刻的主瞬时频率，做法是找每个时间点上幅度最大的频率分量。但直接取max容易被宽频噪声带偏，尤其噪声功率低但分布广时，噪声峰完全可能压过信号峰。所以一般做法是先设一个幅度阈值，例如取全局最大幅度的20%以下的部分置零，再逐列找峰值。

%% 从ST矩阵提取瞬时主频率 A = abs(ST); threshold = 0.2 * max(A(:)); % 全局阈值 20% A_thresh = A; A_thresh(A < threshold) = 0; % 逐列找最大幅度对应的频率 [~, idx_max] = max(A_thresh, [], 1); inst_freq = f_axis(idx_max); % 画主频曲线和时频图叠加 hold on; plot(t, inst_freq, 'k', 'LineWidth', 2);

这里有个更妙的做法：S变换的复数谱其实保留了相位信息，而瞬时频率还可以从相位对时间的导数求出来——对每个时间点，把相邻两个时间点的相位差铺开，除以 (2\pi \Delta t) 就是更精细的瞬时频率估计，分辨率可以远高于频率步长 (df)。这种方法在电机转速分析、电力系统谐波监测这些需要高精度频率估计的场景中非常实用。不过相位解卷绕（unwrapping）是最容易出错的环节，一个跳变的 (2\pi) 毛刺就足以让整条曲线看着像纯随机噪声，下一节专门讲这个坑。

3.3 频段能量占比：定位扰动或故障时刻的最简单指标

还有一种常见分析需求：不需要知道具体频率是多少，只需要判断某个频带内的能量什么时候变强了。典型场景比如输电线路的暂态扰动检测——发生故障瞬间，高频分量会突然爆发。计算某个频带内的S变换幅度平方和随时间的变化，就是一个非常灵敏的检测器：

%% 计算指定频带内的时变能量 band_idx = find(f_axis >= 50 & f_axis <= 150); % 50~150Hz频带 energy = sum(abs(ST(band_idx, :)).^2, 1) / length(band_idx); % 归一化 energy = energy / max(energy); % 找到能量突变的时刻（超过均值+3倍标准差） mu = mean(energy); sigma = std(energy); event_idx = find(energy > mu + 3 * sigma); if ~isempty(event_idx) event_time = t(event_idx(1)); fprintf('检测到能量突变，起始时刻: %.3f s\n', event_time); end

这个方法比直接在原始信号上算滑动方差的好处是它能锁定频带，不会把宽带噪声混进来，也不会因为某个无关的低频缓变分量被误判。把能量曲线画出来叠加在原信号下方，做故障定位时一眼就能看出问题时刻。如果能量曲线频繁超过阈值，通常意味着要么信号里包含周期性冲击，要么频带选得太窄被窄带干扰穿过，这时要回头检查频带边界是否选得合理。

3.4 批量处理多个信号文件：循环与参数统一

做实验时经常遇到几百个信号文件要统一跑S变换。我的习惯是不写复杂的数据导入层，直接把参数统一在一个结构体里，循环处理并保存结果。这样后续调整参数只需要改一处，不会出现文件A和文件B用了不同频率分辨率导致结果没法对比的尴尬。

%% 批量处理多个信号文件 files = dir('data/*.mat'); % 读取 data 文件夹下所有mat文件 params.fs = 500; params.fmin = 0.5; params.fmax = 200; params.df = 0.5; results = struct(); for k = 1:length(files) load(fullfile('data', files(k).name)); % 假设每个 mat 文件里都有一个变量 x results(k).name = files(k).name; results(k).ST = s_transform(x, params.fs, params.fmin, ... params.fmax, params.df); results(k).freq = params.fmin:params.df:params.fmax; results(k).time = (0:length(x)-1) / params.fs; fprintf('已处理 %d / %d\n', k, length(files)); end save('s_transform_results.mat', 'results');

这里注意load之后变量x是直接进到工作区的，如果多个文件的变量名不统一，比如有的存成signal有的存成data，就要在load之后用fieldnames检查再做赋值。批量处理时内存也是一道坎——results结构体里存了完整的ST矩阵，几百个文件下来可能几个GB，如果只需要提取瞬时频率或者频带能量，建议在循环里直接做特征提取，只存标量或者曲线，别存矩阵。这是我踩过最痛的坑之一，跑了一夜的分析因为内存溢出中断，回头看代码全是在存储上扣的细节。

4. 逆S变换与信号重建：能不能把原信号完美拿回来

4.1 逆变换公式与MATLAB实现：从时频图还原信号

S变换有一个非常重要的数学性质：时频矩阵沿时间轴积分回到傅里叶谱，这意味着它存在一个干净的逆变换。逆变换公式是：

[ x(t) = \int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} S(\tau, f) d\tau \right) e^{i2\pi ft} df ]

也就是说，先把S变换矩阵对时间轴求和（积分），得到每个频率对应的傅里叶谱分量，再对这个频谱做逆傅里叶变换就还原了原始信号。这比短时傅里叶变换的逆变换简单——STFT重建需要额外做窗函数归一化甚至重构窗设计，而S变换天然就是精确可逆的。

MATLAB实现代码非常短：

function x_rec = inverse_s_transform(ST, fmin, fmax, df) % 逆S变换：对时间轴积分再逆FFT % 输入：ST - S变换复数矩阵，nf行 × N列 % fmin, fmax, df - 与正变换时的频率轴参数保持一致 % 输出：x_rec - 重建的时域信号（复数） % 沿时间轴积分（求和），等价于对 τ 做积分 X_hat = sum(ST, 2)'; % 得到与频率相关的谱 % 注意：如果正变换时用了频域高斯窗，这里得到的谱就是 x(t) 的FFT x_rec = ifft(X_hat) * (fmax - fmin + df); % 缩放系数 end

这里有个很隐蔽的缩放系数问题。如果正变换按 2.2 节的向量化写法实现，那么逆变换时积分结果和原始FFT之间会有一个与频率步长相关的常数倍关系，必须乘回(fmax - fmin + df)或者df才能让重建信号的幅值恢复到原始水平。这个系数到底取多少，取决于正变换时频率轴是怎么采样的——如果fmin:df:fmax本身不是从0到fs均匀覆盖，逆变换结果会带上系统偏置。我建议的做法是：生成一段已知信号，跑正变换再跑逆变换，用rms(x - x_rec)检查重建误差，然后根据误差量级调整缩放因子，不用纠结理论推导，实测校准最靠谱。

4.2 重建误差的验证方法：RMS与信噪比

在实际使用中，逆变换最常见的用途是滤波。思路是：对S变换矩阵做掩膜处理——把想要保留的时频区域的元素保留下来，其他区域置零，然后逆变换回时域，这样就实现了时频域滤波。比如信号有一段50Hz工频干扰，不需要滤掉的频段以外全部置零，重建出来的信号就是干净的。

%% S变换时频滤波：保留指定区域，其余置零 mask = zeros(size(ST)); % 保留 48~52Hz 全时段（例如工频信号） idx_f = find(f_axis >= 48 & f_axis <= 52); mask(idx_f, :) = 1; % 滤波后信号 ST_filt = ST .* mask; x_filt = inverse_s_transform(ST_filt, fmin, fmax, df); % 重建误差评估 err = rms(real(x_filt) - x); snr_dB = 10 * log10(sum(x.^2) / sum((real(x_filt) - x).^2)); fprintf('重建信噪比: %.2f dB\n', snr_dB);

做这种时频掩膜滤波时，一个常见的坑是：掩膜边界太锐利，导致重建信号在边界处出现吉布斯振铃。改善办法是对掩膜做平滑处理，比如用imgaussfilt或者简单的一个filter2平均滤波器，让掩膜边缘有过渡带。另一个问题是，如果掩膜把某些频率分量完全去掉了，重建信号的直流偏置会漂移，尤其在信号本来就有非零均值时非常明显。这时要么在掩膜之外保留低频零频附近的分量，要么重建后在时域减去均值。

4.3 边界效应：为什么两端总是重建不准

逆变换重建实验做得多了就会发现一个规律：重建信号的头尾总是不如中间准。这不是MATLAB实现的问题，而是S变换本身的边界效应。原因有两个层面。第一，正变换时对信号做FFT默认信号是周期的，但实际信号两端并不连续，频谱泄漏导致两端时频分布异常；第二，高斯窗在时域是无限长的，做离散化时频率越低窗越宽，越容易跨出信号边界，导致边界点附近的有效信息缺失。

想缓解边界效应，有几个土办法但都很好用。最简单的是在正变换之前把信号两端各扩展一段数据（镜面反射扩展），变换完再裁掉；第二种是把信号两端加窗压到接近零，但这样会损失边界处的真实成分；第三种是用pchip插值把序列延长为原来的两倍，算完S变换再只取中间部分。我的个人习惯是镜面反射扩展，代码如下：

%% 镜面反射扩展消除边界效应 ext_len = round(0.1 * length(x)); % 扩展10%长度 x_ext = [fliplr(x(1:ext_len)), x, fliplr(x(end-ext_len+1:end))]; % 对扩展信号做S变换 ST_ext = s_transform(x_ext, fs, fmin, fmax, df); % 只取中间原始长度部分 ST = ST_ext(:, ext_len+1:end-ext_len);

实战下来，这个处理能把边界处的重建误差从几个百分点降到千分之一以下。代价是频率分辨率略有损失（因为有效分析窗变短），但绝大多数场合这点损失完全不影响结论。

5. 参数怎么设：fmin、fmax、df、窗宽的玄学与经验

5.1 频率范围选择：fmax不要等于fs/2，留15%余量

很多人做S变换的第一步就是fmax = fs/2，觉得采样的奈奎斯特频率就是分析上限。这个直觉在FFT频谱分析里是对的，但在S变换里要打个折扣。因为S变换每个频率点用的高斯窗频域宽度和频率本身有关，最高频率点的窗宽会让它实际上采样到超出 (fs/2) 的分量，这些镜像成分折叠回来，会在最高频段造成虚假能量。

我一般取fmax = 0.42 * fs左右，或者直接取信号实际有效带宽的1.2倍。比如采样率fs=1000Hz，信号的主要能量在200Hz以下，那设fmax=240Hz就够用了，没必要冲到500Hz。频率上限越高，矩阵行数越多，计算时间线性增长，但信息增益非常有限——高频段的S变换幅度通常很小，全是噪声和窗函数的边角料。频率下限fmin也不是越小越好的哲学。设fmin=0意味着要分析零频分量，但直流分量对S变换的时频图贡献是一条全时间均匀的细线，没有任何有用信息，白白占一行矩阵。所以fmin一般从1Hz或0.5Hz起步，除非你要分析慢变趋势。

5.2 df的选择：频率分辨率与计算负担的平衡

df直接决定了输出矩阵的列数，也就等于频率点的数量。设0.1Hz意味着输出矩阵有几千行，计算时间变成几十倍，但时频图在频率轴上的平滑度几乎没差别，因为S变换本身的频率分辨率由高斯窗决定，不是由采样步长决定。换句话说，df只是插值密度，不是真正的分辨率。真正的频率分辨率大约等于中心频率除以一个常数，低频段能分辨两三赫兹的差异，高频段差两三赫兹根本分不开。

所以我的经验是：df设为信号预期最低频率分辨率的1/3到1/5就好。分析电力谐波时，基频50Hz，谐波间隔50Hz，那df=5Hz就足够；分析地震波这种低频信号，频率分辨率要精细到0.1Hz的量级，df=0.1Hz才合理。如果你不确定该取多少，一个保守的做法是从df=1Hz开始跑，看时频图够不够细腻，如果频率方向上细节已经糊成一团，再减小df；如果频率方向上细得过分但计算时间明显拉长，那就增大df。不要一上来就拍脑袋设df=0.01，那会让计算时间膨胀到让人怀疑人生的程度。

5.3 高斯窗宽度参数：要不要自定义修改

标准S变换的高斯窗宽度完全由频率决定，它的标准差是 (\frac{1}{|f|})，没有任何可调参数。但有些信号用标准窗分析效果并不理想：比如分析持续时间极短的脉冲信号，标准窗在低频段太宽，脉冲周围的时间定位能力不够；分析宽带噪声背景下的窄带信号，标准窗又显得太窄，频率分辨力不够。这时候就有改进方案，常见的一种是引入一个可调宽度因子 (k)：

[ w(\tau-t, f) = \frac{|f|}{k\sqrt{2\pi}} e^{-\frac{f^2(\tau-t)^2}{2k^2}} ]

当 (k>1) 时窗变宽，频率分辨率提升、时间分辨率降低；(k<1) 时窗变窄，时间分辨率提升、频率分辨率降低。这个参数相当于在小波分析和STFT之间做插值——(k) 极大时接近STFT的固定窗效果，(k) 极小时接近小波变换的高频细节提取。实际实现只需要在代码中把高斯窗的fk^2分母换成k^2再乘回去。如果你想自定义这个参数，相当于做广义S变换，比标准S变换灵活得多。但要注意：修改窗宽后，逆变换的精确重建性质不再自动成立，如果后面要做时频滤波再重建，就得多扛一个误差。

6. S变换的五个高频踩坑记录：现象、原因与解决办法

6.1 时频图有水平方向的条纹噪声，看起来像摩尔纹

现象：画出的S变换时频图在低频段有明显的水平亮条纹，跟电视机的摩尔纹似的，和高频段的细碎结构完全不是一个风格。

原因：90%的情况是MATLAB的imagesc默认对每一行做独立的颜色归一化，或者colormap的动态范围过宽，导致低频部分幅度虽然很小，但相对变化却被色带放大。还有一些情况是真的算错了——circshift的频率平移量没用round而是用了fix，在负频率方向上截断了频谱。

解决：先检查数据本身，画plot(f_axis, mean(abs(ST),2))看每行平均幅度是不是平滑递减，如果曲线波动剧烈就重新检查正变换逻辑；如果数据没问题，那就把imagesc的clim手动设为比如[0, max(abs(ST(:))) * 0.6]，或者改用pcolor加shading interp平滑显示。

6.2 逆变换后信号幅值整体偏移，直流分量全丢了

现象：时频滤波后重建的信号波形形状对，但整体被抬高或压低了一个常量，直流偏置不对。

原因：滤波时把0Hz附近的掩膜置零了。直流分量在S变换中对应的是频率轴第一行（或者fmin以下的部分），如果fmin从1Hz起，很多实现根本不计算零频，逆变换时这个常量成分自然就丢了。

解决：滤波时不要把最底下两行全置零，留出一个1~3Hz宽的低频通道；或者重建后做一次去均值校准，把原始信号均值加回去。做检测场景不受影响，但如果后续要做幅值定量分析，这个问题必须处理。

6.3 高频段全是噪点，信号成分反而看不清

现象：fmax设到接近fs/2，高频段一片密密麻麻的亮点，完全看不出任何与信号相关的结构。

原因：高频段S变换的时间分辨率高，但频率分辨率差，而噪声在时频域里本来就到处都是孤立的高幅值点。加上前面说的镜像折叠，超出奈奎斯特频率的分量也会在这里叠加，噪点密度进一步提高。

解决：把fmax降下来，设到目标信号的最高有效频率乘以1.2到1.5倍就够了。如果确实需要看很高的频段，改用更高采样率重新采集信号，别指望算法把高频信息凭空变出来。噪点实在碍事，可以做一个3×3的中值滤波核先对幅度矩阵做平滑，代价是会轻微降低时间定位精度。

6.4 S变换计算巨慢，跑一次要几分钟

现象：同样的信号长度，STFT眨眼就算完，S变换却要跑半天。

原因：频率点数太多或者代码里用了标量循环没用向量化。最常见的是有人在定义式的时域积分上直接写双重for循环，或者对每个频率点都调用一次fft而且没有预分配输出矩阵。MATLAB里矩阵在循环中动态增长是最伤性能的操作。

解决：按 2.2 节的向量化实现走，一次性计算整行。检查循环前是否做了ST = zeros(nf, N)预分配。更进一步，可以用parfor把频率循环并行化，前提是每个循环体的ifft输入独立，没有共享变量。我实际测过，8核机器上并行化能再快3倍左右。还要检查信号长度是不是2的幂次，如果不是，fft会回退到混合基算法，速度掉一个量级，用nextpow2补零到2的幂很划算。

6.5 提取的瞬时主频率曲线有大量毛刺，时跳时稳

现象：瞬时频率曲线前半段平直，后半段突然上下乱跳，或者整条曲线都像牙签盒倒出来一样。

原因：幅度阈值选得太低或根本没设。噪声频点在某几个时间点上恰好比信号频点幅度大，峰值索引就跳到噪声频率上去了。曲线在低幅区间跳变尤其严重，因为信号本身在那段时间里可能很弱甚至没有，任何噪声都占主导。

解决：设完整的两层筛选：第一层用全局幅度阈值筛掉低幅噪声；第二层用频率连续性约束，如果后一个时间点的主频比前一个点的主频突变超过一个预设值（比如20Hz），就认为这是错误跳变，用前一时刻的频率代替或做滑动中值滤波。曲线提取这种活儿要的就是老老实实做后处理，一步都不能省。

7. 进阶技巧：用S变换做交叉项抑制和信号分离的边界条件

S变换还能干一些不太为人知但很有用的活儿：对多分量信号做交叉项抑制（对比Wigner-Ville分布那种严重交叉项），以及时频域盲信号分离。这两个方向并不需要改S变换本身，只需要用好它的相位信息和矩阵结构。

先说交叉项抑制的逻辑。Wigner-Ville分布对两个叠加的正弦信号会产生严重的交叉项，而S变换的时频矩阵是线性的，没有交叉项。但S变换也有自己的短板——低频段两个靠得很近的分量分辨不开。有个技巧是把S变换矩阵和它的瞬时频率估计结合，对每个时间点做局部峰值聚类，找到三个以上的局部峰时，用相位信息判断这几个峰是否属于同一个源的谐波成分（频率比为整数倍），如果是，把主峰的瞬时频率作为该源的基频，再把对应谐波所在时频区域用掩膜分离出来。这个过程可以循环做到每个源都分离干净。

%% 基于S变换的时频域分离示例：分离两个叠加的chirp信号 % 假设ST已经是两个chirp信号的S变换结果 A = abs(ST); f_axis = params.fmin:params.df:params.fmax; t_axis = (0:size(ST,2)-1) / fs; % 在每个时间点找局部峰值（至少间隔5个频率点） n_peaks = 3; peak_freqs = zeros(n_peaks, size(ST,2)); for col = 1:size(ST,2) amp_col = A(:, col); [peaks, locs] = findpeaks(amp_col, 'MinPeakDistance', 5, ... 'SortStr', 'descend', 'NPeaks', n_peaks); if length(locs) >= 1 peak_freqs(1:length(locs), col) = f_axis(locs); end end % 用连续性追踪两条主频率曲线 % （这里简化为用最大峰和次大峰的频率画两条曲线） primary = peak_freqs(1, :); secondary = peak_freqs(2, :); % 根据两条频率线构造两组掩膜 mask1 = zeros(size(ST)); mask2 = zeros(size(ST)); for col = 1:size(ST,2) if primary(col) > 0 [~, idx1] = min(abs(f_axis - primary(col))); mask1(max(1,idx1-3):min(end,idx1+3), col) = 1; end if secondary(col) > 0 [~, idx2] = min(abs(f_axis - secondary(col))); mask2(max(1,idx2-3):min(end,idx2+3), col) = 1; end end % 分别重建 x1_sep = real(inverse_s_transform(ST .* mask1, fmin, fmax, df)); x2_sep = real(inverse_s_transform(ST .* mask2, fmin, fmax, df));

这段代码的思路是先用findpeaks在每个时间点找到幅度最大的几个频率峰，再根据频率连续性做追踪，最后构造窄带掩膜重建出每个分量。实际操作中，findpeaks的MinPeakDistance和NPeaks参数是玄学，设太小会把一个宽峰的旁瓣都当成独立峰，设太大真的峰值又丢了。我的经验是：先画一列看典型的峰宽度，再回填参数。如果两个源的频率在某个时刻重叠或者靠得太近，这个简单方法会失效，那时可以退化到先验信息——比如知道一个信号是50Hz工频、另一个是频率随时间变化的干扰，就只对工频线做窄带跟踪，剩下的空间全留给另一个信号。

做这类时频域分离有个边界条件必须心里有数：只有当两个分量的时频支撑区域在大多数时间点上不重叠时，分离才干净。如果两个分量频率全程都一样、相位不同，S变换也救不了你——这在数学上本就是欠定问题。所以别指望一个算法通吃所有分离需求，先算一下时频图的支撑重叠率，再决定要不要上这套方法。

另一个进阶方向是把逆S变换用于时频域自适应阈值降噪。做法是：把ST矩阵的幅度做统计建模，假设噪声服从瑞利分布（幅度域的加性高斯噪声），用中值绝对偏差估计噪声水平，然后做软阈值或硬阈值收缩，再逆变换重建。这本质上是一种时频域维纳滤波的简化版。阈值参数lambda的经验公式是 (\lambda = \sigma_n \sqrt{2 \log(N)})，其中 (\sigma_n) 是中值绝对偏差方法估计出来的噪声标准差。MATLAB里一行：

sigma_n = median(abs(ST(:))) / 0.6745; % MAD鲁棒估计 lambda = sigma_n * sqrt(2 * log(N)); ST_denoise = ST .* max(0, 1 - lambda ./ abs(ST)); % 软阈值

软阈值的重建结果比硬阈值平滑，没有恼人的伪吉布斯振荡，这是我在心电信号降噪上反复对比后的结论。但要注意，这个阈值策略对每个时频点独立操作，忽略了邻近点的相关性。如果噪声是时变的，阈值应该跟着噪声水平做二维渐变，而不是全局一个常数。

最后说一个我自己的习惯：每当拿到一个新信号，先快速跑一遍S变换，把时频图最亮的几个区域圈出来，分别逆变换重建，听一听或者看一眼波形，这个过程能极大加深对信号结构的直觉。S变换的矩阵里存着的不只是能量分布，还有相位，相位里藏着频率微变的痕迹。这些细节，光看幅度谱永远是发现不了的。希望这篇实战笔记能帮你在面对非平稳信号时少走弯路——S变换不是万能的，但在该用它的场景里，它确实比STFT和小波都更顺手。

本文还有配套的精品资源，点击获取