简介基于EMD的包络谱故障诊断MATLAB程序实例是一套面向机械设备状态监测与故障诊断学习者的完整示例工程用于对非线性、非平稳振动信号进行经验模态分解和包络解调进而提取滚动轴承等设备的故障特征频率。压缩包共包含15个文件RAR格式仅1.33MB其中3个.m脚本实现EMD分解、包络计算和频谱分析流程1个.mat文件提供实验振动数据7张.jpg为原始信号时域波形、包络谱、IMF频谱等关键结果图另有txt文本和docx文档辅助说明算法原理与使用步骤便于对照学习。已有909人浏览学习适合正在学习Hilbert包络解调、EMD方法或应用MATLAB信号处理工具箱的研究生与工程师。通过运行这些代码可以直观理解从原始信号到IMF分解、再到包络谱故障频率识别的完整链条并可直接替换自己的数据开展诊断实验是一份实用的入门与参考代码包。1. EMD 包络谱到底在解决哪一类故障诊断问题振动诊断里最常见的一类挫败感是频谱图上转频及其谐波清清楚楚可就是找不到故障频率。设备该响的冲击在时域波形里也看得到但直接对原始信号做 FFT冲击被连续旋转分量和随机噪声压住能量散在很宽的频带上故障特征始终出不来。这个问题在滚动轴承早期故障、齿轮局部剥落、往复机械阀座泄漏上反反复复出现。包络解调的思路是把“冲击反复出现”这件事从“高频载波”上剥离出来先做 Hilbert 变换得到解析信号取瞬时幅值作为包络再对包络做频谱分析故障特征频率就会以包络谱上尖锐谱线的形式暴露出来了。而 EMD经验模态分解在这个流程里负责把信号拆成若干个本征模态函数IMF从中挑出冲击主导的分量再进行包络谱分析能显著减少转频、倍频和噪声的干扰。本文面向需要用 MATLAB 做旋转机械故障诊断的工程师和科研人员落地到一套可以复现的程序先构造一段带外圈故障冲击的仿真信号再按“EMD 分解→挑选 IMF→Hilbert 包络→包络谱→定位特征频率”的流程跑通然后讨论采样率、带宽、IMF 筛选和数据长度这些真正决定诊断成败的参数末尾给出一个可直接复用的函数封装。读完你会发现包络谱本身并不神秘难点全在“选对 IMF”和“确认谱峰不是假的”这两件事上。2. EMD 和包络谱的原理为什么包络谱比直接频谱更靠近故障特征2.1 EMD 的筛分过程和 IMF 需要满足的两个条件EMD 的核心假设是一段复杂信号可以分解为若干个频率从高到低排列的本征模态函数IMF加一个残余趋势项。实测振动信号往往是非平稳、非线性的FFT 对这种信号的处理是从全局频率域逼近而 EMD 是自适应地把信号按局部时间尺度逐层剥离。每剥离出一层 IMF相当于从原始信号里抹去一个“平均局部周期”分量剩下的残差继续进入下一轮筛分。筛分过程大致是对当前信号找出全部局部极大值和极小值用三次样条分别包络上下极值点得到上下包络线取均值得到 m1再用原始信号减去 m1 得到候选分量 h1。如果 h1 不满足 IMF 条件就重复这个过程满足条件以后h1 就是第一个 IMF从原始信号中分离对残余信号继续分解直到残余项成为一个单调函数或幅度小于阈值。工程上判断一个分量算不算 IMF看两个条件第一整个数据段内极值点数与过零次数相等或最多相差 1第二对任意时刻上下包络的均值必须接近 0。这两个条件保证了 IMF 是一个窄带分量也就是说它的瞬时频率有意义瞬时幅值能够真正反映调制信息。包络谱分析的全部可靠性都建立在这条性质上只有窄带分量的包络才有明确的调制含义宽带分量取包络会得到一堆混叠的频率尖峰根本没法解释。在 MATLAB 里跑emd(x, MaxNumIMFs, 8)会做一件事最多抽出 8 层 IMF每一层按频率从高到低排列通常第 4 层以后就已经接近周期平滑项。注意EMD 不是变换没有固定基函数每次分解的结果受停止条件、极值样条插值方式和端点延拓方式影响所以别指望两次分解得到逐位相同的输出诊断关注的是谱峰位置而非波形逐点一致。2.2 包络解调从调制信号到包络谱只差两步变换轴承局部故障形成的冲击链本质上是高频结构共振被低频周期冲击所调制。传感器拾取到的信号里高频载波中心频率取决于轴承座、端盖、传感器的结构固有频率而调制频率才是滚动体通过故障位置的重复频率也就是我们关心的故障特征频率。包络解调的过程可以分成两步先用 Hilbert 变换构造解析信号对原始信号 z(t) 做 Hilbert 变换得到其正交分量再合成复解析信号对这个解析信号取模得到的实正序列就是包络。代码上就是abs(hilbert(x))一行事。第二步对包络做一次 FFT包络频谱上的谱峰位置就直接对应调制频率也就是故障特征频率及其谐波。这里能解释一个常见疑问为什么直接看原始信号频谱找故障频率很困难。原始频谱里故障特征频率作为“调制频率”会以边带形式出现在高频共振峰两侧幅值很小而轴转频、齿轮啮合频率等大能量成分占据主导。包络谱相当于做了一次“解调”把调制频率从载波频率搬到低频段信噪比一下子就变了。实际算包络谱时通常还要对包络序列减去直流分量再乘窗函数否则低频段有一个巨大的 0 Hz 分量会掩盖几百赫兹内的真实故障谱峰。2.3 EMD 家族怎么选EMD、EEMD、CEEMDAN 在包络诊断里的分工经典 EMD 有个老毛病叫模态混叠一个物理分量被切成两段分别出现在相邻两个 IMF 里或者两个不同尺度的分量叠加进同一层 IMF。冲击信号特别容易触发模态混叠因为随机冲击会引入非均匀分布的极值点。为了压制这个问题出现了加入白噪声的集合平均方案也就是 EEMD 和 CEEMDAN。日常诊断工作中我一般这样选方法核心改进包络诊断适用场景计算代价EMD直接筛分信号平稳、转速稳定、故障特征明确最小适合批量试算EEMD添加有限幅值白噪声并多次平均冲击成分弱、强噪声背景下找弱故障集合次数越高越贵CEEMDAN每级加入自适应噪声消除残余噪声希望 EMD 分解结果稳定、可复现比 EEMD 更重如果故障冲击很强在时域波形里已经肉眼可见那就是经典 EMD 的主场。如果现场信号噪声很大、故障早期建议先用 EEMD 或 CEEMDAN 跑一轮挑出最接近冲击形态的 IMF 再做包络谱。EEMD 里有三个参数值得固定下来噪声标准差取原始信号标准差的 0.1 到 0.3 倍集合次数至少 50 次白噪声用不相干的独立随机序列生成。记住加了白噪声的 EEMD 不存在“唯一正确分解”它的价值是让统计结果稳定到可以复现。3. 用 MATLAB 从仿真信号跑通 EMD 包络谱的最小实例3.1 首先生成一段带外圈故障冲击的仿真振动信号没有现场数据也能验证整个流程正确性的办法是构造仿真信号。我用一组典型参数生成 1 秒数据采样率 25600 Hz轴转频 30 Hz外圈故障特征频率 96.3 Hz。冲击信号用衰减正弦振荡近似振荡中心频率 3200 Hz 模拟传感器安装位置的结构共振峰再叠加上转频及其二倍频、白噪声模拟真实振动中连续旋转分量和背景噪声并存的情况。fs 25600; % 采样率 T_total 1; % 信号时长 1 秒 N fs * T_total; % 样本点数 t (0:N-1) / fs; % 时间轴 fr 30; % 轴转频 BPFO 96.3; % 外圈故障特征频率 % 生成冲击发生时刻加入 0.5% 量级的随机抖动模拟真实转速波动 n_imp floor(T_total * BPFO); imp_base (0:n_imp-1) * (1 / BPFO) 0.02; imp_times imp_base 0.005 * randn(1, n_imp) / BPFO; fn 3200; % 结构共振频率作为调制载波 zeta 0.04; % 阻尼比决定冲击衰减快慢 x_imp zeros(1, N); for k 1:n_imp idx floor(imp_times(k) * fs) 1; if idx N - 2 L min(N - idx, round(0.15 * fs)); tt (0:L-1) / fs; x_imp(idx:idxL-1) x_imp(idx:idxL-1) ... exp(-zeta * 2 * pi * fn * tt) .* sin(2 * pi * fn * tt); end end % 叠加两个转频谐波和白噪声信噪比约 10 dB x x_imp 0.6 * sin(2*pi*fr*t) 0.3 * sin(2*pi*2*fr*t) 0.2 * randn(1, N);这段代码里前三个变量是流程决策的关键采样率 25600 Hz 是为了覆盖 3200 Hz 共振峰并留出足量频带余量故障特征频率 96.3 Hz 是后续包络谱的检验基准点阻尼比 0.04 决定冲击在时域持续约 3 到 5 个周期。每次冲击用一个 for 循环叠加远比用卷积实现更直观、也更容易修改抖动幅度。信号构造完成后x就是一个振幅在 -2 到 2 之间的随机起伏序列时域上能看到明显冲击链频谱上则表现为 3200 Hz 附近的一团能量。3.2 执行 EMD 并按峭度加相关系数挑出故障 IMF对已经带通滤波或不带通滤波的信号做 EMD输出是一个imf矩阵每行一个分量。仿真心电图式的思维在这里无效分解完之后最要紧的是选哪个 IMF 去算包络谱。两个指标组合使用峭度挑出冲击性最强的分量相关系数保证所选分量与原信号有充分的能量联系。imf emd(x, MaxNumIMFs, 8); % 最多分解 8 层 n_imf size(imf, 1); K zeros(n_imf, 1); % 峭度 C zeros(n_imf, 1); % 与原信号的相关系数 for k 1:n_imf imf_k imf(k, :); m mean(imf_k); s std(imf_k); K(k) mean((imf_k - m).^4) / (s^4) - 3; % 超值峭度 C(k) abs(sum(imf_k .* x) / sqrt(sum(imf_k.^2) * sum(x.^2))); end % 筛选相关系数大于 0.3优先选峭度最大的 cand_idx find(C 0.3); [~, max_k_idx] max(K(cand_idx)); imf_used cand_idx(max_k_idx); fprintf(选用第 %d 个 IMF峭度 %.2f相关系数 %.3f\n, ... imf_used, K(imf_used), C(imf_used));峭度按超值峭度定义计算白噪声分量的超值峭度接近 0而含冲击的 IMF 通常能达到 5 到 20。相关系数取 0.3 作为门槛是为了避免选到那一层幅度很小、只在局部区域与原信号重合的高频残余。实际工作中经常遇到第一个 IMF 峭度也大、相关系数也高但它的包络谱里全是噪声尖峰原因在于它把随机噪声的高频突发也当成冲击来解调了这就是为什么一定要组合第二个指标做筛选。运行这段代码最常被选中的是第 1 或第 2 个 IMF具体峰度值取决于白噪声强度不用追求固定结果。需要提醒的是emd函数在现代版本 MATLAB Signal Processing Toolbox 中直接可用使用前用which(emd)确认当前环境能找得到这个函数。老版本没有自带实现时把第三方 EMD 实现放进搜索路径接口改成imf emd(x)即可后续包络谱代码保持不变。3.3 用 Hilbert 包络谱定位故障特征频率选定 IMF 之后包络谱步骤非常短但对细节敏感。要先用hilbert取解析信号对解析信号取模得到包络序列再减直流、加窗、做 FFT。这一串操作里最容易被省略的一步是减直流不减的情况下 0 Hz 处会耸起一根极高的谱线在绘图时把其他频率谱线全部压成一条水平线。% 对选中的 IMF 做包络解调 env abs(hilbert(imf(imf_used, :))); % 瞬时包络 env env - mean(env); % 去直流 % 加汉宁窗后再补 nfft nfft 8192; win hanning(length(env)); env_win (env .* win) * nfft / sum(win); % 幅度校正 Y fft(env_win, nfft); f (0:nfft-1) / nfft * fs; mag 2 * abs(Y(1:nfft/21)); % 打印 0.5 Hz 500 Hz 范围内前 3 个高峰并和 BPFO 对比 idx_show f 0.5 f 500; [f_peaks, locs] findpeaks(mag(idx_show), SortStr, descend, ... NPeaks, 5); fk_peaks f(idx_show); for i 1:length(f_peaks) fprintf(峰值频率 %.2f Hz幅值 %.4f\n, fk_peaks(locs(i)), f_peaks(i)); end汉宁窗的幅度校正在这里很重要如果不乘上nfft / sum(win)加窗后同一根谱线的幅值会比矩形窗低约一半换窗后看起来像故障变严重了或减弱了容易误导判断。findpeaks的输出参数里locs是相对idx_show的索引因此重建实际频率时要取f(idx_show)再按位置访问这个细节不对的话打印出的频率和图上标注对不上。在这个仿真参数下包络谱上应该在 96.3 Hz 处出现主峰在 192.6 Hz 和 288.9 Hz 处出现递减的谐波。附带说明一点直接对原始信号x做 FFT也能看到 96.3 Hz 这个频率但它的幅值会被背景噪声和调制结构压得极低而经过“EMD 选冲击分量 → 包络解调”之后主峰会高出周围谱线一个量级以上。这就是包络谱的价值它不是生造出新频率而是把已经存在的调制能量集中到一根谱线上。提示如果仿真信号里不加随机抖动也就是冲击间隔严格等间距包络谱会出现一根极窄的主峰加上微小抖动后主峰带宽扩大但峰值仍能锁定在特征频率附近。现实中转速不可能绝对稳定所以诊断时看到主峰附近有一个 1 到 3 Hz 的小带宽是正常现象。4. 把实例改成现场数据要面对的参数与坑4.1 特征频率计算、采样率和数据长度最少取多少仿真里的 BPFO 是直接设置的现场诊断时必须从轴承几何参数和实际转速算出来。四个特征频率公式照着 SKF 或 ISO 的标准写法展开特征频率计算公式说明外圈 BPFO( \frac{Z f_r}{2} (1 - \frac{d}{D} \cos\alpha) )滚动体通过外圈固定损伤位置的频率内圈 BPFI( \frac{Z f_r}{2} (1 \frac{d}{D} \cos\alpha) )损伤随内圈转动通过频率略高于 BPFO滚动体 BSF( \frac{D f_r}{2d} (1 - (\frac{d}{D} \cos\alpha)^2) )滚动体自转一周与内外圈各接触一次保持架 FTF( \frac{f_r}{2} (1 - \frac{d}{D} \cos\alpha) )保持架转速常作为 BPFO 的旁证公式中 Z 是滚动体个数d 是滚动体直径D 是节圆直径α 是接触角f_r 是轴频。算完特征频率后再定采集参数工程上有一个最少条件信号时长至少要覆盖几十个故障周期一般取 20 到 50 个 BPFO 周期以上转速低的设备尤其要注意1 秒数据对 0.5 Hz 的故障周期来说根本不够分辨出特征峰。用频率分辨率公式f_res 1 / T_total反推若想分辨 2 Hz 间隔的边带最少就需要 0.5 秒数据实际要预留 3 到 5 倍余量。采样率的选取不能只看故障特征频率要照顾到冲击激发的高频共振峰。轴承故障冲击的共振峰通常落在 1 kHz 到 8 kHz 之间加速度传感器安装越靠外共振带越高。常见做法是采样率取共振峰上限的 2.56 倍以上16 kHz 或 25.6 kHz 是现场采集器的通用档位。仿真里取 25600 Hz 就是为了覆盖 3200 Hz 共振峰同时给包络解调留下足够带宽。4.2 先带通还是直接 EMD两套流程都要会对仿真信号直接 EMD 也能选出冲击 IMF现场数据里噪声复杂推荐先做一次带通滤波把共振带框出来再做 EMD。这样做有两个收益一是把与故障无关的强低频转频能量滤掉避免 EMD 把转频当作独立分量抽出减少分解层数二是压缩了进入 EMD 的信号带宽模态混叠的概率显著降低。% 带通滤波中心频率根据共振峰位置调整 f_low 2000; f_high 6000; [b, a] butter(4, [f_low, f_high] / (fs / 2), bandpass); x_fil filtfilt(b, a, x);带通范围和共振峰位置怎么确定最简单可靠的办法是给原始信号做功率谱找频谱上有明显能量隆起的频段那就是结构共振带。常见做法是取隆起最高峰频率左右各 30% 到 50% 的频宽。filtfilt是零相位滤波衰减后的信号不会产生相移这对后续包络谱里频率定位至关重要所以不要用普通filter替代。两套流程的选择逻辑如果信号信噪比高、时域冲击看得见直接 EMD 就行少一步滤波能保留更多原始调制细节如果信号背景噪声大先带通后 EMD选 IMF 时也更稳定。也可以反过来对比验证对同一段数据分别跑两套流程如果包络谱主峰都落在同一特征频率附近诊断可信度大幅上升。4.3 模态混叠、端点效应和冲击如何躲进错误 IMF模态混叠最常见的一幕冲击被 EMD 拆成了两段一段出现在第 1 个 IMF另一段出现在第 2 个 IMF。选唯一个 IMF 做包络谱时总觉得谱峰处幅值比预期小一截谐波结构也不规整。解决的办法有两种。第一种简单粗暴把相邻 IMF 的包络谱都画出来人工确认能量分布。第二种更适合工程化用带通预滤波把 2000 到 6000 Hz 之外的内容全部挡在外面让冲击分量只有一个窄带尺度EMD 就不容易把它拆散。端点效应会表现在一个非常具体的迹象上选出的冲击 IMF 在数据段两端出现幅度逐渐增大或呈扫描状的振荡这会直接污染包络谱低频端产生几根假谱峰。处理方式通常是对原始信号做端点延拓后再 EMDEMD 完成后把两端延拓部分裁剪掉。MATLAB 自带实现里对端点有延拓策略但现场数据如果本身很短、冲击间隔又长推荐手动给信号前后各补一段长度为本故障特征频率 5 倍周期的镜像延拓数据。冲击“躲”进错误 IMF 的情况也经常由幅值差异引起强冲击若干次弱冲击若干次弱冲击对应的调制成分可能被分解到相邻 IMF 中。核对方法是比较各 IMF 包络谱在特征频率处的幅值如果两个相邻 IMF 的同频峰值都高于背景就把这两个 IMF 的包络谱做幅值叠加再判断。注意这是幅度叠加不是直接加原始波形原始波形叠加反而会破坏 IMF 的窄带性质。4.4 误诊断的常见来源包络谱里出现特征频率附近有峰真正确认故障之前要过一遍最常见的三种误判来源。第一种是转频谐波直接落入特征频率附近当特征频率恰好接近某个整数倍的转频比如 3 倍频附近时包络谱会出现一个与转频等间距的假峰这与轴承故障无关。核对方法是看主峰左右有没有转频间隔的边带边带间距等于转频说明是轴系不平衡或不对中问题不是轴承故障。第二种是随机冲击造成宽带抬升。现场环境中的敲击、手动敲击管道、异物的随机碰撞都会在包络谱中形成一片连续的隆起带形状像“土丘”而非尖锐谱线。区分办法是把同一段数据分成前后两半分别计算包络谱随机冲击造成的谱峰不会稳定出现在同一频率位置而轴承故障的特征峰两半数据都会保留。第三种是转速波动导致特征频率不是一条直线。变频驱动设备运行时转速在 29.5 到 30.5 Hz 之间摆动包络谱上原本尖锐的谱峰被拉宽成一个平台幅值大幅下降。这种情况下先做角域重采样把信号从时间域映射到角域再做包络谱分析比盲目提高频域分辨率有效得多。5. 把包络谱诊断流程封装成可复用函数现场验证之前先做这些5.1 一个简单可靠的emd_env_spec函数定长调用重复写十次以后封装成函数是必然选择。下面这个函数把“带通 → EMD → 选 IMF → 包络谱”四步压在一起适用绝大多数旋转机械。函数的输入输出都做了最简设计便于嵌入已经建好的数据批处理脚本里。function [f, env_spec, c_params] emd_env_spec(x, fs, opts) % EMD包络谱诊断工具函数 % 输入x 振动信号行向量 % fs 采样率 % opts.f_low, opts.f_high 带通 % opts.nfft 包络谱点数 % 输出f 频率轴env_spec 包络谱幅值c_params 选中的IMF参数 arguments x (1,:) double fs (1,1) double {mustBePositive} opts.f_low (1,1) double 2000 opts.f_high (1,1) double 6000 opts.nfft (1,1) double 8192 end % 1. 带通预滤波 [b, a] butter(4, [opts.f_low, opts.f_high] / (fs / 2), bandpass); x_fil filtfilt(b, a, x); % 2. EMD分解 imf emd(x_fil, MaxNumIMFs, 8); n_imf size(imf, 1); % 3. 用峭度和相关系数挑IMF K zeros(n_imf, 1); C zeros(n_imf, 1); for k 1:n_imf imf_k imf(k, :); m mean(imf_k); s std(imf_k); K(k) mean((imf_k - m).^4) / (s^4) - 3; C(k) abs(sum(imf_k .* x) / sqrt(sum(imf_k.^2) * sum(x.^2))); end cand find(C 0.3); [~, max_k_idx] max(K(cand)); imf_used cand(max_k_idx); % 4. 包络谱 env abs(hilbert(imf(imf_used, :))); env env - mean(env); win hanning(length(env)); env_win (env .* win) * opts.nfft / sum(win); Y fft(env_win, opts.nfft); f (0:opts.nfft-1) / opts.nfft * fs; env_spec 2 * abs(Y(1:opts.nfft/21)); f f(1:opts.nfft/21); c_params.K K(imf_used); c_params.C C(imf_used); c_params.imf_index imf_used; fprintf(选用第 %d 个IMF峭度 %.2f相关系数 %.3f\n, ... imf_used, K(imf_used), C(imf_used)); end函数默认带通范围是 2000 到 6000 Hz这是现场常见共振带的一个保守先验。真实应用中要把opts.f_low和opts.f_high按功率谱预估的共振带替换掉比默认值更重要。输出里c_params结构体保存了选中的 IMF 序号、峭度和相关系数方便批处理时把所有测点信息汇总到一张表中归档。5.2 现场报告里三个必做的验证动作第一个动作把理论特征频率与谱峰频率做比对设定 ±0.5% 的容差窗口。取 30 Hz 转频、BPFO 96.3 Hz 的情况窗口宽度就是约 0.48 Hz在这个范围内找到峰值才算命中。转速稳定时这个容差足够变频工况下把容差放宽到 ±2%并把包络谱主峰能量占比一并列进报告作为“疑似故障可信度”的量化指标。第二个动作至少截取两段不重叠的数据分别分析。数据一取第 0 到第 0.5 秒、数据二取第 0.5 到第 1 秒分别在两段上重跑一次emd_env_spec两个结果在特征频率处都出现稳定谱峰才在报告里写上“确认”。只跑一段数据就下结论是包络谱误诊断的第一大来源。第三个动作把包络谱和原始信号频谱叠在同一张图上。原始频谱的共振带位置用半透明色块标出包络谱的谱峰画出竖线。如果包络谱主峰频率对应的是原始频谱共振带内的某个边带中心诊断逻辑才是自洽的如果主峰落在共振带外先检查包络谱是不是被滤波器通带边缘的振铃污染了。画图时把包络谱横轴上界设在 10 倍特征频率以内不要直接默认到 fs/2这样谱线间距在图上更宽相邻故障谐波的分辨反而更清楚。本文还有配套的精品资源点击获取