S变换原理与MATLAB实现：从STFT到自适应时频分析

简介：这是一份面向电力系统信号处理与故障分析研究的MATLAB源码资源，围绕S变换（Stockwell变换）实现时频分析，重点用于电压暂降、谐波和电压频率检测；S变换在短时傅里叶变换基础上引入尺度因子，可自适应调整时频分辨率，适合分析电压暂降这类非平稳瞬态事件。资源压缩包仅2KB，内含1个.m文件，代码精简、无额外依赖，导入MATLAB即可运行；已有629人学习浏览，说明该主题受到较多关注。通过运行源码，读者可以直接获取S变换公式的实现流程，并在电压暂降数据上提取基频幅值、相位跳变、突变点及频率幅值包络线等关键特征，生成时频可视化结果，为故障定位、谐波源识别和电能质量分析提供有力支撑。整体来看，这份资源兼顾了算法原理与工程应用，适合电气工程、信号处理方向的学生和工程师作为入门与实践参考。

1. S变换：它在什么场景下值得你动手写一遍

做信号时频分析的人，大概率都经历过这个场景：手里一段振动或电能质量波形，既想知道它从哪一瞬间开始变了，又想知道变化时主要是哪些频率在起作用。S变换（S-transform，也叫Stockwell变换）正好卡在这个需求上：它不像FFT那样丢掉时间信息，也不像短时傅里叶那样用固定窗两头将就，而是用一组随频率自动缩放的高斯窗，把频率和时间同时保留在一张二维时频图里。公式不复杂，MATLAB实现也只需要几十行，但它能解决很多实际落地问题。

这篇文章按公式拆解、MATLAB实现、参数定标、常见问题、正确性验证的顺序展开，适合做电力暂态分析、机械故障诊断、地震信号处理、生理信号分析的人直接照着做。下面每一步都给出可复制的代码和参数说明，新手能跟到出图，熟手可以直接跳到第4章看参数边界和第5章看坑。

2. 把S变换公式拆到能写代码的程度：从连续积分到离散求和

2.1 窗函数随频率自适应，是S变换和短时傅里叶的分水岭

短时傅里叶变换（STFT）做的事情很简单：把信号切成固定宽度的小段，每段做一次FFT。这个固定宽度就是它的命门——窗短了，频率分辨率差，两个靠得近的频率会糊在一起；窗长了，时间分辨率差，信号瞬态变化被平均掉。无论你怎么调窗长，全频段用的都是同一个分辨率，这就是STFT和S变换最本质的差别。

S变换的连续公式长这样：

S(τ, f) = ∫ x(t) · w(τ−t, f) · e^(−j2πft) dt

其中窗函数：

w(t, f) = |f| / (√(2π)) · e^(−t²f²/2)

重点全在窗函数上：频率f越高，高斯窗在时间轴上越窄，时间定位越准；频率f越低，窗在时间轴上越宽，频率分辨率越好。从短时傅里叶的“一把尺子量所有频段”，变成高频看细节、低频看整体的自适应尺子，这就是S变换公式里窗函数带来的核心价值。

这里有个容易忽略的点：连续公式里的窗面积对每个频率都是归一化的，所以S变换结果复矩阵沿时间轴积分后能还原出原始信号频谱。很多人在MATLAB里把S变换当黑匣子用，根本没想过它自带一个可验证的性质。这个性质在第6章的验证脚本里会用到，先记住“可逆”两个字。

2.2 离散形式：真正写进m文件的是这个求和式

实际工程里信号是采样得到的离散序列x[n]，连续积分要翻译成离散求和。如果用暴力方式直接算双重循环，慢到没法用。常见做法是先用FFT算出整段信号的频谱X[m]，再对每个频率点n做一次频域搬移和加权，最后逆FFT回到时间域。离散形式的标准写法是：

S[m, n] = Σ(k=0 to N−1) X[(n + k) mod N] · e^(−2π²k²/n²) · e^(j2πkm/N)

这个式子看着复杂，拆开看就三部分：第一部分X[(n+k) mod N]是原始频谱搬移n个位置，相当于把以n为中心的邻域频段取出来；第二部分e^(−2π²k²/n²)是高斯窗在频域的形状，它的宽度随n变化，这就是上一节那个自适应窗的离散化身；第三部分e^(j2πkm/N)是逆变换的旋转因子，把频域拉回时间域。

写代码时需要注意，n对应的是频率索引，真实频率是n·Fs/N，其中Fs是采样率、N是信号长度。另一个很重要的点是n从1开始而不是从0开始，因为n=0时高斯窗指数项的分母为0，直接算出NaN或Inf。凡是S变换一开始就出现一片NaN的人，八成都是没处理直流分量这个问题。

2.3 三种时频工具怎么选：S变换不是万能的

很多初学者把S变换理解成“更好的短时傅里叶”，实际不是。它是另一种选型，有自己的边界和代价。放一张对比表方便你决定当前项目该用谁：

从实际项目选型看：机械故障诊断里想提取某个特征频段的时域波形，S变换逆变换方便，适合；电网暂态分析要同时看电压暂降起止时刻和谐波成分，S变换也顺手；但如果你要做的是音频实时频谱分析，每个频点都要一次FFT，S变换比STFT贵不少，未必划算。下面进入MATLAB实现，把2.2的公式变成能跑的代码。

3. MATLAB实现S变换：主函数代码与调用脚本

3.1 最简版S变换函数：时频卷积法

先给一个可以直接复制的最简实现。它不做任何花哨优化，只忠实还原2.2的离散求和式，方便你理解每一步在算什么。

function ST = s_transform_simple(x, fmin, fmax) % S变换最简实现：对每个频点做频谱搬移 + 高斯窗加权 % 输入 x 为一维信号，行向量或列向量均可 % 输出 ST 矩阵：行 = 频率索引，列 = 时间 x = x(:); % 统一成列向量 N = length(x); X = fft(x, N); % 整段信号的FFT fmin = max(1, round(fmin)); % 0频点公式不可用，从1开始 fmax = min(floor(N/2) + 1, round(fmax)); numF = fmax - fmin + 1; ST = zeros(numF, N); % 预分配结果矩阵 for n = fmin : fmax % 频谱搬移：得到 X[(n+k) mod N]，对应公式第一项 Xshift = circshift(X, -n); % 高斯窗的频域形式，n 是当前频率索引 k = (0 : N-1).'; gauss = exp(-2 * pi^2 * k.^2 / n^2); % 加权后逆FFT，得到该频率下的时间分布 ST(n - fmin + 1, :) = ifft(Xshift .* gauss, N); end end

代码逻辑说明：主函数对整段信号做一次N点FFT得到频谱X。for循环里，对每个频率索引n，用circshift(X, -n)把频谱向左搬移n个位置，等效于取出以n为中心的那一段频谱。然后用该频点对应的高斯窗去加权相邻频带，窗的宽度随n变化，这正是S变换区别于STFT的关键行。最后ifft把加权后的频谱拉回时间域，这一行就是一整个频率点的时间演化。

参数说明：fmin和fmax是频率索引，不是真实Hz。调用前需要按n = round(f / Fs * N)把物理频率换算成索引。fmax建议不要超过floor(N/2)+1，那是奈奎斯特频率对应的索引，更大的索引对应镜像频率，计算了也没有额外物理信息。x(:)这行很重要，它把输入统一成列向量，避免你传一个行向量、函数输出一个转置后的矩阵，最后画图时上下颠倒找不到原因。

3.2 只计算目标频段：把fmin/fmax当作第一性能开关

标准S变换如果从0跑到奈奎斯特频率，每一个频点都要一次完整的搬移、加窗、逆FFT，频点数量大约是N/2个。信号长度N到几万点以后，循环几百上千次，MATLAB跑起来开始有明显延迟。这时候最有效的优化不是上并行，而是减少频点数。

工程上绝大多数场景我们只关心某一段频带。比如采样率10kHz的振动信号，你可能只关心500Hz到1000Hz的边带成分；电网暂态信号你可能只关心50Hz到2500Hz。把fmin和fmax按关注频带设窄，计算量直接按比例下降。上面的函数已经暴露这两个参数，唯一要注意的是输出矩阵的行数变少了，画时频图时y轴坐标需要自己按实际频率索引换算，不能再用默认的1:size(ST,1)。

3.3 调用与绘图：生成一个调频信号并查看结果

为了确认函数能用，先造一个频率随时间上升的chirp信号跑通全流程：

Fs = 1000; % 采样率 1000 Hz N = 2048; % 信号长度 t = (0 : N-1) / Fs; x = sin(2 * pi * (20 + 60 * t) .* t); % 频率从20Hz线性升到80Hz左右 fmin = 1; fmax = 250; % 频率索引范围，对应0.49Hz~122Hz ST = s_transform_simple(x, fmin, fmax); % 画时频图，取模值 imagesc(t, (fmin:fmax) * Fs / N, abs(ST)); set(gca, 'YDir', 'normal'); xlabel('时间 / s'); ylabel('频率 / Hz'); colorbar;

这段脚本里最容易翻车的是imagesc的y轴写法：(fmin:fmax) * Fs / N是把频率索引换算成真实Hz。很多人直接imagesc(abs(ST))然后发现图像上下颠倒，后面用set(gca, 'YDir', 'normal')修正。还有一点是MATLAB版本差异：R2021b到R2023b上这段代码都能直接跑，旧版本里如果你使用在线网页版，注意保存文件时选UTF-8编码，否则中文注释在换平台后变成乱码。导出图片时我建议直接存PNG，导出EPS在Linux上经常遇到字体缺失问题，折腾半天不值得。

4. S变换参数怎么定：频率轴、窗调节系数与输出矩阵方向

4.1 频率轴映射与奈奎斯特约束

S变换里最容易出错的是“频率索引”和“真实频率”之间的换算。结果矩阵ST的第n行，对应的是频率n·Fs/N，不是n本身。我见过不止一次有人把第50行当成50Hz去分析，实际上那可能是49Hz或53Hz，差别虽然不大，但在故障诊断里定位特征频率时这种误差会误导判断。

核心参数建议直接看这张表：

补零是个常见争议点。有人习惯把信号补零到2的幂次再算，这在某些FFT库里有加速意义，但在现代MATLAB里，长度不是2的幂也不慢。补零确实会让频率轴更细，相当于在频谱上做了插值，但它没有增加新的信息，也不会提升真实分辨率。信号长度短时，补零能让时频图看起来更平滑，这点可以用，但要清楚它的作用不是“提高分辨率”。

4.2 高斯窗调节系数p，要不要动

标准S变换的窗是固定的，但工程上会碰到两类情况：低频段窗太宽，时频图上低频成分糊成一片，时间定位差；或者高频段窗太窄，频率方向过于尖锐，带内细节丢失。这时候用广义S变换，在高斯窗指数里加一个调节系数p：

gauss = exp(−2π² · k² / n^(2p))

当p=1时就是标准S变换。p>1时，窗宽随频率增大收缩得更快，高频时间分辨率更好；p<1时，高频窗退化变慢，频率分辨率更好，类似STFT长窗的行为。要注意p的调节范围不需要很大，0.7到1.3已经覆盖了绝大多数场景。从0.5开始乱试的话，窗的形状偏离标准S变换太远，结果解释起来很别扭。

怎么判断当前p合不合适？最可靠的方法是用已知成分的合成信号做验证：造一个包含固定频率和线性调频成分的信号，跑S变换后看时频图上固定频率是否呈现一条平直线，调频成分是否呈现期望的斜线。如果低频段那条线变粗变糊，适当增大p；如果高频段出现额外毛刺，适当减小p。这个过程在第6章会展开写。

4.3 输出矩阵的方向和幅值解释

S变换的结果是复数矩阵。本文给出的函数约定“行=频率，列=时间”，这也是最常用的约定。但接手别人的代码时千万别默认这一点，先size(ST)看一眼再决定要不要转置。我见过两份代码，一份按行存频率，一份按行存时间，混用后画出的时频图横竖颠倒，排查了半天结果只是方向问题。

幅值解释是另一个容易踩的地方。abs(ST)适合做时频能量图，看能量随时间和频率的分布，但不能直接把某个点的幅值当成FFT在那个时刻的幅值。S变换窗函数的幅值本身和频率成正比，所以即使信号是恒定幅值的正弦波，高频段的时频幅值也可能看起来比低频段亮。如果只想做相对对比，可以直接看原始幅值；如果想让整个频带的能量解释更公平，可以做按频率归一化处理。相位信息同样保留在结果里，用angle(ST)可以提取，但多数时频分析场景用模值就够。

5. S变换常见问题排查：5个坑从NaN到注释乱码

坑1：结果出现大片NaN或Inf

现象：跑完S变换，结果矩阵里有大量NaN或Inf，时频图白一块或黑一块。

原因：最常见的两个来源。第一，fmin传入了0，高斯窗指数的分母n²为0，直接算出无穷大；第二，信号本身包含NaN，FFT结果跟着污染。我排查过的最离谱一例，是信号采集卡某一段数据丢了同步，波形里有几十个NaN，FFT一算全盘报废。

解决：在调用前检查信号完整性，并强制fmin至少为1。一行代码的事：

x = fillmissing(x, 'linear'); % 先填补信号里的缺失点 ST = s_transform_simple(x, max(1, fmin), fmax);

坑2：信号两端出现强烈的边缘伪影

现象：时频图的首尾两端出现异常亮的条带，而且是沿频率方向整片亮起来，中间区域正常。

原因：S变换的高斯窗在信号两端覆盖不全，窗只截到了信号的一侧，边缘处能量被低估或产生截断效应。这是S变换的固有边界问题，不是程序写错。

解决：常见做法是先对信号做镜像延拓，变换后再裁剪掉延拓部分。延拓长度取信号长度的10%左右足够：

M = round(0.1 * N); x_ext = [flipud(x(1:M)); x; flipud(x(end-M+1:end))]; ST_ext = s_transform_simple(x_ext, fmin, fmax); ST = ST_ext(:, M+1 : M+N); % 裁剪回原时间长度

坑3：时频图上下颠倒或左右颠倒

现象：图像看起来是沿水平轴翻转的，频率高的在下方，或者时间方向反了。

原因：两个来源，一个是MATLAB的imagesc默认y轴向下，另一个是矩阵约定不统一导致频率轴反向。

解决：画图时固定用set(gca, 'YDir', 'normal')。如果图像左右颠倒，检查t轴生成是从0开始还是从N开始，统一成t = (0:N-1)/Fs。这类问题最让人崩溃的是一开始没发现，等到特征提取时才发现坐标对应不上。

坑4：中文注释乱码，特别是新版MATLAB

现象：代码在Windows上打开，中文注释变成“锟斤拷”一类乱码，通常是MATLAB R2021b之后的版本默认保存编码和无注释的旧版本不一致导致。

原因：新装MATLAB在Windows上默认按GBK读取脚本文件，而脚本本身是UTF-8保存的，两边编码不匹配。这个坑在R2023b上尤其常见，网上搜“matlab 2023的中文注释乱码”能看到大量同款问题。

解决：在脚本开头加一行运行期设置，或者直接把编辑器默认编码改成UTF-8：

feature('DefaultCharacterSet', 'UTF-8');

如果乱码已经出现在文件里，用文本编辑器把文件另存为UTF-8编码格式，再重新用MATLAB打开。我自己的血泪经验是：在Linux机器上写好的脚本拷到Windows笔记本上一打开注释全乱，就是因为两边文件编码习惯不同。现在所有脚本保存前统一UTF-8。

坑5：低频段时频图糊成一片，没法看

现象：时频图高频部分很清晰，低频部分尤其是50Hz以下的能量带，看不出任何时间变化细节。

原因：这是S变换的物理特性，不是bug。低频对应宽窗，时间分辨率天然就差。你改fmin、改窗长都救不回来，因为这是公式本身的特性。

解决：能做的只有三件事。第一，改用广义S变换的p系数，稍微调整窗的缩放速度；第二，把低频段单独提取出来放大画图，至少能看个大概；第三，如果低频时间定位是硬需求，考虑换用小波变换或者直接对信号做包络分析。明白什么时候该换工具，比硬调参数更省时间。

6. 用0.5h做一个程序正确性验证：合成信号+逆变换重建

S变换写完之后，第一件事不是拿到真实信号上去跑，而是用合成信号验证程序正确。验证方法很直接：造一个成分已知的信号，跑S变换，看时频图上能不能清晰指认出已知成分，再用S变换的可逆特性重建原信号。

用下面这段脚本验证：

Fs = 1000; N = 2048; t = (0 : N-1) / Fs; x = 0.8 * sin(2 * pi * 100 * t) + 0.5 * sin(2 * pi * (30 + 50 * t) .* t); ST = s_transform_simple(x, 1, 500); % 对时间维求和还原频谱，再ifft回时域 recon = ifft(sum(ST, 2), N, 'symmetric'); err = max(abs(recon(:) - x(:))); fprintf('重建最大绝对误差：%e\n', err);

验证分成四步看：第一步，时频图上100Hz处应有一条水平亮线，从头到尾稳定；第二步，30Hz到80Hz的chirp应表现为一条斜线，斜率平稳；第三步，两条线能量应接近0.8和0.5的比例关系；第四步，重建误差应在10的负12次方量级，如果误差到了0.01以上，说明代码里有方向或窗处理问题。这四个指认全部通过，你的S变换实现才称得上可靠。

S变换真正好用的地方不止是画图。它还能当滤波器用：在时频平面上做掩膜，把目标频带和时间区域置1，无关区域置0，然后对结果沿频率轴积分再ifft，就能分离出指定时频区域的信号。掩膜边界建议做几个采样点的平滑过渡，否则会产生振铃。我自己每次拿到陌生的振动数据，第一件事永远是先跑一遍上面的合成信号验证，确认时频图里两条已知成分位置对得上，才敢把S变换的结果写进报告。不迈过这一步，后面任何特征提取都是在黑匣子上猜答案。希望帮到你。

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