包络谱分析：轴承故障诊断原理与MATLAB实现方法

简介：面向机械故障诊断与信号处理领域的工程师、研究生及高年级本科生，这份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 深沟球轴承为例，Z=9，节圆直径约 39.04 mm，滚动体直径约 7.94 mm，接触角近似 0°，上面四个频率可以直接换算成转频倍数：

这套系数在转速变化时保持比例关系，所以现场诊断通常把包络谱横轴归一化到 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 配 filter，fir1 的阶数需要根据 fs 调整，固定写 100 阶时，不同采样率下过渡带宽度差异很大；filter 是因果滤波，会产生与阶数相关的群延迟，对包络的瞬时幅值有不可忽略的偏移。我做离线分析时更习惯用 butter 配合 filtfilt：

fc = 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/N，N 越大谱线越密。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.mat；neiquanguzhang.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 Hz，BPFO 约 105.8 Hz，BPFI 约 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 按角度间隔重采样即可实现。

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