S变换Matlab实现：从公式到代码的时频分析指南-北京尧图网络

简介基于Stockwell变换理论的MATLAB源码包内含1个.m脚本适合电力信号处理、科研复现与教学演示。S变换在短时傅里叶变换基础上引入尺度参数使时频分辨率可随频率自适应调节能有效解析电压暂降这类非平稳瞬态信号。借助该脚本用户可提取基频幅值、相位跳变PHASEJUMP、突变点及频率幅值包络线等关键特征据此判断暂降深度、起止时刻与扰动方向还能观察谐波频率成分与幅度分布辅助识别非线性负载造成的谐波污染。当电网出现负荷波动或设备故障时脚本生成的时频图像可直观展示暂降伴随的频率波动趋势帮助精确锁定异常区间并理解电网动态行为。资源包为RAR格式共1个文件、约2KB轻量易用下载后即可在MATLAB中运行便于快速复现算法或嵌入现有分析流程已有629人学习/下载对电力系统电压暂降检测、谐波源定位与后续故障诊断均有直接实用价值。1. S变换是什么把一维信号摊成时间-频率平面做非平稳信号分析的人多半有过这种经历一段振动信号里既有缓慢变化的基频又有突然出现的冲击成分用短时傅里叶变换STFT去切窗窗口短了低频模糊窗口长了高频模糊怎么调都不对劲。S变换S-TransformStockwell变换用一个随频率自动伸缩的高斯窗把这个问题换了个解法低频用宽窗换频率分辨率高频用窄窗换时间分辨率。它的核心公式不复杂Matlab实现也就几十行但能把一维信号变成二维的“时间-频率-幅度”谱在旋转机械故障诊断、电力暂态分析、地震信号处理和脑电信号分析里都很好用。这篇文章就从连续公式讲到离散实现再给出一套能直接跑的Matlab代码和踩坑记录。2. S变换公式不玄从连续定义到离散Matlab实现的一步路要把S变换用起来公式这部分绕不过去。好在这个公式比起小波变换的尺度函数理解成本低不少Matlab里走的也是FFT的老路子。2.1 为什么选S变换而不是STFT或小波三句选型话第一句给STFT它的窗长是固定的时间分辨率和频率分辨率在整张时频谱上被捆绑死。分析频率跨度大的信号高频细节和低频轮廓总要牺牲一个这是窗函数的不确定性原理决定的调参只能缓解不能根除。第二句给小波小波的自适应多分辨率很漂亮但尺度轴和实际频率之间的映射不直观画出来的时频谱要换算半天更麻烦的是小波系数里的相位和原始信号的傅里叶相位关系很绕想做滤波重构还得折腾一番。第三句给S变换它的高斯窗宽度跟着频率走低频窗宽、高频窗窄而且变换结果可以直接逆变换回原信号。S变换的每个复数点都保留着明确的绝对相位信息适合做时频滤波、信号重构这类后续操作。我做旋转机械故障诊断时常用S变换看变频启动过程的时频谱STFT和小波都没它省事。2.2 离散S变换公式拆解三个部分各管什么事连续S变换的定义式长这样S(τ, f) ∫ x(t) · (|f| / √(2π)) · e^(-f²(τ-t)²/2) · e^(-i2πft) dt这里面三个部分分工很清楚x(t)是待分析信号中间那个高斯窗系数保证窗对时间的积分为1最后的复指数负责频率扫描。高斯窗的中心在τ宽度是1/|f|所以频率越高窗越窄时间定位越准频率越低窗越宽频率分辨能力越强。窗面积恒为1这一点很关键它保证了对τ积分后高斯窗自动消失剩下正好是信号的傅里叶谱X(f)。直接按连续公式做数值积分太笨常见做法是利用傅里叶变换把S变换改写成频域形式再离散化S[k, n] Σ_m X[mn] · e^(-2π²m²/n²) · e^(i2πmk/N)这里X是信号x的FFT结果n对应第n个频点m是相对于频点n的循环偏移量e^(-2π²m²/n²)就是频域里的高斯窗最后对m做一次逆FFT就回到时间域。这个形式的好处是彻底避开了对时间的积分只需要FFT、逐点乘法、逆FFT全是Matlab的强项。2.3 最小可运行的S变换代码从零开始写一个不花哨的版本下面这个函数是最朴素的实现不优化、不花哨先保证算得对。输入要求是列向量长度记为N结果是一个N×N的复数矩阵行是时间列是频率。function ST st_basic(x) % st_basic 最小可运行的S变换实现 % 输入 xN×1 列向量信号长度建议为偶数 % 输出 STN×N 复数矩阵行时间列频率 N length(x); X fft(x); % 信号FFT未归一化 ST zeros(N, N); % 预分配结果矩阵 for n 0:N-1 % n频点索引从0直流到N-1 if n 0 % 零频列S变换对时间积分为频谱直流列填均值 ST(:, 1) mean(x); else % 频域高斯窗偏移量m取 -N/2 到 N/2-1 m (-N/2:N/2-1); win exp(-2 * pi^2 * m.^2 / n^2); % 循环平移索引取X中nm位置超出部分绕回 idx mod(n m, N) 1; % 逐点相乘后逆FFT回时间域 ST(:, n1) ifft(X(idx) .* win); end end end调用方式很直接fs 1000; % 采样率 1000 Hz t (0:1023) / fs; % 1 秒数据 x sin(2*pi*50*t) 0.5*sin(2*pi*180*t); % 两个频率叠加 ST st_basic(x(:)); % 强制列向量 size(ST) % 结果 1024×1024这里有几个参数上的讲究。频点n从0到N-1第0列是直流分量单独填均值因为n出现在高斯窗的分母上n0时必须避开不然就是除零。窗宽在高斯窗里由n的平方决定没有额外参数这是标准S变换的定义后续要加调节旋钮我们放到第4章讲。逆FFT用Matlab默认的1/N归一化这个约定换来了一个非常好的性质ST矩阵每一列求和正好等于信号的FFT这为最后一章的验证方法埋了伏笔也是很多人实现S变换时栽跟头的地方。3. 把S变换跑快向量化、矩阵化和限频域三种优化最小版本能跑但别急着上大信号。N1024时双重循环还能等N8192时那个嵌套循环能把整个下午搭进去。这一章讲清楚慢在哪以及常见的三档加速手段。3.1 先看复杂度O(N²)的双重循环在Matlab里怎么拖死你看代码结构就知道外层循环n走N次内层每次做一次长度为N的ifft整体复杂度是O(N² logN)空间上要保存N×N的结果这也是S变换逃不掉的下限。STFT是O(N logN)级别两者差出一个量级所以S变换从来不是处理超长实时信号的方案。N1024时外层循环1024次每次ifft 1024点Matlab跑起来大约一两秒还能接受。N4096时就是16倍的计算量外加128MB的复数矩阵内存开销。我在实际项目里处理过一段60秒的振动信号采样率5120HzN307200直接全套计算完全不现实。这时候要做的第一件事是停下来想想你真的需要全频段全时间的S变换吗绝大多数诊断场景需要的只是某一频带内的时频细节。如果你只是要验证算法正确性N512到2048足够了循环版跑得动结果也好观察。3.2 矩阵化实现用一次ifft算所有频率最小版本的慢主要慢在Matlab的for循环。改进思路很直接既然每个频点n的计算互不依赖可以把所有频点的高斯窗和索引一次性构造成矩阵然后让ifft沿着矩阵维度一次算完。代码量反而更短function ST st_mat(x) % st_mat 矩阵化S变换适合N4096的短信号 N length(x); X fft(x); m (-N/2:N/2-1); % 频率偏移向量N×1 n 0:N-1; % 频点索引1×N [Mn, Nn] ndgrid(m, n); % 生成 N×N 网格 % 高斯窗n0列用1代替避免除零后面会覆盖 safe_n Nn; safe_n(Nn 0) 1; win exp(-2 * pi^2 * Mn.^2 ./ safe_n.^2); % 循环平移后的FFT索引逐点取X idx mod(Nn Mn, N) 1; V X(idx) .* win; % 沿第1维时间做逆FFT一次处理所有频率列 ST ifft(V, [], 1); % 覆盖零频列 ST(:, 1) mean(x); end这个版本的逻辑和最小版完全一致只是把“每列算一遍”换成了“整张矩阵一起算”。ndgrid生成的两张网格里Mn是每个位置的频率偏移量Nn是频点编号两者相除正好就是离散公式里的m/n。idx的处理是精髓mod(nm, N)1把负频率绕回FFT数组的后半段这一步少了高频部分会错得莫名其妙。ifft的第二个参数[]表示不改变变换长度第三个参数1表示沿列方向做得到的结果每一行还是时间轴。矩阵化适合N不超过4096的信号。N4096时win、V、ST三个复数矩阵各占128MB加起来近400MB稍不注意就触发内存不足。所以这个版本要在内存和速度之间做取舍N超过4096回头用循环版逐列算反而更稳。3.3 更挑剔的方案只算关心的频段或分块降采样实际工程里我一般不会全频段跑。以电力暂态分析为例关心的谐波和间谐波集中在0到500Hz采样率却是6400HzN可能上万。这时候把频点n限定到只算关心的范围计算量直接砍到原来的几分之一。function ST st_band(x, n_max) % st_band 限频段的S变换 % n_max需要计算的最高频点索引对应频率 n_max*fs/N N length(x); X fft(x); n_list 0:n_max; ST zeros(N, length(n_list)); for k 1:length(n_list) n n_list(k); if n 0 ST(:, k) mean(x); else m (-N/2:N/2-1); win exp(-2 * pi^2 * m.^2 / n^2); idx mod(n m, N) 1; ST(:, k) ifft(X(idx) .* win); end end end这样算出来的ST不是方阵列数由n_max决定。代价是丢失了高频信息逆变换验证时只能对部分频带成立。另一个常用做法是先把原始信号降采样或分帧每帧N取1024或2048帧间重叠50%逐帧算S变换再拼接。分帧会引入边界效应第5章会专门讲怎么处理。三种提速手段不是互斥的。我的习惯是新信号先降采样到合适范围再限频段N能压到2048以内就上矩阵化压不下来就用循环版逐列算。N超过一万且只关心低频st_band的循环版远优于直接st_mat因为省掉的列越多节省的内存越可观。4. S变换的3个必调参数窗宽因子、频率轴与可视化输出公式和代码跑通后下一个问题是怎么看结果。S变换虽然自称“自适应”但实际使用中可调的参数至少有三个高斯窗的广义宽度、频率轴的选取范围、以及复数结果的可视化策略。这三个调不好时频谱画出来不是一片糊就是一片噪根本读不出信息。4.1 高斯窗宽度γS变换唯一的分辨率旋钮标准S变换的高斯窗宽度由频率f完全决定窗宽σ1/|f|没有自由参数。这是优点也是缺点有些信号高频段太吵时间分辨率过高反而看不出趋势。工程里常见的做法是引入一个宽度因子γ变成广义S变换窗宽变成1/(γ|f|)。γ大于1时窗更窄时间分辨率更好γ小于1时窗更宽频率分辨率更好。离散公式里高斯窗从exp(-2π²m²/n²)变成exp(-2π²m²/(γn)²)。代码改动很小gamma 1.2; % γ1 时间分辨率更好 win exp(-2 * pi^2 * m.^2 / (gamma * n)^2);γ取值对结果的影响可以粗略分三档γ取值窗宽时间分辨率频率分辨率适用场景γ 1变宽变差变好慢变信号、谐波分析γ 1标准标准标准默认首选γ 1变窄变好变差突变检测、冲击定位我的经验是γ不要一上来就调。第一次分析用γ1跑一遍看出明显的“顾此失彼”再动它。比如低频谐波糊成一片分不清几根线就往小调比如0.5到0.8高频冲击定位不准就往大调比如1.5到2。γ变化对结果的影响是连续的不存在“最优值”这种说法只能按任务调。4.2 频率轴到底怎么摆物理频率、FFT bin与负频率S变换算出的ST矩阵第n列对应的是FFT的第n个bin换算成物理频率是f_n n·fs/N。n从0到N/2对应0到fs/2这是正频率段n超过N/2的部分对应负频率。画图时绝大多数场景只画正频率段也就是ST的前N/21列。我见过很多人第一步就是画全频段结果图上一大半是负频率的镜像还纳闷为什么能量上下对称。正频率轴这样构造N size(ST, 1); f_axis (0:N/2) * fs / N; % 从0到fs/2的频率轴 ST_pos ST(:, 1:N/21); % 截取正频率部分实际应用中频率轴还可以进一步裁剪。机械故障诊断关心轴承外圈故障频率附近几十赫兹的边带电力谐波分析关心基波和整数次谐波而这些只需要在f_axis上做逻辑索引比如mask f_axis 40 f_axis 60; ST_band ST_pos(:, mask); f_band f_axis(mask);这样做一方面让后续分析聚焦另一方面内存占用跟着降下来。ST_pos本身仍然是N×(N/21)的复数矩阵单精度时N2048大约占用16MB不算大可以放心截取。4.3 把复数结果画成热力图幅值谱、相位谱的取舍S变换输出是复数矩阵直接画图必须取模或者取相位。幅值谱是最常用的反映的是某个时刻某个频率上成分的强弱。画幅值谱我通常会配合对数压缩或clim截断因为S变换的幅值动态范围很大直流和强谐波会把其他弱成分压制得完全看不见。figure; imagesc(t_axis, f_axis, abs(ST_pos).); set(gca, YDir, normal); % 让频率从低到高显示 xlabel(时间 (s)); ylabel(频率 (Hz)); colorbar; title(S变换幅值谱); clim([0, prctile(abs(ST_pos(:)), 95)]); % 截掉最强5%能量突出弱成分关键点有两个。第一imagesc的第一个参数是x轴、第二个是y轴数据矩阵必须是y轴×x轴所以ST_pos的转置不能少。第二imagesc默认把矩阵第一行放在图像底部加了YDir, normal后低频才会在下面不然整张图上下颠倒分析方向全错。相位谱很少直接画因为低幅值区域相位噪声很大意义不大。如果确实要看相位一个常见做法是先按幅值阈值做掩膜只显示强成分的相位比如只显示幅值超过最大值10%的点。最后提醒一句可视化是给人看的图的参数调整只影响读图观感不改S变换本身的数值。真正影响数值的只有γ、频点范围和输入信号质量别为了图好看反复调clim结果分析却用原始数据两套东西对不上这是我自己走过的弯路。5. S变换Matlab实现避坑实录5条我从代码和信号里踩出来的经验这一章全部来自实际调试过程。每一条我都按“现象→原因→解决”的顺序写很多坑不是S变换本身的问题而是实现细节和Matlab特性叠加出来的但表现形式都像算法写错了。5.1 现象逆变换重构误差大信号整体放大N倍写完第一版S变换后我习惯立刻做重构验证。验证方法是把ST矩阵每列求和再逆FFT也就是x_rec ifft(sum(ST, 1))结果跑出来重构信号的幅值是原信号的几百上千倍波形形状倒是对的。原因出在ifft的归一化上。Matlab的ifft自带1/N归一化也就是说ifft(fft(x))恰好等于x。但很多论文里的离散S变换公式默认使用IFFT时再做一次N倍缩放两种约定混着写我在ifft(tmp)后面额外乘了N导致每列都放大N倍重构自然就放大N倍。解决方法是删掉额外的乘N让ifft保持默认归一化。这样ST矩阵的每一列求和自动等于fft(x)的对应值重构验证能到10的负12次方量级的误差。以后再看到别的实现里乘了N先确认他的逆变换公式是不是也做了对应补偿不要直接照抄。5.2 现象低频段一片竖直条纹高频段全是噪点时频谱画出来后低频部分像斑马线一样一道一道的竖条纹高频部分则散布大量孤立亮点整张图看起来脏兮兮的。很多人的第一反应是代码bug查半天没查出来。原因分两层。低频条纹是S变换本身的特性低频对应宽窗时间分辨率差周围时间点的能量被平均到一大片区域视觉上相邻时间列数值平滑于是形成条纹。高频噪点则是窄窗带来的时间过采样每个时间点几乎独立信号里微小的噪声都被如实放大。解决方法是先做信号预处理去掉趋势项和直流偏置再用对数色标或者按百分比截断动态范围来画图。另外可以先对信号做带通滤波把不关心的频带清掉再算S变换高频段的噪声会干净很多。5.3 现象信号首尾出现异常高能量边界像被“撞”了一下算出来的时频谱时间轴两端总是有一列明显的强能量而且频率越高越明显看起来像信号首尾突然多了一坨东西。这个问题在分帧处理时尤其严重。原因是FFT隐式的周期延拓。S变换的离散公式里idx mod(nm, N)1会自动把超出范围的索引绕回信号开头等于默认信号首尾相连。真实信号首尾一般不是连续过渡的这个人为的跳变会被FFT当成高频成分放大边界处就出现虚假能量。解决方法是计算前对信号做去均值加Tukey窗预处理窗的锥度部分占5%到10%就够或者算完S变换后把开头和结尾各丢弃N×5%的列只读中间区域。分帧处理时帧与帧之间必须有重叠一般50%以上能明显缓解边界效应。还有一个血泪经验验证重构时不要拿加窗后的信号去和原始信号比误差窗函数会改变信号本身重构验证必须用没加窗的原始信号或者反变换后再除以窗函数。5.4 现象N取2048没事N取8192直接内存不足st_mat版本在N2048时运行流畅换成N8192就报Out of Memory。这是S变换的O(N²)空间复杂度撞上了Matlab的复数存储机制。原因在于Matlab默认复数double占16字节st_mat里win、V、ST三个8192×8192矩阵加起来超过3GB。就算机器内存够Matlab也可能因为变量碎片化分配失败。解决方法是分层处理。短信号用st_mat长信号改用st_band循环版只保留需要的频段或者把变量改成单精度复数也就是在乘法前后用single转换内存直接减半虽然精度损失对时频分析基本可忽略但对重构验证有微小影响。另一个方向是把长信号分帧每帧2048点帧间50%重叠算完再拼时间轴这样每帧内存不超过100MB工程上最稳。5.5 现象输入是行向量还是列向量直接决定函数会不会翻车把信号从工作区读进来直接喂给st_basic结果要么报错说维度不匹配要么ST矩阵大小变成1×N而不是N×N。问题出在x是行向量时mean(x)是标量还好但X(idx)是1×N的行向量win也是1×N乘以ifft没问题可初始化zeros(N,N)时把1×N的结果赋给ST(:,n1)是N×1的一列Matlab会直接报错。解决方法是所有函数开头强制统一成列向量一行代码的事x x(:); % 无论行向量还是列向量都转为N×1另一个连带问题是验证时sum(ST, 1)得到的行向量方向与fft(x)的方向可能不一致对比时要用transpose配平。养成在一开始的脚本里就统一方向的习惯后面所有验证代码都不会因为维度问题分心。这不算S变换的数学坑但亲手把函数接口设计得健壮一点能省掉很多调试时间。6. 用“列求和ifft”验证你的S变换五分钟自检一个实现S变换有一个非常好的自检性质利用它能在不依赖任何第三方工具箱的情况下确认你写的函数没算错。原理就是前面提到的窗函数积分为1对S变换的时间轴求和高斯窗被积掉剩下的正好是信号的傅里叶谱。写成验证代码只有三行X_check sum(ST, 1); % ST每列求和N×1 x_rec ifft(X_check); % 逆FFT回时域 err max(abs(x_rec - x)) / max(abs(x)); % 相对误差如果实现正确ST的每一列分别求和后得到的向量应当等于fft(x)。而ifft(fft(x))必然等于x所以只要err在1e-12量级说明从FFT、高斯窗、循环平移到ifft归一化这一整条链路都是自洽的。err如果掉到1e-3量级通常不是浮点误差而是哪个地方把归一化或者索引搞错了回到第5章的坑里逐一排查。这个验证方法还有两个用法。第一检查ST的某一列sum(ST(:, n1))应当约等于fft(x)的第n1个值单独验证某一列可以帮你快速定位问题是出在所有频率列还是特定频率上。第二如果你实现了限频段的st_band版本验证要改成只对比被保留的频点比如X_check zeros(1, N)再把限频段那些位置的X_check对应值填入sum(ST, 1)的结果其余位置补零然后再ifft。这样得到的重构信号只包含关心频带的成分时域上相当于一次滤波误差评估标准也要换成带内容差而不是整体误差。我现在的习惯是每写一个新版本或调一次γ先跑一遍这个五分钟自检确认数值链路没问题再进入参数调优和文章画图环节。S变换这类带“自适应”标签的算法最容易让人困惑的就是不知道哪里算对哪里算错这个自检方法能帮你把置信区间锁死。如果你刚接触S变换建议从最小版本跑通验证开始再逐步换成矩阵化、限频段和广义窗。希望帮到你。本文还有配套的精品资源点击获取

S变换Matlab实现：从公式到代码的时频分析指南

相关推荐

Linux mkdir命令详解：从基础使用到高级技巧

ePass3000GM SDK V1.4 Windows集成避坑指南

LYT-NET：轻量化Transformer在低光照图像增强中的应用

更多资讯

钢化玻璃隔断：办公室与家装空间分区优选

Swift-1.5-Qwen3.8-27B-GGUF 快速上手：如何用 llama.cpp 三步跑起来本地部署

django餐饮公司管理系统12151-计算机课程设计、毕业设计

DLMS/COSEM 蓝皮书解读（五十一）：Measurement data monitoring objects（class_id = 66）—— 抓一段“事故波形“

进程 — 解析计算机运行的精妙

VS Code 差异对比实战指南：4 个对比入口与 4 个窗格搞定 Git 冲突

今日更新

AI智能体开发实战：从技术原理到商业落地

OpenShell 完全定制指南：从安装到深度配置，打造高效 Windows 开始菜单

Superpowers技能体系全解析：从Skills引入到安装的完整指南

本周精选

防幻觉智能体架构：电源硬件设计中的大模型编排与工具化实践

OpenGL面试核心考点：渲染管线、矩阵变换与性能排查全解析

全国计算机二级C语言真题PDF高效刷题指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？