最近在调一维振动信号的时频分析手头能用MATLAB里现成的spectrogram但STFT的窗长一固定低频和高频的分辨率永远顾此失彼换小波变换又面临基函数选择的问题。后来我把频率切片小波变换FSWT的Matlab源码完整跑通才觉得找到了一个真正适合自己的工具它在频域里做切片分析频率可以先验设定还能直接拿系数做分量重构。这篇文章就把FSWT的Matlab源代码、时频图生成流程和调参避坑经验完整分享出来。1. 为什么是FSWTSTFT、CWT的局限和一个更灵活的频域切片思路1.1 时频分析的基本矛盾先说到底在解决什么问题。平稳信号用傅里叶变换就够了但现实里的信号基本都不平稳轴承故障时的瞬时冲击、设备升降速过程中的变频振动、语音里的音节突变频率成分是随时间变化的。这时候需要的是一张二维表——横轴时间、纵轴频率、颜色深浅代表能量大小也就是时频图。STFT的思路很直接把信号切成一段一段分别做FFT。坑在于窗长一旦确定整张图的时间分辨率和频率分辨率就都被锁死了。窗短了时间定位准但频率方向上糊成一片窗长了频率分得清但时间上什么都捕捉不到。更麻烦的是高频成分明明只需要很短的时间窗就能看清低频成分需要更长的时间窗才能分辨出频率STFT却只能用同一个窗。CWT(连续小波变换)其实是想解决这个问题的它通过伸缩母小波自动变化窗宽。但用起来有个很实际的痛点小波基函数的选择太主观。用db4、db8还是cmor3-3结果差异很大而且分析频率在尺度域上不是线性排列的低频处尺度密、高频处尺度稀想单独看某个频段不太顺手。1.2 FSWT的核心公式与直觉FSWT全称Frequency Slice Wavelet Transform中文叫频率切片小波变换。它的核心思路和STFT、CWT都不一样直接在频域里定义一个可移动的“切片窗”沿着频率轴一路扫过去每个位置做一次窄带滤波再逆变换回到时间域得到的就是该中心频率附近成分随时间的变化。公式可以写成这样FSWT(t, ω) (1/2π) ∫ ŝ(u) · φ̂*((u−ω)/σ) · e^{jut} du其中ŝ(u)是原始信号s(t)的傅里叶变换φ̂(·)是频率切片函数σ是该中心频率对应的切片宽度。注意这里的核心是σ和ω是绑定的σ ω/κ。κ称为频率切片参数。这个公式的物理含义非常直白。对每一个想分析的中心频率ω把频谱上以ω为中心、宽度σ范围内的成分“切”出来然后逆傅里叶变换回时间域。切片位置沿着频率轴移动就形成了一张完整的时频矩阵。1.3 κ参数的直观含义κ是整个算法里唯一需要人工调节的关键参数很多人第一次接触FSWT会被它搞懵但它理解起来其实比小波基简单得多。κ越大σ越小频率窗越窄频率分辨率越高能分开相邻很近的两个频率成分但时域对应的等效窗会变长时间定位能力下降瞬态冲击会被抹平。κ越小σ越大频率窗变宽时间分辨率提高但相邻频率成分会混叠在一起时频图上频率方向一团模糊。用生活化的类比κ相当于“眼睛的聚焦程度”——κ大是显微镜模式能看清精细的频率结构但视野窄κ小是广角模式能捕捉瞬态变化但频率细节模糊。后面第4节我会用一组对比实验展示这里先记住这个直觉就够了。2. FSWT的Matlab实现从频域切片到复数时频矩阵2.1 函数输入输出设计写代码之前先把接口设计清楚。我要封装一个fswt函数输入原始信号、采样频率、切片参数κ、最大分析频率fmax和频率轴采样点数M输出复数时频矩阵、频率坐标和时间坐标。function [TFR, f_axis, t_axis] fswt(x, fs, kappa, fmax, M) % FSWT 频率切片小波变换实信号版 % % 输入: % x - 一维实信号行或列向量 % fs - 采样频率(Hz) % kappa - 频率切片参数越大频率分辨率越高 % fmax - 最大分析频率(Hz)默认fs/2 % M - 频率轴采样点数默认256 % 输出: % TFR - 复数时频矩阵行对应频率列对应时间 % f_axis - 频率坐标(Hz) % t_axis - 时间坐标(s)关于输出为什么是复数矩阵这里先卖个关子后面讲单边解析滤波的时候细说。反正结果一定要保留复数相位信息不能只顾着画图就把实虚部丢了。2.2 频率切片函数如何离散化频域切片的实现核心是对FFT之后的频谱逐频点乘上一个以当前中心频率ω为中心的窗函数。我选用高斯函数作为频率切片函数原因有两个一是高斯函数在频域和时间域都是高斯形态不会产生旁瓣泄漏二是高斯窗的形状只由方差决定编程实现极其简单。实际中也有用汉宁窗、甚至自定义窗的但初次跑通建议直接用高斯。离散化的时候要特别注意频率轴的单位换算。FFT输出的频点索引对应的是fb (0:N-1)*fs/N单位Hz。对每个中心频率ω计算每个FFT频点u对应的相对偏移量η (u−ω)/σ然后代入高斯窗表达式φ̂(η) exp(−0.5·η²)这个值作为复数权重乘到频谱上。注意这里用了共轭形式conj(phi)这是为了保证时域重建时相位正确不要省掉。2.3 完整代码fswt.mfunction [TFR, f_axis, t_axis] fswt(x, fs, kappa, fmax, M) % FSWT 频率切片小波变换实信号版 % % 输入: % x - 一维实信号 % fs - 采样频率(Hz) % kappa - 频率切片参数通常0.5~20 % fmax - 最大分析频率(Hz)默认fs/2 % M - 频率轴采样点数默认256 % 输出: % TFR - 复数时频矩阵 % f_axis - 频率坐标(Hz) % t_axis - 时间坐标(s) if nargin 4 || isempty(fmax), fmax fs/2; end if nargin 5 || isempty(M), M 256; end x x(:).; N length(x); t_axis (0:N-1) / fs; % 分析频率网格从第2个点开始避开0Hz f_axis linspace(0, fmax, M); f_axis(1) []; X fft(x); fb (0:N-1) * fs / N; % FFT各个bin对应的频率 halfN floor(N/2); TFR zeros(length(f_axis), N); for k 1:length(f_axis) w f_axis(k); sigma w / kappa; % 切片宽度随中心频率变化 % 高斯频率切片函数在FFT频点上的采样值 phi exp(-0.5 * ((fb - w)/sigma).^2); % 单边加权只保留正频率部分 Xw zeros(1, N); Xw(1:halfN1) X(1:halfN1) .* conj(phi(1:halfN1)); % 除直流外正频率能量乘以2构造解析信号 if mod(N, 2) 0 Xw(2:halfN) 2 * Xw(2:halfN); % 偶数长度Nyquist点不乘2 else Xw(2:halfN1) 2 * Xw(2:halfN1); % 奇数长度没有独立Nyquist点 end % 逆FFT得到该中心频率附近分量的解析信号 TFR(k,:) ifft(Xw); end end这段代码有几个地方值得展开讲。2.4 关于单边解析滤波的设计注释第一为什么频率网格要从第2个点开始因为ω0的时候σω/κ0高斯窗分母为零直接出NaN。低频分析下限可以从fs/N、或者fmax/M这样的值开始具体取多少取决于你要分析的最低频率。第二为什么要把负频率清零、正频率加倍这是为了构造解析信号。原始实信号的FFT是共轭对称的如果直接在整个频谱上乘高斯窗逆变换后只能得到实信号幅值减半相位信息也混在一起。先把负频清零再把正频非直流部分乘2逆变换得到的就是解析信号——实部是原始带通分量虚部是希尔伯特变换结果。这样做的直接好处是abs(TFR)就是瞬时幅度包络angle(TFR)就是瞬时相位后续做解调和重构都很方便。第三注意循环里每个频率点都对整个频谱做一次加权。复杂度是O(M·N·logN)M是频率点数。M256、N1024时很快但如果信号长度到十万点M再取到512MATLAB跑起来就会比较吃力。优化思路后面第6节单独讲。3. 时频图生成从系数矩阵到可读的信息图3.1 幅度压缩为什么不能直接画abs(TFR)拿到复数时频矩阵TFR之后最直观的想法是imagesc(abs(TFR))直接画。实际画出来往往是一张几乎看不清的图——大能量成分把色标动态范围占满弱成分被压成一片深蓝。信号不同成分的能量往往差好几个数量级线性色标根本表现不了这种动态范围。正确做法是先归一化再做对数压缩P abs(TFR); P P / max(P(:)); % 归一化到[0,1] Pdb 20 * log10(P eps); % 转成dBeps避免log10(0)20·log10的转换就是dB人耳听觉和振动分析里都习惯看dB谱因为弱信号在dB标度下也能看出结构。之后绘图时把色标下限设在−60或−50dB上限0dB低于下限的值都显示为背景色这样动态范围不足的小成分依然能显示。3.2 测试信号构造与FSWT时频图脚本下面用一个人造信号完整演示时频图生成。信号包含三个成分一个线性扫频chirp、一个固定频率正弦波、一个短时高斯脉冲。这三个成分分别代表了连续变频、稳态窄带和瞬态冲击三类典型特征。% demo_fswt.m clc; clear; close all; fs 1024; N 1024; t (0:N-1) / fs; % 三成分测试信号 s1 sin(2*pi*(50*t 60*t.^2)); % 线性chirp频率从50Hz扫到约170Hz s2 0.6 * sin(2*pi*220*t); % 220Hz固定正弦 s3 0.8 * exp(-((t-0.65)/0.008).^2) .* sin(2*pi*380*t); % 380Hz短时冲击 x s1 s2 s3; % FSWT参数 kappa 2; fmax 500; M 256; [TFR, f_axis, t_axis] fswt(x, fs, kappa, fmax, M); % dB压缩 P abs(TFR); P P / max(P(:)); Pdb 20 * log10(P eps); % 绘图 figure; imagesc(t_axis, f_axis, Pdb); axis xy; clim([-60 0]); % 老版本MATLAB用 caxis([-60 0]) xlabel(时间 (s)); ylabel(频率 (Hz)); title(FSWT时频图 (\kappa2)); colorbar; colormap(jet);如果环境是R2022a之后的MATLABclim直接可用再老一些的版本需要换成caxis。建议统一写成caxis([-60 0])兼容性更好我在自己的机器上两个都试过效果一致。3.3 读图三条能量脊对应三种成分跑完上面的代码你会在时频图上看到三条清晰的能量脊一条从左下角向右上方倾斜的亮带是chirp信号起始频率约50Hz随时间线性上升一条大约在220Hz处水平的亮线是正弦分量能量稳定贯穿整个时间轴一条在大约0.65s、380Hz处的竖直亮斑是瞬态冲击时间定位非常准。这就是FSWT时频图最典型的输出形态。对比STFT你会明显感觉chirp那条线更细、更连续而冲击成分的起止时刻也更接近真实位置。原因在于FSWT在每个中心频率都按σω/κ自适应地调整了切片窗高频处窗更宽、时域定位更好低频处窗更窄、频率分辨率更高正好踩中了“高频用短窗、低频用长窗”的理想模式。4. κ参数怎么定用一组对比实验理解时频分辨率博弈4.1 κ与切片带宽的定量关系在动手调参之前先建立定量的概念。高斯窗在频域的−3dB带宽大致是Δf≈2.355·σ代入σω/κ得到Δf ≈ 2.355·ω/κ也就是说中心频率越高绝对带宽越大但相对带宽Δf/ω≈2.355/κ是常数。这正好对应“恒Q滤波”特性也是FSWT区别于STFT的重要特征。如果在100Hz处想要5Hz的分辨率κ≈47在1000Hz处想要50Hz的分辨率κ≈47。所以κ本质上是一个相对频率分辨率的控制旋钮和你分析的具体频段无关。时域对应的等效窗长约等于带宽的倒数。Δf越大时域窗越短时间定位越好。想提高时间分辨率就往小调κ想提高频率分辨率就往大调κ两个方向永远互相压制没有绝对最优值。4.2 四组κ值对比实验用第3节的混合信号做同一个FSWT分别取κ0.5、2、8、20画成2×2的子图对比kappas [0.5, 2, 8, 20]; figure; for k 1:4 [TFRk, fk, tk] fswt(x, fs, kappas(k), fmax, M); Pk 20*log10(abs(TFRk)/max(abs(TFRk(:))) eps); subplot(2,2,k); imagesc(tk, fk, Pk); axis xy; caxis([-50, 0]); title([\kappa , num2str(kappas(k))]); xlabel(时间 (s)); ylabel(频率 (Hz)); end colormap(jet);实际跑完你会看到非常鲜明的差异κ0.5频率方向几乎糊成一片chirp和220Hz正弦分不清边界但0.65s处那个冲击的起止时刻非常锐利时间定位一流κ2三条成分都能看清chirp线和220Hz线可以分辨冲击的时间宽度稍微变肥了一点整体最均衡κ8频率分辨率显著提升chirp线很窄很锐但冲击变成了一小片“蝌蚪”时间上拉长了不少κ20频率方向极其锐利两个正弦成分分得清清楚楚但瞬态冲击在时间方向上被拉得很宽甚至和chirp能量脊尾部混在一起。如果你关心的对象是冲击类故障κ0.5~2比较合适如果你是做多分量稳态信号的频率分离κ可以拉到8以上混合信号里既有瞬态又有窄带2~5是比较稳妥的区间。4.3 选κ的实战经验我自己的习惯是先用大κ比如10快速扫一眼频率上有哪些峰值确定主要频率成分的位置再针对关心的频段用小κ重算观察瞬态的时间定位。两个方向各算一遍基本上就清楚了。如果算力充足还可以在κ1~10里做一次小网格扫描看时频图哪个值下“最清爽”——能量脊最连续、没有明显染色和模糊就选哪个。另外提醒一点κ不要取非整数太小的值比如κ0.2此时每个切片窗都非常宽多个成分混在一个窗里时频图会出现互调伪影看起来像多出了一些交叉亮带但那不是真实的信号成分。5. 选择性重构把混合信号中的单个分量“抠”出来5.1 为什么FSWT适合选频重构FSWT比STFT和CWT更适合做选频重构原因是频率切片函数在频域有很好的局部性。对某个中心频率ω切片之后系数序列本质上就是该频率附近成分的窄带解析信号。那么把一段频率范围内的所有切片系数叠加起来就能从混合信号中还原出该频段的分量。这个思路在工程上非常实用齿轮箱故障诊断里要提取啮合频率及其边带轴承故障里要提取共振频段的包络信号语音处理里要把某个频段的共振峰单独拿出来分析都可以用FSWT的选频重构完成。5.2 重构步骤与Matlab代码还是用第3节的三成分信号目标是把220Hz正弦分量从混合信号里单独提出来。步骤分三步先算FSWT再按频率范围选择系数行最后对选定行求和取实部。% 提取220Hz正弦分量 idx find(f_axis 190 f_axis 250); x220 real(sum(TFR(idx, :), 1)); figure; subplot(3,1,1); plot(t, x, k); title(混合信号); xlabel(时间 (s)); subplot(3,1,2); plot(t, s2, b); hold on; plot(t, x220, r); legend(原始220Hz分量, FSWT重构结果); title(220Hz分量对比); xlabel(时间 (s)); subplot(3,1,3); plot(t, x220 - s2, g); title(重构误差); xlabel(时间 (s));跑完你会看到重构波形和原始正弦分量基本重合误差主要出现在信号两端和瞬态冲击附近。频率范围选得越准重构效果越好。5.3 重构精度的影响因素几个直接影响重构精度的细节必须留意。第一频率选择边界不要太硬。直接find截段相当于在频域乘了一个矩形窗会在重构信号两端引起吉布斯振铃波形边缘出现毛刺。更稳的做法是给选中的系数行加一个平滑的权向量比如高斯权重或余弦渐变让边界过渡柔和一些。我实际测试过加了20%余弦过渡之后误差峰值能降一半以上。第二系数累加的幅值会受到高斯窗加权的影响。每个切片系数在做频域加权时乘了高斯窗累加之后模拟出的是经过高斯滤波的分量幅值会有微微衰减。如果对幅值精度有要求可以用最小二乘或峰峰值对齐做一次校正。简单起见用x220 x220 * (rms(s2)/rms(x220))这类定标方式也能对齐。第三重构前最好先对原始信号做去均值处理。直流分量在ω0附近会污染低频切片系数如果你选的重构频带恰好很低直流的泄漏会明显加大误差。我的习惯是在进入fswt之前先执行x x - mean(x)。6. 实战中踩过的坑与排查建议6.1 六个高频问题汇总写这个函数的过程中我在各种信号上测试过也遇到过不少奇怪现象。把最常踩的坑集中列在下面。现象根因解决办法时频图顶端出现倒影亮带fmax取得太接近fs/2高频切片窗泄漏到负频镜像fmax不要超过0.9*fs/2留出过渡带低频段出现NaN或空白ω0时σ0除零频率网格从fmax/M或fs/N开始避开0Hz时频图动态范围差、弱成分几乎看不见直接画abs(TFR)没有做dB压缩归一化后转dB色标下限设为−50~−60dB重构波形两端大幅振荡频域选择加了矩形硬边界用平滑权向量或余弦过渡κ调大后瞬态冲击被拉成宽条大κ导致频窗窄、时窗长时间定位下降瞬态分析用κ0.5~2频率分离用大κ信号很长时计算非常慢频率点数M和信号长度N都很大会让循环变得极慢M控制在256以内必要时先降采样/滤波6.2 复信号、短信号与计算量的边界条件上面给出的fswt函数是为实信号写的。如果输入是复信号比如通信里的基带信号单边清零的逻辑就不成立了——复信号的频谱本来就不对称负频也承载着独立信息。这种场景需要改成双侧加权版本直接Xw X .* conj(phi)然后TFR(k,:) ifft(Xw)不要清零也不要加倍。短信号也要特别注意。N128、fs1024时FFT频率分辨率只有8Hz而低中心频率对应的σω/κ可能还不到几个Hz。这时高斯窗宽小于一个FFT bin切片操作几乎选不到有效能量时频图上会出现断断续续的“马赛克”。解决办法有两个一是给信号补零到足够长度再做FFT二是适当增大κ或调整频率网格让σ不小于2~3个频率分辨率。计算量方面最慢的部分是M个频率点各做一次N点IFFT。N1024时完全无压力N4096时稍微等两三秒N16384以上、M还是256的话建议先降采样或者把信号分段处理。另一个思路是矩阵化预先构造一个M×N的切片权矩阵用矩阵乘法一次性出结果代码稍微复杂但速度提升明显感兴趣的可以自己试着优化。6.3 几个建议如果你准备把这套代码用到自己的研究或项目里我的建议是先把第4节的κ对比实验跑一遍用自己的信号去找“手感”不要直接照抄别人论文里的参数。然后保留TFR的复数输出不要只存abs(TFR)——后续做相位分析、重构、瞬时频率估计都要用到复数信息。最后单独维护一个参数记录表把每次分析的采样率、fmax、M、κ、预处理方式记下来等到需要复现结果的时候你会感谢当时的这个习惯。我后来在轴承内圈故障数据上反复试kappa固定在3~5、fmax设在4000配合选频重构出的包络谱基本一次能打中故障特征频率。如果你也在做类似的时频分析建议先把上面的对比实验完整跑一遍理解每个参数在图上留下的“指纹”再上自己的数据会少走很多弯路。