简介:一份面向信号处理、故障诊断与工业智能运维方向研发人员的MATLAB项目实例,围绕短时傅里叶变换与支持向量机结合,解决非平稳振动信号下的多类别故障识别问题。内容覆盖原始振动信号采集与预处理、STFT时频图构造、关键频带能量特征提取、SVM模型训练与交叉验证、混淆矩阵评估以及单样本时频图与分类结果联动展示,并给出完整程序、GUI设计和代码详解,便于复现与扩展至旋转机械、电力系统、轨道交通等状态监测场景。压缩包内共1个docx文件,约120KB,集中承载项目说明、模型架构、参数配置与工程部署建议,目录按数据准备、时频分析、样本集划分、模型训练和可视化等模块组织,方便按流程查阅。已有132人学习,适合具备一定MATLAB基础、希望深入理解时频分析与机器学习结合应用,并关注特征工程与故障机理关联的工程师和科研人员参考。
1. 振动信号里看不出故障时,STFT-SVM 该怎么切入
电机驱动端轴承内圈出现早期剥落时,加速度传感器采回的时域波形幅值往往只比正常状态高一点点,包络谱上的边频带被背景噪声盖住。靠肉眼看波形,或者只截一段做 FFT 看幅值谱,很容易判成正常。真正能区分状态的证据是:冲击在时间轴上周期性出现,同时能量集中在某几个高频共振带——这是时域和纯频域都说不清楚的第三类信息。
短时傅里叶变换(STFT)把一维信号切成许多短帧、逐帧做 FFT,把「什么时候」和「哪个频率」同时保留下来;支持向量机(SVM)再从这张时频矩阵里学出一个能把正常、内圈、外圈、滚动体分开的边界。整套流程在 MATLAB 里从 matlab安装步骤走完就能跑通,数据处理、训练、画图、GUI 全在一个环境内闭环,既适合做轴承故障诊断的入门项目,也适合把已有的故障诊断代码整理成能交付给同事的小工具。
2. STFT 时频特征提取:从原始振动信号到 SVM 能吃的特征矩阵
SVM 的输入必须是一个定长的数值向量,而 STFT 输出的是「频率 × 时间」的二维复数矩阵。这一步的核心工作,就是把二维时频矩阵压成一个既能保留故障信息、又不会让维度爆炸的特征向量,同时保证正常样本和故障样本走的是完全相同的处理链路。
2.1 窗长、重叠率与 FFT 点数:三个参数决定特征质量
STFT 的本质是「加窗截断 + 逐帧 FFT」。窗长决定了你在时间轴上看得多细,FFT 点数决定了频率栅格有多密,重叠率决定了帧与帧之间会不会漏掉冲击。
以常见的轴承数据为例,采样率取 12 kHz。窗长 256 点时,单帧覆盖约 21.3 ms;FFT 点数补零到 512,频率分辨率为 12000/512 ≈ 23.44 Hz,足以分辨出轴承外圈故障特征频率附近的调制边带。重叠取窗长的 75%,即每次只前进 64 点,保证一次短促冲击至少落在两帧里,不会因为分帧位置恰好错过而被抹平。
下表是这套参数在不同场景下的取舍关系,可以按自己的数据替换:
| 参数 | 取值 | 影响 | 调大后的代价 |
|---|---|---|---|
| 窗长 | 128 / 256 / 512 | 时间分辨率与频率分辨率的平衡 | 窗越长,频率越清晰,冲击定位越模糊 |
| 重叠率 | 50% / 75% / 90% | 短时冲击的捕获率 | 帧数线性增长,训练样本与耗时同步上升 |
| FFT 点数 | 512 / 1024 | 频率栅格密度 | 特征维度翻倍,小样本下容易过拟合 |
| 窗函数 | hann / hamming | 抑制频谱泄漏 | 主瓣变宽,相邻谱线区分度下降 |
注意:如果用的是 R2019b 之前的 MATLAB,函数名要换成
spectrogram,参数写法不同,但含义一致。
2.2 MATLAB 批量提取 STFT 特征的完整代码
单条信号处理完不够,要的是把所有样本拼成特征矩阵。下面这段代码把一段原始信号变成一行归一化特征向量:
fs = 12000; % 采样频率,与采集设备配置保持一致 winLen = 256; % 窗长,约 21.3 ms overlap = 192; % 重叠点数,占窗长 75% nfft = 512; % FFT 点数,决定频率栅格密度 segLen = 4096; % 单样本长度,按工况切分 function feat = stft_feature(x, fs, winLen, overlap, nfft) x = x - mean(x); % 去直流,避免 0 Hz 分量主导能量 w = hann(winLen); % 汉宁窗,抑制频谱泄漏 [S, ~, ~] = stft(x, fs, 'Window', w, ... 'OverlapLength', overlap, 'FFTLength', nfft); P = abs(S).^2; % 转成功率谱 P = P(1:nfft/2+1, :); % 只取单边谱 band = mean(P, 2).'; % 沿时间轴取均值,得到频带能量分布 feat = band / (sum(band) + eps); % 总量归一化,抵消载荷带来的幅值漂移 end featMat = zeros(0, nfft/2+1); label = strings(0, 1); for k = 1:numel(files) x = readmatrix(files{k}); for s = 1:floor(numel(x)/segLen) seg = x((s-1)*segLen+1 : s*segLen); featMat(end+1, :) = stft_feature(seg, fs, winLen, overlap, nfft); %#ok<SAGROW> label(end+1, 1) = classNames{k}; %#ok<SAGROW> end end逻辑说明:先按segLen把长信号切成等长样本,每条样本独立做 STFT,再沿时间轴对功率谱取均值——这一步把「时间」维度消掉,得到的是每个频点的平均能量。归一化用总量除法而不是 z-score,是因为幅值本身受载荷影响,而能量在各频带的相对分布才是故障特征。
参数说明:overlap必须小于winLen;nfft建议不小于winLen,否则频率栅格比窗的主瓣还粗,等于白补零;segLen至少要覆盖若干个冲击周期,太短会丢失周期性。
2.3 特征归一化与标签构造:避开量纲和类不平衡两个坑
257 维特征里,低频段的能量往往比高频段高两三个数量级。SVM 用的是欧氏距离,如果不做处理,它会几乎只盯着低频那几个点,高频共振带的信息被完全淹没。解决办法有两个:一是像上面那样做总量归一化,二是交给templateSVM的Standardize参数做 z-score。两者选一个即可,重复使用会削弱特征本身的物理含义。
标签构造要注意类不平衡。实际工况下正常运行样本占绝大多数,如果直接按采集比例切分,SVM 会把「全部判为正常」当成一个损失很小的解,准确率看着有 90%,实际上一个故障都抓不出来。常见做法是对每类限制相同的样本数,或者在fitcecoc里通过Prior参数指定先验。切分样本时按时间顺序连续切,不要随机打乱后切,否则同一段信号的相邻帧会同时落进训练集和测试集,交叉验证的损失会虚低。
提示:特征矩阵保存成
.mat时把featMat、label、fs、winLen一起存,后面 GUI 加载模型时要用到这些参数做一致性校验。
3. SVM 分类器训练与参数寻优:核函数、C 与 sigma 怎么定
特征准备好之后,剩下的事就是选一个分类边界。故障类别通常是四类以上,线性不可分的情况居多,直接上 RBF 核的多分类结构是最省事的起点。
3.1 为什么多分类故障诊断用 fitcecoc 而不是 fitcsvm
fitcsvm只解决二分类。四类故障要做一对一或一对多拆分,自己写循环再投票,代码量大且容易在概率输出上出错。fitcecoc是 MATLAB 自带的多分类封装,内部默认用一对一编码,把 C(4,2)=6 个二分类器拼起来,输出时用投票加置信度排序。它的接口和fitcsvm一致,可以通过templateSVM传入同一套核参数,后期想换核函数只改一处。
load('featureData.mat'); % 含 featMat (N x 257) 与 label (N x 1) cvp = cvpartition(label, 'KFold', 5, 'Stratify', true); % 分层,保证每折类别比例一致 t = templateSVM( ... 'KernelFunction', 'rbf', ... % 高斯核,处理非线性边界 'BoxConstraint', 10, ... % 惩罚系数 C 'KernelScale', 'auto', ... % 核宽 sigma,'auto' 用启发式估计 'Standardize', true); % 内部做 z-score,抵消量纲差异 mdl = fitcecoc(featMat, label, 'Learners', t, ... 'Coding', 'onevsone', 'CVPartition', cvp); loss = kfoldLoss(mdl); fprintf('5 折交叉验证损失: %.4f\n', loss);逻辑说明:cvpartition加Stratify是关键,它保证每一折里各类样本比例和总体一致,否则某一折可能完全没有滚动体故障样本,损失值会剧烈波动。CVPartition传进fitcecoc后,返回的是分区模型对象,kfoldLoss直接给出平均损失,不用自己写循环。
参数说明:Coding选onevsone时训练 6 个二分类器,样本量小时比onevsall更稳;样本量上万时onevsall训练更快。
3.2 RBF 核两个关键参数的物理含义与初始取值
RBF 核只有两个参数,但它们的组合决定了边界的形状。
BoxConstraint(C)控制对错分样本的惩罚力度。C 小,边界平滑,容错高,容易欠拟合;C 大,边界贴着训练点走,容易过拟合。KernelScale(sigma)控制单个样本的影响半径。sigma 小,每个点只影响自己周围一小片,边界碎裂;sigma 大,影响范围扩散,边界退化成近似线性。
| 参数 | 常用范围 | 偏小的表现 | 偏大的表现 |
|---|---|---|---|
| BoxConstraint | 1e-2 ~ 1e3 | 欠拟合,各类混在一起 | 过拟合,训练集满分测试集掉分 |
| KernelScale | 1e-2 ~ 1e2 | 边界破碎,噪声被当特征 | 边界过平滑,故障类无法分离 |
实战里有个省事的起点:先用'KernelScale','auto'让 MATLAB 按特征维度的启发式值估计,只调 C,从 1、10、100 各试一遍,观察交叉验证损失的变化趋势,再决定要不要联合寻优。很多时候 C 取 10 左右配上 auto 的 sigma 就已经能到 95% 以上,继续调收益有限。
3.3 用 K 折交叉验证与贝叶斯优化找参数组合
手调两个参数效率低,用bayesopt做自动搜索更快。它的目标函数就是交叉验证损失,搜索空间取对数尺度:
vars = [ ... optimizableVariable('BoxConstraint', [1e-2, 1e3], 'Transform', 'log'), ... optimizableVariable('KernelScale', [1e-2, 1e2], 'Transform', 'log')]; objfun = @(p) kfoldLoss(fitcecoc(featMat, label, ... 'Learners', templateSVM('KernelFunction', 'rbf', ... 'BoxConstraint', p.BoxConstraint, ... 'KernelScale', p.KernelScale, ... 'Standardize', true), ... 'Coding', 'onevsone', ... 'KFold', 5)); % 每轮重新分层,避免过拟合到固定划分 results = bayesopt(objfun, vars, ... 'MaxObjectiveEvaluations', 30, ... 'IsObjectiveDeterministic', false, ... 'Verbose', 1); bestC = results.XAtMinObjective.BoxConstraint; bestSigma = results.XAtMinObjective.KernelScale;逻辑说明:目标函数里用'KFold',5而不是固定的CVPartition,让每次评估重新分层划分。如果复用同一个划分,30 轮搜索会把参数往这个特定划分上过拟合,换一批数据后性能会掉。代价是每轮耗时更长,30 轮在几千样本规模下通常十几分钟能跑完。
参数说明:Transform','log'让搜索在对数空间均匀取点,否则大数值区间的采样点会过于稀疏;IsObjectiveDeterministic设为 false,因为交叉验证划分有随机性,同一组参数两次评估结果可能不同,贝叶斯优化需要按随机目标处理。如果手上装的是 matlab优化工具箱,也可以用fmincon在外层套交叉验证循环,但需要自己处理对数变换和边界约束,代码量更大。
4. MATLAB GUI 设计:把 STFT-SVM 诊断流程做成可交互工具
训练脚本跑通之后,现场同事不会愿意改代码里的文件名再点运行。把「加载信号 — 看时频图 — 出诊断结果」这条链路做成界面,才算把故障诊断代码变成能用的工具。
4.1 App Designer 与 GUIDE 的取舍
GUIDE 在较新版本里已经不推荐用于新项目,生成的.fig+.m双文件结构在版本管理里也不好处理。App Designer 用单个.mlapp文件,布局用网格自动排布,控件属性和回调函数在同一个代码视图里,导出和打包更顺。手上有历史 GUIDE 项目要维护可以继续用,新建一律建议 App Designer。
界面布局建议分三块:左上放「加载信号」和采样率、窗长输入框,右侧放两个坐标轴分别显示原始波形和 STFT 时频图,下方放核参数输入、训练按钮和结果文本区。matlab下载安装教程里那些工具箱勾选项,这里要确认信号处理工具箱和统计与机器学习工具箱都勾上了,缺一个都会在调用stft或fitcecoc时报未定义。
4.2 控件布局与回调函数写法
核心是三个回调:加载、特征提取与绘图、训练与预测。
% 加载按钮回调:读取信号并显示原始波形 function LoadButtonPushed(app, ~) [file, path] = uigetfile({'*.mat;*.csv;*.txt', '信号文件'}); if isequal(file, 0), return; end app.RawSignal = readmatrix(fullfile(path, file)); plot(app.UIAxes, app.RawSignal); title(app.UIAxes, '原始振动信号'); xlabel(app.UIAxes, '采样点'); ylabel(app.UIAxes, '幅值'); end % 特征提取回调:画时频图并存下特征向量 function ExtractButtonPushed(app, ~) fs = app.SampleRateEditField.Value; wl = app.WinLenEditField.Value; x = app.RawSignal - mean(app.RawSignal); [S, f, t] = stft(x, fs, 'Window', hann(wl), ... 'OverlapLength', round(wl*0.75), 'FFTLength', 512); P = abs(S(1:257, :)).^2; band = mean(P, 2).'; app.FeatVec = band / (sum(band) + eps); imagesc(app.UIAxes2, t, f, 10*log10(P + eps)); % 转 dB 显示,动态范围更清楚 axis(app.UIAxes2, 'xy'); xlabel(app.UIAxes2, '时间 / s'); ylabel(app.UIAxes2, '频率 / Hz'); colormap(app.UIAxes2, 'jet'); colorbar(app.UIAxes2); end % 训练/预测回调:加载模型并对当前特征分类 function TrainButtonPushed(app, ~) if isempty(app.FeatVec) app.ResultTextArea.Value = '请先提取特征'; return; end S = load('svmModel.mat'); % 里面存着 mdl 和特征归一化参数 [pred, score] = predict(S.mdl, app.FeatVec); [maxScore, ~] = max(score); app.ResultTextArea.Value = sprintf('诊断结果:%s\n置信度:%.3f', ... string(pred), maxScore); end逻辑说明:加载和特征提取分开成两个按钮,是为了让中间结果可检查——时频图长什么样,直接决定后面分类能不能做。10*log10是把功率谱转成分贝,线性刻度下高频微弱成分几乎看不见,转 dB 后共振带的位置一目了然,这也是 matlab图像处理里时频图渲染的常规做法。predict返回的score是各类的置信度,最大值低于 0.6 时建议在界面上给出「结果存疑」的提示,而不是硬报一个类别。
参数说明:控件命名用UIAxes、SampleRateEditField这类可读名,自动生成的app.UIAxes2之类在后面加控件时容易混淆;colorbar只加一次,重复调用会在图上叠出多条色条。
4.3 时频图渲染与诊断结果面板的联动
坐标轴和结果区之间的联动是体验的关键。可以在ExtractButtonPushed末尾顺便清空结果区,避免用户换了信号却还看着上一次的结论。时频图坐标轴要固定色标范围,否则每条信号自动缩放,正常和故障的图看起来颜色差不多,反而误导判断。
| 控件 | 回调 | 关联变量 | 注意点 |
|---|---|---|---|
| 加载按钮 | LoadButtonPushed | RawSignal | 读取后校验长度是否大于一个样本段 |
| 采样率输入 | 无 | SampleRate | 必须与采集配置一致,否则频率轴全错 |
| 提取按钮 | ExtractButtonPushed | FeatVec | 同时刷新时频图与清空结果区 |
| 训练按钮 | TrainButtonPushed | mdl | 加载.mat模型,校验特征维度一致 |
| 结果文本区 | 无 | — | 置信度低时给提示而非直接下结论 |
模型保存时把训练用的winLen、nfft、fs一起写进.mat。GUI 提取特征前先比对当前参数与模型参数,不一致就提示重新提取,能挡掉大部分「维度对不上」的报错。打包分发用 Application Compiler 生成安装包,依赖的运行时会在安装时一并处理,接收方不需要单独装 MATLAB。
5. 诊断准确率上不去时的排查顺序与可解释性技巧
交叉验证损失降不下来,先别急着换模型,九成问题出在特征或数据划分上。按下面的顺序查,比盲目调参快得多。
| 现象 | 优先排查 | 处理方式 |
|---|---|---|
| 所有类都判成样本最多的那一类 | 类不平衡 | 每类限同一样本数,或设Prior |
| 训练集接近满分、测试集很差 | 划分泄漏 | 按时间段切分,不要随机打乱后切 |
| 各类损失都高且接近 | 特征无区分度 | 检查窗长是否过长,冲击被平均掉 |
| 某两类反复混淆 | 特征频带重叠 | 按故障特征频率重划频带 |
| 换一批数据性能骤降 | 幅值未归一化 | 加总量归一化或Standardize |
频带重划是最容易被忽略、也最有效的一招。轴承外圈、内圈、滚动体的故障特征频率可以由几何参数算出,例如外圈通过频率 BPFO ≈ (n/2)·fr·(1 − (d/D)·cosθ),其中 n 为滚动体数,fr 为转频,d 为滚动体直径,D 为节圆直径,θ 为接触角。把 257 维的全频段能量改成只保留基频及其 2~3 次谐波附近的若干窄带,维度能降到 30 维左右,小样本下泛化明显更好,而且每个维度都能对应到一个物理量,可解释性远强于端到端的黑箱模型——这也是当前从信号到语义的时序解析大模型与可解释故障诊断方向反复强调的一点:特征层可解释,比事后加注意力图更可靠。
如果还想再压维度,可以在这批窄带特征上做 PCA,保留累计贡献率 95% 对应的主成分,通常落在 15~25 维。注意 PCA 的变换矩阵必须只用训练集拟合,再应用到测试集,否则同样构成泄漏。验证阶段除了看整体准确率,一定要用confusionchart把混淆矩阵画出来,重点看误判是发生在「内圈↔外圈」这种物理上本就相近的类别之间,还是发生在「正常↔故障」这种不该出错的地方——前者说明特征频带选得不够细,后者说明数据或标签本身有问题。
最后补一个工程上的小技巧:把 GUI 里的时频图坐标轴色标范围固定为训练集所有样本 dB 值的 5% 到 95% 分位数,正常与故障的图放在同一色标下对比,视觉差异和分类器的判断依据就能对上,现场排查误判时不用再回头翻代码。
本文还有配套的精品资源,点击获取