尧图网站建设 
简介这是一份面向电力系统信号处理与故障分析研究的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] Σ(k0 to N−1) X[(n k) mod N] · e^(−2π²k²/n²) · e^(j2πkm/N)这个式子看着复杂拆开看就三部分第一部分X[(nk) mod N]是原始频谱搬移n个位置相当于把以n为中心的邻域频段取出来第二部分e^(−2π²k²/n²)是高斯窗在频域的形状它的宽度随n变化这就是上一节那个自适应窗的离散化身第三部分e^(j2πkm/N)是逆变换的旋转因子把频域拉回时间域。写代码时需要注意n对应的是频率索引真实频率是n·Fs/N其中Fs是采样率、N是信号长度。另一个很重要的点是n从1开始而不是从0开始因为n0时高斯窗指数项的分母为0直接算出NaN或Inf。凡是S变换一开始就出现一片NaN的人八成都是没处理直流分量这个问题。2.3 三种时频工具怎么选S变换不是万能的很多初学者把S变换理解成“更好的短时傅里叶”实际不是。它是另一种选型有自己的边界和代价。放一张对比表方便你决定当前项目该用谁工具时频分辨率频率含义逆变换计算量短时傅里叶窗长固定全频段一致直接用有成熟ISTFT算法低连续小波随尺度缩放尺度到频率需要换算有CWT基重建中等S变换随频率缩放高频时间好直接用对时间积分即可还原频谱中等偏高从实际项目选型看机械故障诊断里想提取某个特征频段的时域波形S变换逆变换方便适合电网暂态分析要同时看电压暂降起止时刻和谐波成分S变换也顺手但如果你要做的是音频实时频谱分析每个频点都要一次FFTS变换比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[(nk) 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差别虽然不大但在故障诊断里定位特征频率时这种误差会误导判断。核心参数建议直接看这张表参数建议取值影响说明fmin1或关注频段最左侧索引0频点高斯窗分母为0必须避开fmax不超过floor(N/2)1超过奈奎斯特频率对应镜像谱无物理意义NFFT长度直接用length(x)补零能细化频率轴但不增加真实信息窗调节系数p默认1.0动手范围0.7~1.3控制窗宽随频率缩放速度见4.2补零是个常见争议点。有人习惯把信号补零到2的幂次再算这在某些FFT库里有加速意义但在现代MATLAB里长度不是2的幂也不慢。补零确实会让频率轴更细相当于在频谱上做了插值但它没有增加新的信息也不会提升真实分辨率。信号长度短时补零能让时频图看起来更平滑这点可以用但要清楚它的作用不是“提高分辨率”。4.2 高斯窗调节系数p要不要动标准S变换的窗是固定的但工程上会碰到两类情况低频段窗太宽时频图上低频成分糊成一片时间定位差或者高频段窗太窄频率方向过于尖锐带内细节丢失。这时候用广义S变换在高斯窗指数里加一个调节系数pgauss exp(−2π² · k² / n^(2p))当p1时就是标准S变换。p1时窗宽随频率增大收缩得更快高频时间分辨率更好p1时高频窗退化变慢频率分辨率更好类似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直接算出无穷大第二信号本身包含NaNFFT结果跟着污染。我排查过的最离谱一例是信号采集卡某一段数据丢了同步波形里有几十个NaNFFT一算全盘报废。解决在调用前检查信号完整性并强制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-M1:end))]; ST_ext s_transform_simple(x_ext, fmin, fmax); ST ST_ext(:, M1 : MN); % 裁剪回原时间长度坑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-8feature(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变换的结果写进报告。不迈过这一步后面任何特征提取都是在黑匣子上猜答案。希望帮到你。本文还有配套的精品资源点击获取 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] Σ(k0 to N−1) X[(n k) mod N] · e^(−2π²k²/n²) · e^(j2πkm/N)这个式子看着复杂拆开看就三部分第一部分X[(nk) mod N]是原始频谱搬移n个位置相当于把以n为中心的邻域频段取出来第二部分e^(−2π²k²/n²)是高斯窗在频域的形状它的宽度随n变化这就是上一节那个自适应窗的离散化身第三部分e^(j2πkm/N)是逆变换的旋转因子把频域拉回时间域。写代码时需要注意n对应的是频率索引真实频率是n·Fs/N其中Fs是采样率、N是信号长度。另一个很重要的点是n从1开始而不是从0开始因为n0时高斯窗指数项的分母为0直接算出NaN或Inf。凡是S变换一开始就出现一片NaN的人八成都是没处理直流分量这个问题。2.3 三种时频工具怎么选S变换不是万能的很多初学者把S变换理解成“更好的短时傅里叶”实际不是。它是另一种选型有自己的边界和代价。放一张对比表方便你决定当前项目该用谁工具时频分辨率频率含义逆变换计算量短时傅里叶窗长固定全频段一致直接用有成熟ISTFT算法低连续小波随尺度缩放尺度到频率需要换算有CWT基重建中等S变换随频率缩放高频时间好直接用对时间积分即可还原频谱中等偏高从实际项目选型看机械故障诊断里想提取某个特征频段的时域波形S变换逆变换方便适合电网暂态分析要同时看电压暂降起止时刻和谐波成分S变换也顺手但如果你要做的是音频实时频谱分析每个频点都要一次FFTS变换比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[(nk) 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差别虽然不大但在故障诊断里定位特征频率时这种误差会误导判断。核心参数建议直接看这张表参数建议取值影响说明fmin1或关注频段最左侧索引0频点高斯窗分母为0必须避开fmax不超过floor(N/2)1超过奈奎斯特频率对应镜像谱无物理意义NFFT长度直接用length(x)补零能细化频率轴但不增加真实信息窗调节系数p默认1.0动手范围0.7~1.3控制窗宽随频率缩放速度见4.2补零是个常见争议点。有人习惯把信号补零到2的幂次再算这在某些FFT库里有加速意义但在现代MATLAB里长度不是2的幂也不慢。补零确实会让频率轴更细相当于在频谱上做了插值但它没有增加新的信息也不会提升真实分辨率。信号长度短时补零能让时频图看起来更平滑这点可以用但要清楚它的作用不是“提高分辨率”。4.2 高斯窗调节系数p要不要动标准S变换的窗是固定的但工程上会碰到两类情况低频段窗太宽时频图上低频成分糊成一片时间定位差或者高频段窗太窄频率方向过于尖锐带内细节丢失。这时候用广义S变换在高斯窗指数里加一个调节系数pgauss exp(−2π² · k² / n^(2p))当p1时就是标准S变换。p1时窗宽随频率增大收缩得更快高频时间分辨率更好p1时高频窗退化变慢频率分辨率更好类似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直接算出无穷大第二信号本身包含NaNFFT结果跟着污染。我排查过的最离谱一例是信号采集卡某一段数据丢了同步波形里有几十个NaNFFT一算全盘报废。解决在调用前检查信号完整性并强制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-M1:end))]; ST_ext s_transform_simple(x_ext, fmin, fmax); ST ST_ext(:, M1 : MN); % 裁剪回原时间长度坑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-8feature(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变换的结果写进报告。不迈过这一步后面任何特征提取都是在黑匣子上猜答案。希望帮到你。本文还有配套的精品资源点击获取