简介：这份MATLAB轴承故障诊断代码包面向机械振动信号处理与设备状态监测的学习者，围绕轴承故障特征提取与模式判别，整合Hilbert包络谱、Haar小波、数学形态学、时域无量纲参数及FFT等分析思路。压缩包共30个文件，约1.04MB，包含5个m脚本、11个dat数据文件及14张jpg结果图，脚本用于复现分析流程，数据与图像便于对照验证和结果展示。已有6610人学习下载，说明其在故障诊断入门与课程实践中具有一定参考价值。读者可据此理解不同方法在早期微弱故障识别、噪声抑制、周期性冲击提取及故障特征频率定位中的差异，并借助现成脚本完成振动数据处理、参数计算与频谱、包络谱观察，形成从理论方法到MATLAB实现的完整链路，适合作为课程设计、毕业设计或工程预研的基础素材。

1. 从一段"看不出问题"的振动信号说起

一台减速机换了新轴承，三个月后测振仍然报异常，可 FFT 幅值谱上干干净净，一根像样的峰都没有。问题不在轴承，而在于只看幅值谱根本抓不到早期冲击。这套轴承故障 matlab 代码要解决的就是这件事：先用 FFT 与希尔伯特包络谱把故障特征频率找出来，再叠加一层 haar 小波分解，把埋在噪声里的周期性冲击放大到能判读的程度。它面向手里已经有振动数据、想快速验证判据的工程师和做课题的学生，也适合当作诊断算法的基线工程——目录干净、链路完整、每一步都能单独拎出来替换。下面按我拆包时的顺序，从目录结构、信号链路一直讲到参数怎么改、判据怎么加。

2. 工程结构拆解：数据、算法、主脚本三层怎么分工

拿到这个包，最先要做的不是急着run，而是看清它把什么东西放在了哪里。这类诊断工程的翻车点，八成不是算法写错，而是数据、参数、调用顺序三者没对齐。我拆的时候习惯先画一张"谁调用谁"的图，再逐个文件读。

2.1 目录结构与文件职责

包里的组织方式很典型：数据单独放，算法函数按处理环节拆开，主脚本只负责串流程。实际文件命名可能略有出入，但职责分工基本是下面这几块。

这种"函数 + 主脚本"的结构，好处是每个环节都能单独拿一段信号喂进去试。坏处是 MATLAB 的工作目录机制很容易坑人：脚本在src/里跑，data/的相对路径就找不到。常见做法是在main.m开头把路径补上，而不是靠手动切目录。

% main.m 开头：把 src 和 data 都挂到搜索路径上，避免"函数未定义" rootDir = fileparts(mfilename('fullpath')); addpath(fullfile(rootDir, 'src')); addpath(fullfile(rootDir, 'data'));

mfilename('fullpath')拿到的是当前脚本的绝对路径，再往上取一层就是工程根目录。这样无论从哪个目录启动，路径都不会错。参数说明上，addpath加的是目录不是文件，所以新增函数只要丢进src/就自动可用。

2.2 数据格式与采样参数约定

诊断结果的可靠性，一半压在采样率和转速这两个数上。包里对数据格式的约定比较宽松，但有两件事必须自己确认：第一列到底是什么，采样率到底是多少。我见过太多人把带表头的 CSV 直接readmatrix，结果第一行全变成NaN，后面 FFT 出来一片噪声。

function [x, fs] = loadSignal(filePath, fs) % loadSignal 读取单通道振动信号并做最基础的前处理 % filePath : 数据文件路径，支持 csv / txt / mat % fs : 采样率(Hz)，由 config 传入，绝不从文件里猜 raw = readmatrix(filePath); % 纯数值读取，表头会被解析成 NaN x = raw(:, end); % 取最后一列：单列时就是它，两列时是幅值 x = x(:); % 强制列向量，防止后续维度对不上 bad = isnan(x); if any(bad) warning('检测到 %d 个无效点，已用线性插值补齐', sum(bad)); x(bad) = interp1(find(~bad), x(~bad), find(bad), 'linear', 'extrap'); end x = detrend(x, 'constant'); % 去直流分量，否则 0 Hz 会压满量程 x = x - mean(x); % 再减一次均值，双保险 x = x / max(abs(x) + eps); % 归一化到 ±1，便于跨数据比较阈值 end

这段代码里有三个决策值得说。取raw(:, end)而不是raw(:, 1)，是因为很多采集软件导出时第一列是时间戳，取最后一列更稳；如果数据只有一列，end就等于 1，两种情况都覆盖。detrend(x,'constant')和减均值看着重复，实际差别在于前者对缓慢漂移更有效，两个都做不会出错。最后归一化这一步，是为了让峭度、峰值这些无量纲指标在不同量程的数据之间可比——如果你要做绝对幅值判断（比如 ISO 振动烈度），这行得去掉。

fs明确作为入参传进来，不从文件里读。这是血泪经验：靠文件头或靠"点数除以时长"反推采样率，只要数据段被截过，算出来的fs就是错的，后面所有频率轴全部偏移。

2.3 主流程脚本的调用顺序

main.m的职责就一件事：按正确顺序把函数串起来，并把中间结果画出来。顺序错了，后面全白搭。

%% main.m 主流程 cfg = config(); % 1. 全部参数集中在这里 [x, fs] = loadSignal(cfg.file, cfg.fs); % 2. 读数据 + 去直流 ff = faultFreqs(cfg.bearing, cfg.rpm); % 3. 先算理论特征频率 [f, S] = envelopeSpectrum(x, fs, cfg.band); % 4. 带通 + 包络谱 [c, l] = haarDecompose(x, cfg.haarLevel); % 5. Haar 多层分解 d3 = wrcoef('d', c, l, 'haar', 3); % 取第 3 层细节重构 kurt = kurtosis(d3); % 6. 峭度指标 [peakFreq, snr] = findFaultPeak(S, f, ff.BPFO); % 7. 特征频率邻域找峰 fprintf('BPFO=%.2f Hz, 实测谱峰=%.2f Hz, 峭度=%.2f\n', ... ff.BPFO, peakFreq, kurt);

调用顺序上，第 3 步算特征频率必须排在第 7 步找峰之前，而第 3 步又依赖cfg.rpm——转速是唯一无法从信号里直接得到的量，必须人工给。这里有个容易忽略的点：faultFreqs接收的是转速 rpm，内部要除以 60 转成转频 Hz，单位换算放在函数内部比放在config里更不容易乱。

wrcoef('d', c, l, 'haar', 3)里的3是层号，1 是最细的那层。想拿哪一层的冲击成分，就把这个数字换掉，后面第 4 章会讲怎么定这个层数。整条链路里，只有cfg是需要频繁改的，其他函数基本一次写定。

提示：第一次跑不要急着看结论，先把x的时域波形和S的包络谱画出来，确认波形不是一条直线、谱图不是一片糊，再谈诊断。

3. FFT 包络谱：从原始振动到故障特征频率

轴承故障的物理本质是"冲击"：滚动体碾过剥落坑，每次都会激起一个宽带脉冲，脉冲重复的频率就是故障特征频率。问题是这个重复频率往往在几十到几百赫兹，而冲击能量集中在几千赫兹的共振带上——直接对原始信号做 FFT，看到的是共振带被调制后的宽包，特征频率本身被"藏"在调制里。包络解调做的就是把这层调制拆出来。

3.1 去直流、去趋势与带通滤波

带通滤波这一步决定了成败。滤波带不是随便选的，它的目标是把共振频带留下来、把低频的转频谐波和高频噪声切掉。选窄了，冲击的边带被削掉，包络失真；选宽了，其他部件的干扰混进来，谱峰被淹没。

function [b, a] = designBand(fs, band, order) % designBand 设计零相位带通滤波器 % fs : 采样率 % band : [fLow fHigh]，共振频带，单位 Hz % order : 滤波器阶数，4 阶够用，太高会引起数值不稳定 if band(2) >= fs/2 error('带通上限 %.1f Hz 超过奈奎斯特频率 %.1f Hz，检查采样率', ... band(2), fs/2); end wn = band / (fs/2); % 归一化到 0~1 [b, a] = butter(order, wn, 'bandpass'); end

butter返回的是 IIR 系数，4 阶带通实际是 8 阶原型，滚降足够陡而且不会像高阶那样在低频段数值发散。wn归一化必须除以fs/2，这一步写错是新手最常见的错误之一——写成band/fs，滤波带就直接偏移一倍。

真正执行滤波时用filtfilt而不是filter：

[b, a] = designBand(fs, cfg.band, 4); xf = filtfilt(b, a, x); % 零相位滤波，正反各滤一次

filtfilt会把相位延迟抵消掉，代价是边界附近有暂态。做包络分析时相位一致性很重要，因为后面要用冲击的时间位置，filter引入的群延迟会让冲击"挪位"。如果你的数据段很短（不到滤波器阶数的 10 倍），边界暂态会吃掉可观的一段，那就掐头去尾再分析。

带通中心怎么定？常见做法是先看原始信号的功率谱，找 2 kHz～8 kHz 区间里那个凸起的宽包；如果谱图很平，就试三组带（比如 1–3k、2–5k、4–8k），哪组的包络谱峰最干净用哪组。更严谨的办法是谱峭度图，但那属于额外工程量，先跑通再说。

3.2 希尔伯特包络与包络谱

包络的数学定义是解析信号的模。MATLAB 里hilbert返回的已经是解析信号（复数列），取abs就得到包络，不需要自己造虚部。

function [f, S, env] = envelopeSpectrum(x, fs, band) % envelopeSpectrum 带通 -> 希尔伯特包络 -> 单边幅值谱 [b, a] = designBand(fs, band, 4); xf = filtfilt(b, a, x); % 零相位带通 env = abs(hilbert(xf)); % 解析信号取模 = 包络 env = env - mean(env); % 去掉包络直流，关键一步 env = env .* tukeywin(numel(env), 0.05); % 轻微加窗，压边带泄漏 N = 2^nextpow2(numel(env)); % 补零到 2 的幂次，加速 FFT Y = fft(env, N) / numel(env); % 幅度归一化 f = (0:N/2-1)' * fs / N; % 单边频率轴 S = 2 * abs(Y(1:N/2)); % 单边谱，幅值乘 2 end

参数上几个细节：2^nextpow2补零不改变频率分辨率，真正的分辨率还是fs/numel(env)，补零只是让 FFT 算得快、曲线插值更密。所以别误以为补零能"提高精度"，想提高分辨率只能加大数据长度。

env - mean(env)这行是整个函数里最容易被忽略、后果最严重的一行。包络恒为正，均值很大，如果不去掉，FFT 之后 0 Hz 处会出现一根比故障峰高两三个数量级的巨峰，纵轴一压缩，什么细节都看不见了。我第一次跑的时候就栽在这，盯着图看了半天以为算法有 bug。

tukeywin加 5% 的锥形窗，是为了让包络两端平滑过渡到零，减少频谱泄漏。窗太狠（比如 0.5）会把有效数据压掉太多，5% 是经验值。

拿到S之后，判读方法很直接：在理论故障特征频率f0附近 ±2% 的窗口里找最大值，再算这个峰值相对于窗口附近背景中位数的比值，作为信噪比。比值大于 3 就认为有峰，小于 2 就当没检出。

function [peakFreq, snr] = findFaultPeak(S, f, f0, tol) % 在 f0 邻域内找谱峰，并给出相对背景的信噪比 if nargin < 4, tol = 0.02; end win = f >= f0*(1-tol) & f <= f0*(1+tol); % 目标窗口 guard = f >= f0*0.5 & f <= f0*1.5 & ~win; % 参考背景，避开目标 [pk, k] = max(S .* win); % 窗口外置零后取最大 peakFreq = f(k); snr = pk / (median(S(guard)) + eps); % eps 防零除 end

S .* win这个技巧很实用：逻辑索引当掩码乘上去，窗口外直接变 0，max就只在窗口内取值。snr用中位数而不是均值做分母，是因为背景里可能混着其他谐波峰，均值会被拉高。

3.3 故障特征频率计算与谱峰比对

理论特征频率由轴承几何决定，四个公式必须记牢，写反一个就全盘皆输：

function ff = faultFreqs(brg, rpm) % brg : 结构体，字段 d(滚动体直径) D(节圆直径) n(滚动体数) alpha(接触角,度) % rpm : 轴转速，转/分 fr = rpm / 60; % 转频 Hz，别把 rpm 直接当转频用 r = brg.d / brg.D * cosd(brg.alpha);% 无量纲比值，Haar 分析里也用得上 ff.BPFO = brg.n/2 * fr * (1 - r); % 外圈故障频率 ff.BPFI = brg.n/2 * fr * (1 + r); % 内圈故障频率 ff.BSF = brg.D/(2*brg.d) * fr * (1 - r^2); % 滚动体自转故障频率 ff.FTF = fr/2 * (1 - r); % 保持架故障频率 end

cosd用的是角度制，如果参数表给的是弧度就要换cos，这个坑很隐蔽——接触角通常写 0°、15°、25°，用cos会得到 0.999 而不是 1，误差小到看不出来，但深沟球轴承的alpha=0时结果应该完全等于n/2*fr，一旦对不上就说明单位用错了。

判读的时候，外圈故障看 BPFO 及其 2 倍频，内圈看 BPFI，同时要留意转频fr的边带——内圈故障因为随轴旋转，调制更明显，BPFI 两侧常有BPFI ± fr的边频。如果只看到fr的谐波而没有任何特征频率成分，那多半是不平衡或不对中，不是轴承问题。

一个实用习惯：把ff里四个频率竖线画到包络谱上，一眼就能看出哪根峰落在哪条线上。这一步比看数字快得多，也更容易发现"峰位整体偏移"这种系统性错误。

4. Haar 小波分解：把早期微弱冲击从噪声里挖出来

包络谱对中晚期故障很好用，但早期剥落面积小、冲击弱，包络谱上可能只是一个勉强过阈值的小包。这时候小波变换的价值就出来了：它在时频两个维度上同时定位，能把"短时、周期性、宽频带"的冲击成分单独拎出来。

4.1 为什么在这套链路里选 Haar

小波基的选择是个老生常谈的问题。db4、sym8、morl 都有人用，这套代码里选 Haar（等价于 db1），理由有三个，都跟"冲击检测"这个具体任务有关。

第一，Haar 的支撑长度最短，只有两个点，时域定位精度最高。冲击信号的特征就是"突变"，用长支撑的基函数去匹配，等于用一个平滑窗去糊一个尖峰，冲击的起跳点会被抹开。第二，Haar 是正交基，分解后的各层细节系数能量可以直接相加，不用考虑基函数之间的相关性，做能量占比统计时省事。第三，Haar 的滤波器系数只有[1/√2, 1/√2]和[1/√2, -1/√2]，计算量极小，批量跑几百个文件时差距明显。

代价是频域局部性差，频率划分是"二分"的、不平滑，相邻层之间有重叠。所以它适合做冲击检测，不适合做精细的频率估计——频率估计交给 FFT 那一路。

4.2 多层分解与单支重构

分解本身一行就够，关键是层数怎么定。层数决定了每一层对应的频带宽度：第k层细节系数覆盖的是fs/2^(k+1)到fs/2^k这个区间。要让目标频带落在某一层里，层数就得按采样率反推。

function [c, l] = haarDecompose(x, level) % haarDecompose 多层一维小波分解 % level : 分解层数，需满足 numel(x) >= 2^level * 滤波器长度 wname = 'haar'; maxLev = wtdemm('max', numel(x), wname); % 允许的最大层数 if level > maxLev warning('层数 %d 超过上限 %d，已自动截断', level, maxLev); level = maxLev; end [c, l] = wavedec(x, level, wname); % c 是全部系数，l 是各段长度 end

wavedec返回的c是把近似系数和各层细节系数首尾拼起来的一维长向量，l记录每段的长度。直接操作c很容易下标错位，所以取某一层时用detcoef或wrcoef，别自己切片。

function d = haarDetailBand(c, l, level) % haarDetailBand 重构第 level 层细节系数，返回与原信号等长的分量 d = wrcoef('d', c, l, 'haar', level); end

wrcoef('d', ...)做的是单支重构：把第level层细节系数保留，其余全部置零，再逆变换回原尺度。输出的长度和原信号一致，可以直接和原信号叠画对比。要特别注意，单支重构后的幅值不等于原始信号里该频带的真实幅值——它只是一个"该层系数对应的时域波形"，做阈值判断只能用相对值（比如相对整段信号的峭度、相对包络谱峰值），不能拿绝对值去和 ISO 标准比。

层数怎么选？我一般这么算：目标共振带取 2–5 kHz，采样率 20 kHz，那么第 2 层细节覆盖 2.5–5 kHz，第 3 层覆盖 1.25–2.5 kHz。两层都试，看哪一层的重构信号里周期性冲击更清楚（时域上能看到等间隔的尖峰簇）。这比死记硬背公式靠谱。

4.3 峭度与频带能量作为辅助判据

小波分解出来的东西，最终要变成能比较的数字。两个指标最常用：峭度和频带能量占比。

峭度衡量分布的"尖峭程度"，正态信号峭度约为 3，冲击成分越多峭度越高。它对早期冲击特别敏感，缺点是也容易被偶发的单点干扰骗到。

function [kurt, ratio] = bandIndicators(x, d) % kurt : 重构分量的峭度 % ratio : 该分量能量占原信号能量的比例 kurt = kurtosis(d); % 正态信号基准值 3 ratio = sum(d.^2) / (sum(x.^2) + eps); end

kurtosis默认用的是有偏估计（不减去 3），所以判断基准是 3 而不是 0。如果你用的是kurtosis(d) - 3那种超额峭度，基准就变成 0，两种写法混用会让阈值全部失效——这类"基准不一致"的错误在多人协作时特别常见。

ratio用来说明"这层到底有没有料"。如果某层能量占比不到 1%，即使峭度很高也要怀疑是噪声尖峰；能量占比合理（5%–30%）且峭度明显高于 3，才算是可信的冲击成分。

实际判据常见做法是"双门限"：包络谱上特征频率处的信噪比大于 3，同时对应小波层的峭度大于 4，两个条件都满足才报"疑似故障"。单靠一个指标误报率会明显上升，尤其是在变工况数据上。

5. 跑起来与环境排查：五个让结果失真的坑

代码能跑出图，不等于结果可信。诊断类工程最大的特点是"错了也出图"，一条平滑的谱线看着挺专业，实际上可能全是伪影。这一章把我踩过的、以及帮别人看代码时反复见到的几类问题摊开讲。

5.1 最小跑通路径与环境

先确认环境。这套代码只用到 Signal Processing Toolbox 和 Wavelet Toolbox，butter、filtfilt、hilbert属于前者，wavedec、wrcoef、wtdemm属于后者。缺任何一个，报错信息会是"未定义函数"而不是"缺少工具箱"，容易误判成路径问题。

% 一键检查依赖是否就位 deps = {'butter','filtfilt','hilbert','wavedec','wrcoef','kurtosis'}; missing = deps(cellfun(@(f) exist(f, 'file') == 0, deps)); if isempty(missing) disp('依赖齐全'); else fprintf('缺失: %s\n', strjoin(missing, ', ')); end

跑通的最短路径是四步：把数据文件放进data/，在config.m里改cfg.file、cfg.fs、cfg.rpm和cfg.band，运行main.m，看图。cfg.rpm必须填真实转速，填个大概值（比如 1500 而实际 1480）会让特征频率偏 1.3%，刚好落在 ±2% 的搜索窗边缘，找峰结果时对时错——这类问题最难查，因为"有时候能出结果"。

数据长度也有要求。做包络谱至少要有 10 个以上的故障冲击周期，也就是fs/BPFO * 10个采样点；采样率 20 kHz、BPFO 100 Hz 的话，至少要 2000 点，实际建议 1 万点以上。数据太短，频率分辨率不够，谱峰展宽到几个 bin，max取到的位置随机性很大。

5.2 五条踩坑记录

第一坑：包络谱上 0 Hz 一根巨峰，右侧什么都看不见

现象是画出来的谱图左边一根顶到画布顶部，其余部分几乎是平的。原因是包络信号恒为正值，均值远大于交流分量，FFT 之后直流项占了绝大部分能量，纵轴自动缩放把细节全压没了。解决很简单，求 FFT 之前减掉包络均值，再顺手高通一下（截止 5 Hz 足够）。判断是否已经处理过的方法：看env的均值是否接近 0。

第二坑：特征频率对不上，整体偏了几十赫兹

现象是谱峰位置和理论 BPFO 差得离谱，但换个数据又能对上。原因通常是转速单位搞错，把 rpm 直接当转频用了。rpm = 1500对应的转频是 25 Hz，如果直接拿 1500 去乘，算出来的 BPFO 会大 60 倍，完全对不上。还有一种情况是采样率填错，频率轴整体缩放，这时候不光是特征频率，所有峰的位置都不对。解决办法是先用一段已知频率的标定信号（或者转频本身）验证频率轴——转频峰应该出现在rpm/60处，对不上就说明fs有问题。

第三坑：换一段数据结果全乱，峭度从 5 跳到 20

现象是同一台设备换个时间段采的数据，指标差了好几倍。原因往往是固定了分解层数却没管采样率变化。第 3 层细节在 20 kHz 采样下覆盖 1.25–2.5 kHz，换成 50 kHz 采样就变成 3.1–6.25 kHz，完全跑到另一个频带去了。解决办法是把层数从"写死的常数"改成"由fs和目标频带算出来"，写个小函数pickLevel(fs, bandCtr)返回最匹配的层号。

第四坑：小波重构后的幅值比原信号小一大截

现象是把wrcoef的结果和原信号叠画，发现幅值只有原来的几分之一，怀疑代码算错了。原因这是单支重构的正常行为：只保留一层系数，其余置零，能量自然只剩那一小部分，幅值小是必然的。解决办法是别用绝对幅值做判断，改用峭度、能量占比这类相对指标，或者同时对多层重构后求和再比较。

第五坑：CSV 首行有表头，第一行数据变成 NaN

现象是readmatrix读进来第一行全是NaN，detrend之后整段信号变成NaN，图上一片空白。原因是导出软件默认带列名。解决办法是读完之后做一次isnan检查和插值补齐，或者在readmatrix里用'NumHeaderLines'参数跳过表头。前者更稳，因为它能同时处理中间偶发的缺采样点。

注意：这五类问题的共同点是"不报错、只出图"。所以每换一个数据集，先看时域波形形状对不对、再看转频峰在不在正确位置，最后才看特征频率，这个检查顺序能挡掉大部分误判。

6. 二次开发：把参数抽干净，把判据换成自己的

这套代码当基线用完，接下来必然要改。改之前先做一件事：把所有会变的量从函数里搬到config.m。我见过太多人直接在envelopeSpectrum.m里把band写成硬编码的[2000 5000]，换个采样率就全崩，而且过两周自己也忘了改过哪。

6.1 参数集中与配置抽离

config.m至少要有这几项：文件路径、采样率、转速、带通上下限、小波层数、轴承几何参数、判据阈值。用结构体返回，主流程里只出现cfg.xxx，不出现任何魔法数字。

function cfg = config() % 所有可变参数集中在此，换设备只改这一个文件 cfg.file = fullfile('data', 'vib_01.csv'); cfg.fs = 20000; % 采样率 Hz cfg.rpm = 1480; % 实测转速，别填铭牌值 cfg.band = [2000 5000]; % 共振频带，由功率谱观察确定 cfg.haarLevel = 3; % 目标频带落在第 3 层 cfg.bearing.d = 7.94e-3; % 滚动体直径 m cfg.bearing.D = 39.04e-3; % 节圆直径 m cfg.bearing.n = 9; % 滚动体个数 cfg.bearing.alpha = 0; % 接触角，度 cfg.th.snr = 3.0; % 包络谱信噪比门限 cfg.th.kurt = 4.0; % 峭度门限 end

单位统一用米和赫兹，别混着毫米用。d和D在公式里是比值，单位一致就抵消了，但一旦有人把d写成毫米、D写成米，比值就差 1000 倍，特征频率直接飞到天上。我一般会在faultFreqs里加一句合理性检查：BPFO应该在3*fr到12*fr之间，超出就报错。

6.2 批量跑与判据验证

单文件跑通之后，下一步是批量。把main.m拆成runOne(cfg)返回一个结构体结果，外层套一个循环扫目录，最后汇总成表。这样几百个文件跑下来，异常的那几条一眼就能挑出来。

function res = runOne(cfg) % 单文件诊断，返回结构化结果，便于批量汇总 [x, fs] = loadSignal(cfg.file, cfg.fs); ff = faultFreqs(cfg.bearing, cfg.rpm); [f, S] = envelopeSpectrum(x, fs, cfg.band); c = haarDecompose(x, cfg.haarLevel); d = wrcoef('d', c, [], 'haar', cfg.haarLevel); [res.peakFreq, res.snr] = findFaultPeak(S, f, ff.BPFO); res.kurt = kurtosis(d); res.verdict = res.snr > cfg.th.snr && res.kurt > cfg.th.kurt; end

改判据的时候有个验证技巧：拿一段确认正常的数据，人为在时域上按 BPFO 的周期叠加几个窄脉冲（幅值取背景噪声标准差的 2–3 倍），再跑一遍流程，看snr和kurt能不能被推过门限。这个半仿真测试能告诉你门限定得是高还是低——如果 2 倍标准差注入都检不出，说明门限太严；如果纯噪声也能过，说明太松。比直接拿现场数据调参可靠得多。

wrcoef的第二个参数传[]是允许的，只要wavedec时用的是默认模式，很多版本会自动沿用，但更稳妥的写法是把l一起传出来。这一点在不同 MATLAB 版本上表现不一致，属于典型的兼容性坑。

最后是扩展方向。想加倒频谱，在envelopeSpectrum后面接一层real(ifft(log(abs(fft(env)))))就能看出谐波族的间隔；想加谱峭度选带，把cfg.band换成自动搜索结果，接口不用动；想上深度学习，runOne输出的结构体正好当特征表用。这套代码的价值不在算法多新，而在于每个环节都能单独替换、替换后还能验证。

从那以后我每次拿到一批新数据，都强制先跑一遍"转频峰定位"——只看fr = rpm/60那根线在不在正确位置。在不在，决定了后面所有特征频率的判断有没有意义。希望帮到你。

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