简介面向机械振动信号处理与设备状态监测学习者的轴承故障诊断 MATLAB 代码包围绕滚动轴承振动数据的故障特征提取与类型判别展开适合具备一定信号处理基础、希望把课堂理论落到代码实践的高年级本科生、研究生及工程技术人员参考。压缩包共 30 个文件约 1.04MB以 5 个 m 脚本承载 FFT、Hilbert 包络谱、Haar 小波与数学形态学等分析流程11 个 dat 文件提供振动数据14 张 jpg 用于结果图或谱图对照便于按“数据—算法—图形”线索快速梳理。内容覆盖时域无量纲参数、快速傅里叶变换、包络解调、小波多分辨率分解和形态学滤波去噪可帮助读者理解内圈、外圈及滚动体故障特征频率的提取思路。已有 6610 人学习适合用作课程设计、实验复现和故障诊断入门练手通过运行与修改脚本积累信号预处理、特征构造和结果判读经验。1. 一段 12k 采样的振动波形摆在面前这套 MATLAB 轴承诊断代码在算什么轴承故障诊断落到工程上其实就一句话从一段看似杂乱的振动波形里把那个反复出现的冲击周期找出来。拿 matlab 打开现场采回来的 12k 采样数据屏幕上是一片看不出规律的毛刺直接 FFT 只见一堆转频谐波真正管用的往往是包络解调加 haar 小波去噪这一套组合拳。这份代码包做的就是这件事——读.mat、去噪、算频谱、算包络谱、对照理论故障特征频率给出结论一条流水线跑完。它适合两类人手上已经有振动数据、想快速验证“到底是内圈还是外圈”的现场工程师以及刚接触故障诊断、需要一份能跑通的最小可用骨架的学生和转岗开发者。代码量不大但把数据装载、分段、频谱、包络、特征提取、批量落表这几段都占全了属于那种“跑一遍就懂、改两行就能用在别处”的结构。2. 工程结构拆解从 .mat 原始数据到诊断结论的完整链路这份代码最容易被低估的地方是目录组织。很多人拿到压缩包第一反应是打开main.m从头读到尾结果被几十行参数和绘图代码糊住反而看不清数据是怎么流的。把它按“读入—预处理—变换—判决”四段拆开看脉络就很清楚。2.1 目录长什么样每个文件负责哪一段解压后大致是这样一棵树bearing_fault_diag/ ├── main.m # 单文件入口一条信号跑完全流程 ├── batch_run.m # 批量跑整个 data/ 目录输出特征表 ├── data/ # 原始振动数据每个工况一个 .mat │ ├── normal_1730rpm_0hp.mat │ ├── inner_1730rpm_0hp.mat │ ├── outer_1730rpm_0hp.mat │ └── ball_1730rpm_0hp.mat ├── func/ │ ├── load_signal.m # 统一读取兼容不同变量名 │ ├── segment_signal.m # 按点数切段、加窗 │ ├── fft_spectrum.m # 单边幅值谱 │ ├── envelope_demod.m # 带通 希尔伯特包络谱 │ ├── haar_denoise.m # haar 小波阈值去噪 │ ├── fault_freq.m # 四个理论故障特征频率 │ ├── feat_extract.m # 时域/频域特征 │ └── plot_helper.m # 绘图封装统一坐标轴 ├── result/ # 特征表、图片输出 └── README.mdfunc/下每个文件只做一件事函数签名统一以“一段信号 采样率”开头返回“频率轴 幅值”或者“一个特征结构体”。这种切法的好处是二次开发时替换其中一环不会牵连其它模块比如你想把 haar 换成 db4只动haar_denoise.m一行字符串就行。常见做法是再补一个config.m把参数集中起来避免main.m和batch_run.m各写一份转速。目录/文件职责改动频率data/原始振动信号高换数据就换这里func/可复用的算法单元中调参或换算法main.m单条信号演示流程低batch_run.m批量统计与落表中加特征列result/输出产物不手动改2.2 cfg 参数配置那几个数决定结论对不对整个流程的准确性一半压在几个物理参数上。转速错了后面算出来的理论故障频率全偏包络谱上明明有峰你也对不上号。% main.m 顶部的参数区batch_run.m 里也有一份务必保持一致 cfg.fs 12000; % 采样频率 Hz按采集卡实际设置填 cfg.rpm 1730; % 实测转速用转速计别抄铭牌额定值 cfg.n_ball 9; % 滚动体个数 cfg.d_ball 7.94e-3; % 滚动体直径 m cfg.D_pitch 39.04e-3; % 轴承节径 m外圈滚道中心到内圈滚道中心 cfg.contact 0; % 接触角深沟球轴承取 0 cfg.seg_len 4096; % 每段点数决定频率分辨率 cfg.overlap 0.5; % 相邻段重叠比例 cfg.level 5; % haar 分解层数 cfg.band [2000 5000];% 共振解调带通范围 Hzfs从采集卡设置里来不是猜的。rpm也一样——铭牌 1730 rpm 是额定值实际带载后可能差几十转这一差在包络谱上就是好几赫兹的偏移用眼睛找峰时够你难受半天。seg_len取 4096是因为 12000/4096 ≈ 2.93 Hz 的频率分辨率刚好能把 100 Hz 附近的外圈特征频率和它的边频分开取 1024 点分辨率只有 11.7 Hz边频直接被抹平。提示d_ball和D_pitch查轴承手册单位统一成米。手册上给的是毫米忘了换算会让特征频率错三个数量级这是新手最常踩的一脚。2.3 load_signal 与信号分段维度、直流、长度三件套原始数据五花八门变量名可能叫vib、DE、x、data长度也可能不是 2 的整数次幂。load_signal.m干的活就是把这些差异吞掉。function s load_signal(fp) % 统一读取 .mat兼容常见的几种变量命名 raw load(fp); fn fieldnames(raw); keys {vib,DE,x,data,signal}; name ; for i 1:numel(keys) hit strcmpi(fn, keys{i}); % 不区分大小写匹配 if any(hit) name fn{hit}; break; end end if isempty(name), name fn{1}; end % 兜底取第一个变量避免报错中断 v double(raw.(name)); % 统一转 double防止 int16 参与 fft v v(:); % 压成列向量消除行列不一致 v v - mean(v); % 去直流否则 0 Hz 处一个大峰 s.vib v; s.src fp; end三行处理的每一行都有原因。double是因为很多采集设备存的是int16直接送进fft会参与整型运算幅值溢出或者被截断v(:)是防止有的文件存成行向量、有的存成列向量拼起来的时候维度对不上去均值则是为了压掉 0 Hz 的直流分量不然频谱第一根柱子能顶到天上纵轴一自动缩放其它频率全被压成一条直线。分段在segment_signal.m里逻辑简单但别省function segs segment_signal(x, seg_len, overlap) % 按固定点数切段末尾不足一段的直接丢弃 step round(seg_len * (1 - overlap)); n floor((numel(x) - seg_len) / step) 1; segs zeros(n, seg_len); w hann(seg_len).; % 汉宁窗减小频谱泄漏 for k 1:n idx (k-1)*step (1:seg_len); seg x(idx).; segs(k,:) (seg - mean(seg)) .* w; % 每段单独去均值再加窗 end end加窗会带来幅值衰减汉宁窗的相干增益是 0.5严格做幅值标定时要乘 2 补回来。做故障有无的判断时通常不做这个补偿因为看的是峰的位置和相对高低不是绝对幅值。这个取舍要先想明白否则换了加窗函数发现幅值腰斩会误以为代码坏了。2.4 main.m 主流程一条信号走完五个处理站入口脚本本身很薄真正的工作都在函数里。这种“薄 main”的写法值得学读的人一眼能看清数据流向。%% main.m —— 单条信号诊断演示 clear; clc; close all; cfg get_cfg(); % 参数集中管理见 2.2 fp fullfile(data,outer_1730rpm_0hp.mat); s load_signal(fp); x s.vib.; [xd, ~, ~] haar_denoise(x, cfg.level); % 1. 小波去噪 [f, P1] fft_spectrum(xd, cfg.fs); % 2. 单边幅值谱 [fe, eP] envelope_demod(xd, cfg.fs, cfg.band); % 3. 包络谱 ft fault_freq(cfg); % 4. 理论特征频率 feat feat_extract(xd, cfg.fs); % 5. 特征向量 fprintf(转速 %.0f rpm转频 %.2f Hz\n, cfg.rpm, ft.fr); fprintf(理论 BPFO %.2f Hz / BPFI %.2f Hz / BSF %.2f Hz\n, ... ft.bpfo, ft.bpfi, ft.bsf); fprintf(时域峭度 %.3f峰值因子 %.3f\n, feat.kurt, feat.crest); plot_helper(f, P1, fe, eP, ft, s.src); % 三联图频谱 包络谱 标注按 1730 rpm、9 个滚动体、直径比 7.94/39.04 算下来转频 28.83 Hz外圈通过频率约 103.4 Hz内圈约 156.1 Hz滚动体自转约 67.9 Hz保持架约 11.5 Hz。这几个数就是后面找峰的“标准答案”。包络谱上如果在 103 Hz 附近出现明显峰并且 206 Hz、309 Hz 处有它的谐波基本可以判定外圈有问题如果峰出现在 156 Hz 且两侧带 28.8 Hz 的边频那就是内圈——内圈随轴转故障点进出载荷区天然带转频调制。注意plot_helper里把包络谱纵轴做了归一化不同文件放一张图上才有可比性。如果你自己改代码去掉了归一化几张图叠在一起会完全看不出峰。3. FFT 与包络解调把埋在共振带里的故障频率挖出来新手最容易犯的错是拿原始信号直接 FFT然后盯着 100 Hz 附近找峰找不到就怀疑数据有问题。实际上轴承的局部损伤产生的是一串冲击冲击激励的是轴承座和壳体的结构共振能量集中在几千赫兹的高频段。故障特征频率本身并不直接出现在频谱上它只以“重复频率”的形式调制着那个高频载波。要把它提取出来必须走包络解调这条路。3.1 直接 FFT 为什么看不到故障频率一个局部剥落的滚道每被滚动体碾过一次就产生一次冲击。冲击的频谱很宽能量铺在几百到上万赫兹而冲击的重复频率——也就是我们关心的 BPFO 或 BPFI——只是这个宽频带的包络变化频率。原始频谱里高频共振区的幅值被这些冲击抬高但你要找的那个 100 Hz混在转频谐波、不平衡、不对中这些低频成分里被压得几乎看不见。更麻烦的是信噪比。现场数据里齿轮啮合、电机电磁噪声、结构共振都会盖过轴承信号。所以合格的流程一定是“先去噪、再带通选共振频带、然后做包络、最后对包络做 FFT”。少一步很多时候就是从“能看出来”变成“看不出来”。血泪经验是不带通直接做包络包络谱上往往只剩转频和它的倍频看着像模像样其实什么都没诊断出来。3.2 单边幅值谱的正确写法fft_spectrum.m是整个包里被调用次数最多的函数也是最容易写错的一个。写错的典型症状是幅值差 2 倍或者频率轴对不上。function [f, P1] fft_spectrum(x, fs) % 返回单边幅值谱频率轴单位 Hz x x(:).; % 强制行向量 x x - mean(x); % 去直流重复一次以防调用方忘了 N numel(x); Y fft(x); P2 abs(Y) / N; % 双边谱幅值先除以点数 nh floor(N/2) 1; % 正频率部分取到奈奎斯特 P1 P2(1:nh); if mod(N, 2) 0 P1(2:end-1) 2 * P1(2:end-1); % 偶数点去掉 0 和 Nyquist 后翻倍 else P1(2:end) 2 * P1(2:end); % 奇数点Nyquist 不单独占一格 end f fs * (0:nh-1) / N; % 频率分辨率 fs/N end两个细节值得掰开说。第一个是除以Nfft出来的结果是“和”不是“平均”不除N的话点数越多幅值越大换个段长结论就变了。第二个是乘 2实信号的负频率部分是正频率的镜像能量被劈成了两半取单边谱就得把正频率那半翻倍补回来。至于偶数点和奇数点的区别只在奈奎斯特那一格——偶数点时它是独立的最后一格不能翻倍翻了就重复计一次能量。f fs*(0:nh-1)/N这一行决定了你能分辨多细的频率间隔。12000 Hz 采样、4096 点分辨率 2.93 Hz。你做一次 ±1% 的特征频率搜索在 103 Hz 附近就是 ±1 Hz 都不到还没一个频率分辨率宽搜索窗口至少开到 ±3 个分辨率也就是 ±9 Hz否则真实峰偏离一点就漏了。3.3 希尔伯特包络解调带通、包络、再 FFTfunction [f, envP] envelope_demod(x, fs, band) % 带通选共振频带 - 希尔伯特包络 - 包络谱 x x(:). - mean(x); if nargin 3 ~isempty(band) Wn band / (fs/2); % 归一化到奈奎斯特 [b, a] butter(4, Wn, bandpass); % 4 阶足够再高容易数值不稳 x filtfilt(b, a, x); % 零相位滤波避免包络整体时移 end env abs(hilbert(x)); % 解析信号的模即瞬时包络 env env - mean(env); % 去掉包络的直流否则 0 Hz 处压过一切 [f, envP] fft_spectrum(env, fs); envP envP / (max(envP) eps); % 归一化便于多组数据同图对比 endband选哪儿是这套流程里最依赖经验的参数。常用的判断办法是先看原始频谱找那种“鼓起来一片、但没有明显单根谱线”的频段那就是结构共振带。3.1 小节的例子取 2000~5000 Hz在 12k 采样下奈奎斯特是 6000 Hz这个范围不越界。选低了会把转频谐波漏进包络选高了贴近奈奎斯特滤波器过渡带被压扁反而滤不干净。filtfilt而不是filter是因为前者前后各滤一遍相位抵消包络不会整体滞后。做单次分析时这点延迟无所谓但如果你要拿包络去做冲击定位、算时间差相位一偏结论就错了。butter的阶数取 4再往上在filtfilt里容易出现数值不稳定特别是低频截止离 0 很近的时候会看到滤波结果开头一段剧烈振荡。3.4 理论故障特征频率四个公式和一次对谱function ft fault_freq(cfg) % 输入几何参数与转速输出四个理论故障特征频率 fr cfg.rpm / 60; % 转频 Hz r cfg.d_ball / cfg.D_pitch * cosd(cfg.contact);% 直径比接触角为 0 时为 d/D n cfg.n_ball; ft.fr fr; ft.bpfo n/2 * fr * (1 - r); % 外圈通过频率 ft.bpfi n/2 * fr * (1 r); % 内圈通过频率 ft.bsf cfg.D_pitch/(2*cfg.d_ball) * fr * (1 - r^2); % 滚动体自转频率 ft.ftf fr/2 * (1 - r); % 保持架频率 end这四个式子的物理含义不难理解滚动体经过外圈某一点每转一圈经过n/2 * (1-r)次乘上转频就是每秒通过次数。内圈因为滚动体还要“追”着内圈滚道跑通过次数更多所以是1r。滚动体自转最快因为它既绕轴公转又绕自身轴自转。保持架最慢只有转频的一半左右。对谱时不要死盯着一个点。正确的找法是在bpfo ± 1%范围内搜包络谱最大值再检查2×bpfo、3×bpfo处是否也有峰以及峰的两侧有没有fr间隔的边频。谐波和边频同时成立才叫“对上了”。只有一根孤零零的峰很可能是电磁噪声或者某阶共振别急着下结论。4. Haar 小波去噪与特征提取从一条波形到一行特征向量去噪这一步在整条链路里的位置有点微妙。理论上包络解调本身就有抑制低频干扰的作用但真实数据里冲击往往被宽带噪声埋得只剩一点点直接做包络包络谱上那根特征频率的峰可能只比底噪高两三倍靠肉眼判断就变成玄学了。小波去噪把冲击保留、把高斯白噪声压下去峰凸显出来判断就有了底气。4.1 为什么这个场景下 haar 够用小波族很多db4、sym8、coif5 各有拥趸为什么这份代码默认用 haar原因不玄乎haar 是最简单的正交小波滤波器长度只有 2分解重构没有重叠带来的边界拖尾处理几千点一段的短信号时边界效应最小。它的频域过渡带确实粗糙频带划分不锐利但你做的是“把噪声压掉、把冲击留住”不是做精细的频带分离粗糙一点无所谓。另一个好处是可解释性强。haar 的一层分解就是把序列两两平均近似系数和两两求差细节系数人脑能跟上。调试时你可以直接把第一层细节系数画出来冲击在哪儿一目了然不像 db4 那样一堆卷积系数看得头晕。当然如果你的信号里冲击很缓、频率很低haar 的阶跃特性会显得不够平滑那时候换 db4 更合适——改一行字符串的事。4.2 分解、阈值、重构三步实现function [xd, c, l] haar_denoise(x, level, mode) % haar 小波 level 层分解对细节系数做阈值收缩后重构 if nargin 3, mode s; end % s 软阈值h 硬阈值 x x(:).; % 1. 长度截断到 2^level 的整数倍避免 waverec 后长度对不上 n 2^level * floor(numel(x) / 2^level); x x(1:n); % 2. 分解 [c, l] wavedec(x, level, haar); % 3. 用第一层细节系数鲁棒估计噪声标准差 d1 c(end - l(end) 1 : end); sigma median(abs(d1)) / 0.6745; thr sigma * sqrt(2 * log(n)); % 通用阈值VisuShrink % 4. 只对细节系数收缩近似系数原样保留 idx l(1) 1 : numel(c); c(idx) wthresh(c(idx), mode, thr); % 5. 重构 xd waverec(c, l, haar); xd xd(1:numel(x)); end第 3 步的median(abs(d))/0.6745是这套方法的核心。高斯白噪声的绝对值中位数和标准差之间有固定比例关系用中位数估标准差不会被少数几个大幅值的冲击带偏——这正是我们要的因为冲击恰恰是要保留的信号。要是用std(d1)去估冲击一大估出来的 sigma 跟着变大阈值水涨船高最后把冲击一起削平了这就是典型的“去噪去过头”。第 4 步的l(1)1 : numel(c)是把近似系数跳过、只动细节系数。近似系数里装的是信号的主体轮廓和低频成分动了它波形就变形。l是wavedec返回的长度向量第一个元素是最高层近似系数的长度后面的依次是各层细节系数长度顺序是从高层到低层所以第一层细节系数在数组末尾。关于阈值模式软阈值把超过阈值的系数整体收缩重构出来更平滑但冲击的幅值会被削掉一块硬阈值只把小于阈值的置零幅值保得住代价是可能引入一点振铃。诊断场景下我一般先用软阈值看到冲击被压得太扁再换硬阈值。这个切换只要把mode参数从s改成h。4.3 十个特征指标feat_extract.m 里的时域与频域量只看谱形做判断换个工况就容易翻车。稳妥的做法是再抽一组数值特征既能做阈值判断也能直接喂给分类器。function feat feat_extract(x, fs) % 时域 9 个 频域 1 个返回结构体方便 struct2table x x(:). - mean(x); ax abs(x); feat.rms sqrt(mean(x.^2)); % 有效值反映整体能量 feat.peak max(ax); % 峰值 feat.p2p max(x) - min(x); % 峰峰值 feat.kurt kurtosis(x); % 峭度对冲击最敏感 feat.skew skewness(x); % 偏度反映波形不对称 feat.crest feat.peak / feat.rms; % 峰值因子 feat.shape feat.rms / mean(ax); % 波形因子 feat.impulse feat.peak / mean(ax); % 脉冲因子 feat.margin feat.peak / mean(sqrt(ax)).^2; % 裕度因子 feat.sk spectral_kurtosis(x, fs); % 谱峭度指示共振带位置 end function sk spectral_kurtosis(x, fs) % 用短时傅里叶变换近似算谱峭度取最大值 win round(0.01 * fs); hop round(win/2); [S, ~] spectrogram(x, hann(win), win-hop, 256, fs); P abs(S).^2 eps; sk mean(P.^2, 2) ./ (mean(P,2).^2) - 2; % 复高斯噪声理论值为 0 sk max(sk); end峭度是这个列表里最值得盯的一个。正常轴承的振动接近高斯分布峭度在 3 附近一旦出现局部损伤冲击把分布的尾巴拉长峭度能飙到 5 以上甚至十几。它的好处是不依赖转速和轴承几何参数坏处是对偶然的野点也敏感所以一般配合峰值因子一起看。谱峭度则是帮你选带通范围的工具。它对每个频点算“这个频带的信号有多像冲击”值最大的频段就是共振带所在。有了它3.3 节的band就不用凭经验试了直接取谱峭度曲线的峰值位置往两边扩一段。这也是我在实际项目里最常用的一个组合先谱峭度定带再包络解调对谱。裕度因子peak / mean(sqrt(ax))^2分母上那个平方容易写错写成mean(sqrt(ax))^2和(mean(sqrt(ax)))^2是一回事但如果你写成mean(ax)就退化成了峰值因子。这类细节不影响程序跑只影响特征的可比性改的时候对着公式核一遍。4.4 批量跑数据与结果落表单条信号看波形是演示真正干活是整目录批处理。batch_run.m的结构大概是files dir(fullfile(data,*.mat)); rows struct(file,{},rms,{},kurt,{},crest,{}, ... bpfo_ratio,{},bpfi_ratio,{},label,{}); for k 1:numel(files) s load_signal(fullfile(data, files(k).name)); xd haar_denoise(s.vib., cfg.level); [fe, eP] envelope_demod(xd, cfg.fs, cfg.band); ft fault_freq(cfg); % 特征频率附近 ±2% 的峰值与全局中位数之比作为“有没有峰”的量化 m1 fe ft.bpfo*0.98 fe ft.bpfo*1.02; m2 fe ft.bpfi*0.98 fe ft.bpfi*1.02; rows(k).file files(k).name; rows(k).rms rms(xd); rows(k).kurt kurtosis(xd); rows(k).crest max(abs(xd)) / rms(xd); rows(k).bpfo_ratio max(eP(m1)) / median(eP); rows(k).bpfi_ratio max(eP(m2)) / median(eP); rows(k).label regexp(files(k).name,^[a-z],match,once); end T struct2table(rows); writetable(T, fullfile(result,features.csv)); disp(T);bpfo_ratio这个量比“看有没有峰”靠谱得多。分母用包络谱的中位数当底噪基准分子用特征频率附近的峰值比值超过 5 基本可以认为是真峰。这样即使换了设备、换了传感器灵敏度判据依然成立不用每次重新调阈值。label从文件名前缀抠出来是为了后面做有监督的验证命名规范在这里的价值就体现出来了——文件名乱起这批数据基本就废了一半。5. 避坑与排查五个让诊断结论跑偏的细节代码能跑通和结论能信是两件事。下面这几条是实际调试里反复撞上的按“现象—原因—解决”记下来省得下次再从头查一遍。5.1 包络谱上只有 0 Hz 一根巨峰后面全是平的现象画出包络谱第一根柱子顶到天上纵轴一自动缩放100 Hz 往后的部分被压成一条贴近横轴的直线什么都看不见。原因包络信号的均值没有被去掉。abs(hilbert(x))出来的包络恒为正均值是一个不小的直流分量做 FFT 后全落在 0 Hz 那一格幅值比交流成分大一到两个数量级。MATLAB 自动缩放纵轴其余部分自然被压没了。解决包络之后必须再减一次均值env env - mean(env)。如果减完还是有个矮一点的 0 Hz 峰那就是带通滤波器的低频残留把band的下限往上提个几百赫兹或者直接画图时把频率轴从 5 Hz 开始。5.2 换一台设备所有特征频率都对不上现象同样的代码、同样的轴承型号在 A 设备上峰正好落在理论值上搬到 B 设备峰偏移了 4~6 Hz 甚至更多怎么调都对不齐。原因一是转速用了铭牌额定值。带载之后实际转速可能比额定低几十转转频一变四个特征频率按比例全偏。二是滚动体存在打滑实际通过频率通常比理论值低 1%~2%这是滚动轴承的固有特性不是误差。解决转速用转速计实测别信铭牌对谱窗口从 ±1% 放宽到 ±3%同时把理论值乘 0.98 的修正系数再搜一遍。如果两个位置都有峰那基本可以确认是打滑导致的偏移按实测峰为准反推转速更可靠。5.3 waverec 重构后长度和原信号对不上画图报维度错误现象haar_denoise返回的xd和输入x长度差一两个点后面做减法画对比图时报“矩阵维度必须一致”。原因小波分解要求信号长度是 2^level 的整数倍wavedec内部会做边界延拓waverec出来的长度等于延拓后的长度不保证和原始长度相同。层数越高差得越多。解决进函数前先截断n 2^level * floor(numel(x)/2^level); x x(1:n);出函数后再保险一次xd xd(1:numel(x));。如果原始信号本来就短于 2^level比如只有 512 点却设了 10 层分解那floor直接得 0整段信号被扔掉——这种情况要么降层数要么补零别硬上。5.4 去噪之后冲击没了峭度反而比去噪前更低现象加了一层 haar 去噪波形确实光滑了但峭度从 6.2 掉到 3.1包络谱上的特征峰也跟着矮下去等于把有用的东西一起抹了。原因典型的是层数给多了或者阈值估偏了。层数越高被收缩的细节系数越多冲击这种宽带的瞬态成分在高层细节里也有分布一并被压掉。另外如果用std而不是中位数去估噪声冲击会污染 sigma阈值被抬得过高。解决层数从 3 层起步试看峭度保留情况再往上加一般 4~5 层封顶噪声估计坚持用median(abs(d1))/0.6745软阈值压得太狠就换硬阈值。判断标准很简单——去噪后峭度至少要保持去噪前的 80%掉太多说明方法用过头了。5.5 分类准确率 99%换一批数据就崩到 60%现象用自己采的数据训练一个分类器交叉验证准确率漂亮得不像话换个工况、换台设备重新测准确率断崖式下跌。原因同一批数据是按段切分的相邻段之间重叠 50%训练集和测试集里存在大量高度相关的样本等于把答案漏进了考卷。再加上特征没有归一化rms动辄零点几、kurt是几、bpfo_ratio能到几十量纲差两个数量级距离-based 的分类器会被大数值特征完全主导。解决切分时按“文件级”而不是“段级”划分训练测试集同一个.mat里的段只能进同一个集合特征统一做 z-score 标准化报告结果时给出混淆矩阵而不是单一准确率。这三件事做了准确率往往会掉十几个点但那才是真实水平。6. 二次开发换轴承型号、换采样率、接分类器的三条改造路径拿到这份代码多数人的下一步都是“用到自己的设备上”。改动其实集中在三个地方而且都不需要动主流程。第一处是参数表化。把cfg里的轴承几何参数抽成一张表按型号查表取值function cfg get_cfg(model) % 常见型号的几何参数查表找不到就报错别默默用默认值 T readtable(bearing_params.csv); % 列model,n,d,D,contact row T(strcmpi(T.model, model), :); assert(~isempty(row), 型号 %s 不在参数表中, model); cfg.fs 12000; cfg.rpm 1730; cfg.n_ball row.n; cfg.d_ball row.d*1e-3; cfg.D_pitch row.D*1e-3; cfg.contact row.contact; end第二处是采样率与频带。fs一换band必须跟着换因为共振带是按绝对频率选的。稳妥的做法是用 4.3 节的谱峭度自动定带把band从手工输入改成[f_peak*0.6, f_peak*1.6]跟着数据走而不是跟着经验走。同时确认新采样率下的奈奎斯特频率高于band上限否则butter会直接报归一化频率越界。第三处是接分类器。特征表已经落成 CSV直接喂给分类模型就行用fitcecoc做多分类、按文件级划分训练测试集、输出混淆矩阵这三步是验证的底线。我一般会把bpfo_ratio、bpfi_ratio、峭度、谱峭度这四个作为主特征其余作参考——特征不是越多越好冗余特征在样本量小的时候只会让模型学到噪声。从那次 99% 准确率翻车以后我每次做诊断验证都强制走一遍“按文件划分 特征标准化 看混淆矩阵”这三步少一步都不敢往报告里写。这套代码的价值不在算法多新而在于它把从原始波形到判决结论的每一环都摆在了明面上你想改哪一环改完立刻能在图上看到后果。希望帮到你。本文还有配套的精品资源点击获取