
简介面向机械故障诊断与信号处理领域的工程师、研究生及高年级本科生这份Matlab源码包聚焦包络谱轴承故障诊断提供一套可直接运行的完整分析流程。资源共8个文件压缩包仅1.14MB包含7个.m脚本和1个.mat实测振动数据脚本按功能划分既有希尔伯特变换、包络谱计算、FFT频谱分析等通用函数也有内圈故障、外圈故障等专题分析代码便于按需调用与二次开发。通过加载.mat数据并运行主脚本即可自动完成滤波去噪、希尔伯特解调、包络谱绘制以及滚动体、内圈、外圈故障特征频率的提取与对比直观展现正常与故障状态的频谱差异有助于深入理解包络谱诊断机理。同时脚本中保留了关键参数与注释通过修改转速、采样率等参数还可适配不同工况信号为科研验证、课程设计或工程排故提供良好起点。目前已有2187人学习下载适合具备基础Matlab操作能力、希望快速复现包络谱诊断方法的读者使用。1. 包络谱轴承早期故障诊断里成本最低的一步解调轴承外圈剥落一小块振动总能量几乎不变直接看频谱只能见到载波附近几根细边带很容易漏判。把同一段信号做 Hilbert 变换、取模、再对包络做 FFT冲击重复频率会变成低频段清晰的主峰这就是包络谱也叫共振解调。我处理电机、泵、轧机振动数据时怀疑轴承早期故障的第一动作就是取包络谱。这篇内容适合设备监测工程师也适合想用 MATLAB 验证信号处理流程的研究生。它不需要深度学习模型一个 butter 滤波器加 hilbert 函数就能跑通但带通频段选错包络谱上会什么都看不见。下面把特征频率、参数设置和误判点一次讲清。2. 轴承故障频率模型与包络谱的理论边界2.1 外圈、内圈特征频率公式与典型轴承参数轴承局部缺陷与滚道接触会产生周期性冲击冲击重复频率由缺陷所在位置决定。设 Z 为滚动体个数fr 为转轴转频d 为滚动体直径D 为轴承节圆直径α 为接触角四个特征频率的简化表达式如下外圈故障频率 BPFO Z/2 · fr · (1 − d/D · cosα)内圈故障频率 BPFI Z/2 · fr · (1 d/D · cosα)滚动体故障频率 BSF D/(2d) · fr · (1 − (d/D · cosα)²)保持架故障频率 FTF 1/2 · fr · (1 − d/D · cosα)内圈频率系数大于外圈是因为内圈随轴旋转单个缺陷在承载区内经过的次数更多。以公开数据集中最常见的 6205-2RS 深沟球轴承为例Z9节圆直径约 39.04 mm滚动体直径约 7.94 mm接触角近似 0°上面四个频率可以直接换算成转频倍数故障位置频率代号fr 倍数6205-2RS 参考外圈BPFO3.585内圈BPFI5.415滚动体BSF2.323保持架FTF0.398这套系数在转速变化时保持比例关系所以现场诊断通常把包络谱横轴归一化到 fr 的倍数而不是绝对频率。换轴承型号时不能照抄系数必须把实测的 d、D、α 代回公式重算接触角一般取轴承手册标称值多数深沟球轴承在 0° 到 15° 之间。2.2 为什么直接 FFT 看不到早期故障轴承振动可以近似看成调幅信号冲击以重复频率 f0 激励起轴承座某个结构共振频率 fc振动波形的表达近似为 x(t) A(1 m·cos 2πf0t)·cos 2πfct。直接对这个信号做 FFT能量集中在 fc 和 fc±f0 的边带上。早期故障的调制深度 m 很小边带只比噪声高几个 dB而且 fc 经常落在 2 kHz 以上的结构共振区不在常规频谱分析的重点频段内所以很难确认峰值。包络谱换了个思路先用带通滤波器把载波所在的共振频带单独滤出来再用 Hilbert 变换得到瞬时包络 A(1 m·cos2πf0t)对包络做 FFT 后f0 及其谐波变成低频段的主峰信噪比远高于直接频谱。这里有一个容易忽略的点包络检波本身是强非线性操作对 |cos(2πf0t)| 做频谱天然含有 f0 的 2 倍、3 倍整倍谐波所以包络谱里出现多倍频是正常现象不能把 2×BPFO 误判成另一个独立故障。提示对包络做 FFT 前必须去直流否则 0 Hz 的巨大分量会把低频段谱线压得看不清。2.3 包络谱的两个关键参数共振频带与频率分辨率包络谱质量由两个参数决定。第一是带通滤波范围它必须覆盖轴承冲击激励起的结构共振峰通常在 2 kHz 到 15 kHz 之间具体数值跟轴承座刚度、传感器安装方式高度相关需要看频谱图或做谱峭度搜索来确定。第二是 FFT 频率分辨率 Δf fs/N想区分间隔为转频 fr 的边带数据时长至少要有 1/fr一般取 30 到 50 个转轴周期。举个例子转频 30 Hz 时要分辨 1 Hz 以内的谱线数据长度至少 1 秒要做到 30 圈以上统计稳定采集时长要 5 秒以上。如果现场采集系统的采样率只有 2 kHz共振频带根本采不进来包络谱对早期故障基本无效这是采集方案阶段就要确认的硬条件。3. MATLAB 实现Hilbert 变换、带通滤波与包络谱脚本3.1 从 .mat 读取振动数据与预处理这个源码包里振动数据存放在 zczdI7.mat 中从文件名看是某个测点的振动采集结果。读取时先用 whos 确认变量名不要直接把变量名写死load(zczdI7.mat); whos % 查看数据变量名和维数 x data(:); % 统一转为列向量避免行列方向问题 fs 12000; % 采样率按采集系统实际设置修改 x x - mean(x); % 去直流分量load 之后的具体变量名以文件实际内容为准常见命名有 data、acc、vib 三种。强制转成列向量是因为 hilbert 和 fft 对行列的处理结果一致但后面滤波、拼接时行向量容易出维度错误。去直流是必要预处理否则 FFT 后 0 Hz 处会堆一个无关的大分量压缩整个纵轴动态范围。如果现场数据不是 .mat 而是 CSV用 readmatrix(data.csv) 读入后走完全相同的流程。3.2 带通滤波零相位滤波器锁定共振带Hilbert 变换前必须先带通滤波。源码包里常见组合是 fir1 配 filterfir1 的阶数需要根据 fs 调整固定写 100 阶时不同采样率下过渡带宽度差异很大filter 是因果滤波会产生与阶数相关的群延迟对包络的瞬时幅值有不可忽略的偏移。我做离线分析时更习惯用 butter 配合 filtfiltfc 5000; % 共振带中心频率需根据谱图调整 bw 2000; % 共振带带宽 fL (fc - bw/2) / (fs/2); fH (fc bw/2) / (fs/2); [b, a] butter(4, [fL fH], bandpass); x_f filtfilt(b, a, x);butter 的第二个参数必须是归一化频率除以奈奎斯特频率 fs/2 是最容易漏掉的细节fs 必须与采集卡实际采样率一致。filtfilt 对信号先正向滤波再反向滤波零相位失真适合离线分析4 阶巴特沃斯在阻带衰减和计算量之间比较平衡。带宽设太窄会把冲击信号的边带能量削掉设太宽会引入邻近齿轮啮合分量一般先按 2 kHz 带宽扫一遍找到清晰峰值后再收窄。3.3 Hilbert 变换、包络计算与幅值修正MATLAB 的 hilbert(x_f) 返回解析信号实部是原信号虚部是原信号的希尔伯特变换取绝对值就得到包络。对包络做 FFT 时要自己完成去直流、单边谱截取和幅值修正env abs(hilbert(x_f)); % 包络信号 env env - mean(env); % 去直流防止 0 Hz 大分量 N length(env); Y fft(env); P abs(Y(1:ceil(N/2))) * 2 / N; % 单边幅值谱 f_ax (0:ceil(N/2)-1) * fs / N; % 频率轴 plot(f_ax, P); xlim([0 500]); grid on; xlabel(频率 (Hz)); ylabel(包络谱幅值);abs(Y(1:ceil(N/2))) * 2 / N 完成单边谱修正FFT 结果在正负频率对称分布取一半后乘 2 再除以 N峰值幅度才和时域信号真实幅度一致。频率轴每格 fs/NN 越大谱线越密。xlim 设为 500 Hz 是因为包络谱只关心低频调制分量轴承故障特征频率通常在 500 Hz 以下高频部分画出来反而干扰判断。3.4 一个可复用的包络谱脚本骨架把上面几步串联起来就是源码包里 shili.m 的组织方式。给一个不依赖具体文件名的完整骨架方便替换数据直接跑load(zczdI7.mat); x data(:); fs 12000; x x - mean(x); [b, a] butter(4, [4000 6000]/(fs/2), bandpass); x_h filtfilt(b, a, x); env abs(hilbert(x_h)); env env - mean(env); NFFT 2^nextpow2(length(env)); % 补零到 2 的幂加速 FFT P abs(fft(env, NFFT)); P P(1:NFFT/2) * 2 / length(env); f (0:NFFT/2-1) * fs / NFFT; findpeaks(P, f, MinPeakHeight, max(P)*0.25, MinPeakDistance, 2);NFFT 取 2 的幂只影响运算速度和谱线插值密度不会改变峰值位置。findpeaks 的 MinPeakHeight 设为最大峰值的四分之一避免把噪声小峰当故障特征MinPeakDistance 至少设 2 Hz防止同一个谱峰的两根相邻谱线被识别成两个峰。源码里的 neiquanguzhang.m 和 waiquanguzhang.m 内部结构与此一致调试时分别在滤波器后、包络后打印一段数据就能快速定位问题在滤波还是 FFT。4. 源码结构拆解内圈与外圈故障脚本怎么协作4.1 源码文件分工与调用关系压缩包里 7 个文件的分工比较清晰shili.m 是主入口负责加载 zczdI7.matneiquanguzhang.m 处理内圈故障waiquanguzhang.m 处理外圈故障两者都会调用 hilbertbianhuan.m 完成希尔伯特变换用 baoluopu.m 绘制包络谱用 shiyufft.m 生成直接 FFT 的对照图emdfj.m 是可选的 EMD 预处理模块。整体依赖关系是 shili.m 按故障位置分流到内圈或外圈脚本两个脚本内部共享希尔伯特变换和包络谱绘图函数shiyufft.m 的结果用于证明“直接频谱看不出故障而包络谱能看出”。我一般会把两个故障脚本里相同的带通和峰值提取部分抽成公共函数只保留特征频率判据的差异这样改一次滤波参数内圈外圈同时生效避免“内圈改了外圈忘了改”的低级错误。4.2 内圈与外圈故障的包络谱判据内圈故障的包络谱特征不止一个 BPFI 主峰那么简单。内圈随轴旋转缺陷相对承载区的位置周期性变化冲击幅度被转频二次调制所以 BPFI 两侧会出现 BPFI±fr、BPFI±2fr 的边带族。外圈固定在轴承座上缺陷相对承载区位置固定包络谱通常是干净的 BPFO 基频加高次谐波没有明显的 ±fr 边带。判读时看三点主峰频率除以转频得到 5.4±0.1 附近判定内圈3.6±0.1 附近判定外圈主峰两侧有间隔等于 fr 的对称边带进一步确认为内圈故障如果高频段出现一对距离很近的对称边带优先怀疑齿轮啮合而不是立刻下轴承结论。峰值搜索不能只看最高峰。假设转频 29.5 HzBPFO 约 105.8 HzBPFI 约 159.7 Hz电网 50 Hz 的 3 次谐波 150 Hz 很容易与 BPFI 混淆。区分的办法只有一个变转速验证。把转速从 1750 r/min 调到 1450 r/min故障频率按比例漂移电源相关频率纹丝不动在 MATLAB 里对两组数据分别取包络谱横轴按 fr 比例缩放后叠加一锤定音。4.3 EMD 预处理的定位与坑emdfj.m 在这个流程里不是必须的。EMD 把信号分解成多个 IMF对冲击型信号可以分离出与故障相关的高频固有模态再对选定的 IMF 做包络谱能在共振带重叠时改善频谱清晰度。但 EMD 有两个工程陷阱端点效应会使数据两端产生虚假振荡包络谱两端出现假峰模态混叠会让一个 IMF 混入多个时间尺度故障频率被平均掉后用包络谱也搜不出来。我的使用边界是直接带通包络谱能看清峰值时完全不用 EMD只有现场多个冲击源叠加、包络谱一片糊时才先跑 emdfj.m取与原始信号相关系数最高的 IMF 再分析。对 IMF 做包络谱前把数据两端各截掉 10%避开端点效应污染区imf emdfj(x_h); % 常见输出为 IMF 矩阵每列一个模态 env_imf abs(hilbert(imf(:,1))); env_imf env_imf(round(0.1*end):round(0.9*end)); % 截掉两端IMF 选第几列不要固定用两个指标排峭度越大说明冲击成分占比越高与原始信号相关系数不能太低否则选到的是噪声模态。如果 emdfj.m 输出的是元胞数组把 imf(:,1) 改成 imf{1}运行前先看函数返回值类型这类接口不统一问题在开源脚本里很常见。5. 共振频带自动搜索与现场诊断的几个快决策5.1 用谱峭度粗搜共振频带的降级实现包络谱的带通中心频率靠人眼在频谱上找既慢又不稳定。谱峭度工具包能自动定位冲击能量最集中的频带但不是 MATLAB 内建函数装起来麻烦。这里给一个五分钟能写完的粗搜版本把 0 到 fs/2 分成若干带通区间逐段滤波、取包络、算峭度峭度最大的频带就是冲击激励最明显的共振带。band 1000; step 500; f_edges 0:step:fs/2-band; kurt_values zeros(size(f_edges)); for k 1:numel(f_edges) fL f_edges(k) / (fs/2); fH (f_edges(k)band) / (fs/2); [b, a] butter(4, [fL fH], bandpass); yf filtfilt(b, a, x); kurt_values(k) kurtosis(abs(hilbert(yf))); end [~, idx] max(kurt_values); fc_best f_edges(idx) band/2;band 是搜索带宽step 是步进step 小于 band 时相邻频带重叠搜索更平滑但计算量增大。对 10 秒 12 kHz 采样的数据这个循环在普通桌面 CPU 上不到 30 秒。峭度对冲击极其敏感早期故障的包络峭度明显高于正常状态这也是为什么用滤波输出算峭度而不是对原始信号直接算。5.2 现场诊断的几条经验转频用频谱峰值搜索获得不要用铭牌转速异步电机满载和轻载的转差不同用 30 Hz 附近最大峰对应的频率作为 fr 更可靠。数据长度少于 20 个转轴周期时边带糊成一片采集时优先保证时长采样率满足共振带覆盖即可。包络谱峰值幅值不能直接当故障严重度要和同一测点同一工况的历史基线比上涨 3 倍以上才有明确恶化信号。数据文件要脱离 MATLAB 查看时用 Python 的 scipy.io.loadmat 读取keys() 查看变量名不需要打开 MATLAB 就能确认数据结构。5.3 把包络谱横轴归一到转频倍数最后说一个现场对比的好用做法把频率轴除以转频 fr包络谱变成阶次谱内圈故障固定出现在 5.415 阶外圈固定 3.585 阶转速变化时峰位不漂移。这样不同转速下的历史包络谱可以直接叠加比较不用每次重新标注理论频率。order_axis f_ax / fr; % fr 从包络谱低频段转频峰获得 plot(order_axis, P); xlim([0 10]);这个方法只在转速波动小于 1% 时成立转速大幅波动时阶次峰会展宽甚至折叠需要先做角度域重采样再回到帧内包络分析MATLAB 里用 resample 按角度间隔重采样即可实现。本文还有配套的精品资源点击获取