简介面向机械故障诊断与信号处理学习者的EEMD集成经验模态分解MATLAB程序包围绕非线性、非平稳振动信号分解与轴承故障特征提取展开。压缩包共11个文件约761KB以9个.m脚本为主包含EEMD主流程及极值点、包络、IMF筛选、频谱与显著性分析等辅助函数另附1个.mat振动数据便于直接试跑1个.asv为编辑器备份文件。已有2034人学习。程序通过向原始信号加入小幅白噪声并多次执行EMD再取平均缓解模式混叠与噪声干扰输出若干IMF分量便于观察轴承异常振动对应的频率成分。使用者可替换自有数据完成分解、分量选择、包络与频谱可视化并据此理解滚珠磨损、内圈裂纹等故障特征的提取思路适合具备MATLAB基础、希望把EEMD用于设备状态监测与故障预测的研究人员和工程师参考。1. 从一段轴承振动信号说起EEMD 为什么值得单独拆一遍拿一段内圈剥落的轴承振动信号直接做包络解调频谱上常常是一堆糊在一起的峰故障特征频率被共振带压得看不见。原因不复杂传感器拾到的是多个激励源、多条传递路径叠出来的非平稳信号直接分析几乎必然失败。EEMD集成经验模态分解就是用来做这一步的往原始信号里成组加入高斯白噪声反复做经验模态分解再取平均把混叠的模态按尺度分开得到一组物理意义更清楚的 IMF 分量。这份 MATLAB 程序包做的正是轴承故障诊断里「信号预处理 → 模态分离 → 包络谱找故障特征频率」这条链路。适合两类人手上有一堆实测振动数据、想跑通一次完整诊断流程的以及只想拿一份能改的 EEMD 参考实现、塞进自己项目里的。2. 包内结构与信号链路rar 里那几个 .m 文件各自管什么.rar解开之后不会是一个巨型脚本正常情况是按「入口脚本 算法核心 工具函数 数据」四层摊开的。先看骨架再看代码能省掉大量来回跳转的时间。这个项目的价值在落地它把从原始振动数据到一张能看懂的包络谱图中间的每一步都摆出来了哪一步改了会出什么后果代码里基本都能追到。2.1 典型目录构成与文件职责不同人打包习惯不一样但这类程序包的文件职责划分大同小异常见做法是下面这套。你解压之后先对着表把文件名对一遍名字对不上没关系看每个文件开头几行在干什么归类是一样的。文件 / 目录承担的职责读代码时该盯的地方eemd_main.m入口脚本读数据、设参数、调分解、画图默认参数都写在这里改参数先改它eemd_core.m集合循环加噪、调用单次 EMD、累加平均NE、Nstd的实际用法和单位emd_single.m单次 EMD 筛分返回 IMF 矩阵与残余分量停止准则、最大层数、残余判据extrema_envelope.m找局部极值、三次样条拟合上下包络端点怎么处理有没有延拓imf_select.m按相关系数 / 峭度挑出用于解调的分量指标计算口径归一化没归一化env_spectrum.mHilbert 包络、FFT、标注故障特征频率横轴是归一化频率还是 Hzdata/.mat或.txt振动数据采样率记在哪个变量或注释里这份拆法的好处是每一层都能单独拿出来测。我一般先把emd_single.m单独跑一个已知的仿真调幅信号看它拆出来的分量对不对再去动eemd_core.m的集合参数。如果直接从入口脚本一路跑到黑中间哪一层有问题你根本分不清只能靠猜这就成了玄学调试。2.2 EEMD 的算法链路噪声集合、单次 EMD 与 IMF 平均EEMD 的核心思想一句话能说完单次 EMD 在极值点分布不均匀的地方会出现模态混叠那就人为地把极值点补均匀。往信号里加一组有限幅度的高斯白噪声噪声在整个时频平面上均匀分布等于给每个尺度都补上了参照极值点分解出来的 IMF 就不会在不同尺度之间跳来跳去。单次加噪的结果带有随机性所以要做NE次最后对应层取平均噪声在平均里互相抵消剩下的是稳定的分量。链路是四步读入振动信号x估计其标准差sigma第e轮生成Nstd * sigma * randn(size(x))的白噪声叠到x上得到x_e对x_e跑一次完整的 EMD得到该轮的 IMF 集合NE轮结束后按层号对齐求平均得到最终 IMF 和残余分量。这里面有一个常被忽略的定量关系平均之后残留的噪声误差约等于Nstd / sqrt(NE)。想要残留误差减半NE得翻四倍。这就解释了为什么大家默认NE取 100 到 200 而不是 20——取 20 的时候那个残留噪声足以在包络谱上压出假峰。理解了这个平方根关系你调参就不会盲目往上堆集合次数了。2.3 从分解结果到故障特征频率包络谱在哪一步接上轴承局部损伤的典型表现是每转到一个损伤点就产生一次冲击冲击激起轴承座和传感器的高频共振。共振把故障信息「搬运」到了高频段在原始频谱上表现为一个共振带故障特征频率本身反而淹没在里面。标准的处理顺序是先用 EEMD 把高频共振那一层 IMF 单独分出来再对这条 IMF 做 Hilbert 包络相当于把高频载波剥掉只留下冲击的包络最后对包络做 FFT谱峰就落在故障特征频率及其谐波上。所以env_spectrum.m里那几行是整个诊断流程的判决环节。横轴的换算必须用真实采样率f (0:N-1) * fs / N。这一步错了后面所有峰值判读全部作废而且错得很隐蔽——图上明明有漂亮的谱峰就是和理论频率差着一个固定倍数。用峭度或相关系数选 IMF 也在这条链路上峭度对冲击敏感相关系数保证选出来的分量和原信号还有足够的关联两个指标一起看比只看峭度稳。3. 核心代码拆解主循环、加噪幅度与 IMF 筛选算法核心其实就两个函数。eemd_core.m管集合emd_single.m管单次筛分。前者短但坑多后者长但逻辑是标准的。下面按我习惯的写法拆跟原始文件的变量名可能对不上结构是一致的。3.1 EEMD 主循环骨架function [imf, res] eemd_core(x, NE, Nstd, maxIMF, stopCri) % x : 振动信号行/列均可 % NE : 集合次数常用 100~200 % Nstd : 噪声幅度相对信号标准差的倍数常用 0.1~0.3 % maxIMF : 最大分解层数上限 % stopCri: 单次筛分的标准差停止阈值 x x(:); % 统一成行向量后面 size 才不会打架 N numel(x); sigma std(x); maxIMF min(maxIMF, floor(log2(N)) - 1); % 层数上限防止过分解出虚假分量 imfSum zeros(maxIMF, N); cnt zeros(maxIMF, 1); % 每层实际被累加了几次 for e 1:NE w Nstd * sigma * randn(1, N); % 每轮必须重新生成不能复用 xe x w; [imfE, ~] emd_single(xe, maxIMF, stopCri); m size(imfE, 1); % 本轮实际分解出的层数 m min(m, maxIMF); imfSum(1:m, :) imfSum(1:m, :) imfE(1:m, :); cnt(1:m) cnt(1:m) 1; end cnt(cnt 0) 1; imf imfSum ./ cnt; % 按层分别平均而不是统一除以 NE res x - sum(imf, 1); % 残余趋势项 end先说cnt这个数组。它不是可有可无的保险而是这个实现里最关键的一行设计。不同噪声下 EMD 分解出来的层数经常不一样某一轮只出了 6 层另一轮出了 8 层如果你图省事写成imfSum imfSum imfE再统一/ NE第 7、8 层就会被少加几次平均后的幅值整体偏低包络谱上表现为某一层能量莫名很弱。分开计数之后每一层只除以它自己实际参与的次数这才是正确的平均。再说randn(1, N)的位置。它必须写在循环体内每轮重新采样。我见过把噪声生成提到循环外面、NE次都用同一段噪声的写法那样平均之后噪声根本抵消不掉EEMD 退化成「一次加噪的 EMD」模态混叠一点没改善白跑一百轮。参数上还有两个可改的地方。maxIMF用floor(log2(N)) - 1封顶是个经验做法N取 10000 时上限大概 12 层实际用不到这么多轴承信号一般取 5 到 8 层就够。stopCri传进来交给单次筛分用默认给 0.2给太小筛分次数暴涨但分量未必更好。如果要跑批量数据把for换成parfor就行注意imfSum和cnt是累加变量parfor不允许这样归约改法是用parfor只做每轮的 EMD把结果塞进 cell循环外再统一平均。这个改造在最后一章展开。3.2 单次 EMD 的筛分过程与三次样条包络function [imf, r] emd_single(x, maxIMF, stopCri) r x; imf []; for k 1:maxIMF h r; for s 1:20 % 单层筛分迭代上限 [idxMax, idxMin] local_extrema(h); if numel(idxMax) 3 || numel(idxMin) 3 break; % 极值点太少没法拟合包络 end t 1:numel(h); emax spline(idxMax, h(idxMax), t); % 上包络 emin spline(idxMin, h(idxMin), t); % 下包络 m (emax emin) / 2; hPrev h; h h - m; sd sum((hPrev - h).^2) / (sum(hPrev.^2) eps); % 标准差准则 if sd stopCri break; end end if isempty(h) || all(~isfinite(h)) break; end imf [imf; h]; r r - h; [iMax, iMin] local_extrema(r); if numel(iMax) numel(iMin) 4 % 残余已是单调趋势停止分解 break; end end end这段是 EMD 的标准筛分。逻辑上就是把信号减去上下包络的均值反复减直到剩下的部分满足「极值点数目与过零点数目相差不超过 1、上下包络均值为零」这两个 IMF 条件。工程上没人真去严格判定这两条都用一个类似柯西收敛准则的sd阈值代替sd越小筛得越狠。几个参数的实际含义内层20是单层筛分次数硬上限防止某层筛不下来死循环stopCri一般取 0.2 到 0.3取 0.05 会让筛分次数翻好几倍分解出来的 IMF 更「纯」但容易出现幅值调制失真。外层那个「残余极值点少于 4 个就停」的判据很重要没有它程序会把一条单调趋势硬拆成两三层虚假 IMF这些假分量能量很小但会在包络谱低频谱段制造假峰。spline是这段代码里最脆弱的一环。三次样条在数据两端没有约束会向外发散。信号首尾几个采样点如果刚好不是极值点包络线在两端就会飞出去减完之后 IMF 的两端出现大幅摆动这就是常说的端点效应。轴承信号做循环平稳分析时这个问题尤其明显包络谱上会出现一片低频抬升。处理办法下一章讲。3.3 参数怎么定NE、Nstd、maxIMF 与阈值的取值和后果参数常用取值作用调大的后果调小的后果NE集合次数100 ~ 200平均次数压残留噪声耗时线性增长收益按平方根递减残留噪声大IMF 幅值抖动包络谱假峰Nstd噪声幅度0.1 ~ 0.3 倍信号标准差补充各尺度极值点噪声淹没弱冲击分量产生伪 IMF极值点补不够模态混叠改善不明显maxIMF最大层数5 ~ 8限定分解深度末尾几层是虚假分量与端点效应高频共振层没分出来信息丢失stopCri筛分阈值0.2 ~ 0.3单层筛分的收敛判据筛分不充分IMF 混叠筛分过度幅值失真、耗时暴涨Nstd和NE要一起看。残留噪声误差约等于Nstd / sqrt(NE)如果你把Nstd从 0.2 提到 0.4 想改善混叠等于把残留噪声翻倍这时候NE得从 100 提到 400 才能回到原来的噪声水平。我一般的做法是Nstd固定在 0.2先跑NE 50看趋势对不对确认没问题再提到 100 或 200 出正式结果。这个两步走的习惯救过我不少次——参数设错的时候50 轮跑 20 秒就能发现200 轮跑三分钟才发现就有点亏。还有一点Nstd是相对倍数代码里乘了sigma。如果你的数据预先做过归一化或者量纲是 mV 而不是 m/s²忘记这一点直接传绝对幅值加进去的噪声要么可以忽略要么直接把信号埋掉。4. 跑通流程从加载振动数据到出包络谱4.1 环境准备与路径约定MATLAB 版本上没太多讲究R2016b 之后基本都能跑。核心分解部分只用到spline、randn、矩阵运算这些基础函数不需要额外工具箱包络那一步如果用了hilbert就需要 Signal Processing Toolbox。手上没有这个工具箱也没关系解析信号可以用 FFT 自己构造把正频率部分乘 2、负频率置零、再 IFFT效果一样。目录上建议保持原包的结构把data/和几个.m文件放在同一层入口脚本里用相对路径读数据。MATLAB 的当前文件夹和搜索路径是两回事双击脚本运行时当前文件夹会跟着变如果你的数据路径写成了绝对路径换台机器就报「文件不存在」。我一般会在入口脚本开头加一行cd(fileparts(mfilename(fullpath)))把工作目录强制切到脚本所在位置这个习惯能省掉一半的路径类问题。数据格式上包内示例大概率是.mat里面存一个行向量或列向量。读进来第一件事是确认三样东西长度、采样率、量纲。采样率有时写在.mat的另一个变量里有时只写在脚本注释里找不到就问数据来源别猜。4.2 主脚本分步执行与结果判读%% eemd_main.m 分步执行版 clear; clc; close all; cd(fileparts(mfilename(fullpath))); fs 12000; % 采样频率 Hz必须与采集设置一致 S load(data/bearing.mat); x S.x(:); % 统一成行向量 N numel(x); t (0:N-1) / fs; NE 100; Nstd 0.2; maxIMF 8; stopCri 0.2; tic; imf eemd_core(x, NE, Nstd, maxIMF, stopCri); fprintf(EEMD 用时 %.1f s得到 %d 层 IMF\n, toc, size(imf,1)); k imf_select(imf, x); % 综合相关系数峭度挑一层 fprintf(选中第 %d 层 IMF 做包络解调\n, k); env_spectrum(imf(k,:), fs, [79.7, 120.3, 9.96, 58.9]); % 后四个是理论特征频率跑完之后重点看三处输出。第一是耗时NE100、数据长度 1 万点时通常在几十秒到几分钟量级如果跑了半小时还没结束八成是maxIMF或筛分次数设太大。第二是 IMF 层数正常在 5 到 8 层之间出来 12 层说明上限没封住或者信号里有大量噪声。第三是imf_select选中的层号轴承信号通常落在第 2 到第 4 层也就是高频共振带对应的那几层如果选中了最后一两层多半是峭度指标被端点效应污染的虚假分量骗了。包络谱图出来之后判读的顺序是先看低频段有没有明显的谱峰族再看这些峰的间距是否等于某个理论特征频率最后看有没有以该频率为基频的谐波序列。单看一个峰不算数真实故障一定伴随谐波而且内圈故障还会有转频调制形成的边频带。这些判据不看图是说不清的跑一遍对着图看最直观。4.3 换成自己的数据采样率、转速与故障特征频率换算换成自己的数据唯一需要动脑子的是理论特征频率的计算。轴承的四个特征频率由几何参数和转速决定我一般先按下面这套公式算一遍把结果作为标注传给env_spectrum。特征频率计算公式说明外圈 BPFO(n/2) · fr · (1 - (d/D)·cosα)冲击落在固定位置最容易被检出内圈 BPFI(n/2) · fr · (1 (d/D)·cosα)随轴旋转常带转频调制边频滚动体 BSF(D/(2d)) · fr · (1 - ((d/D)·cosα)²)幅值通常较弱保持架 FTF(fr/2) · (1 - (d/D)·cosα)频率低易与转频谐波混淆式中n为滚动体个数fr为轴转频Hzd为滚动体直径D为节圆直径α为接触角。举个例子n8、d7.94 mm、D39.04 mm、接触角为 0、转速 1500 r/min 即fr25 Hz算出来 BPFO 约 79.7 Hz、BPFI 约 120.3 Hz、BSF 约 58.9 Hz、FTF 约 9.96 Hz。换成自己的数据时要改的只有三处fs改成真实采样率fr按实际转速换算几何参数按轴承型号查。前两个错了包络谱直接对不上第三个错了你会去找一个根本不存在的频率。另外注意实际测到的特征频率通常比理论值偏 1% 到 2%这是滚动体打滑造成的判读时留一点容差别死磕小数点。5. 常见问题排查IMF 层数错位、端点飞起与包络谱假峰这一节列的都是实际跑的时候真会遇到的。每条按「现象 → 原因 → 解决」写对着症状找。5.1 运行报错与维度类问题现象Matrix dimensions must agree或者循环里某一行报错位置每次都不同。原因行向量和列向量混用。入口脚本读进来的x是列向量函数内部按行向量处理x w直接维度不匹配另一种情况是imfE和imfSum的层数不同却没做截断。解决在eemd_core入口处统一x x(:)在累加处用m min(size(imfE,1), maxIMF)显式截断别指望 MATLAB 自动对齐。现象结果是全NaN或者一片Inf。原因Nstd传的是绝对幅值而不是倍数噪声把信号淹了之后sd分母除零或者采样率参数写成 0后面fs/N除以零。解决eemd_core里加一句断言assert(Nstd 1 Nstd 0)把量纲问题挡在入口。这个断言我建议一直留着。5.2 结果类问题模态混叠没改善、IMF 数量对不上现象跑完 EEMD 发现和前一次跑的结果差别很大IMF 波形抖动。原因NE太小残留噪声Nstd/sqrt(NE)太大平均没压住随机性或者噪声生成写在了循环外。解决先确认randn在循环体内再把NE提到 100 以上。判断标准很直观把NE从 100 提到 200如果 IMF 波形明显稳定下来了说明之前就是集合次数不够。现象平均之后个别 IMF 层的幅值明显偏低或者出现整层接近零。原因不同轮次分解出来的层数不同累加时没按层分别计数用了统一除数NE。解决就是 3.1 里那个cnt数组按层分别除。这是这份程序包里最容易写错、也最难从图上发现的一处。现象IMF 两端大幅飞起图上一看就是发散的。原因三次样条没有端点约束首尾极值点外推发散。解决在local_extrema之后、spline之前做一次镜像延拓把信号首尾各镜像一段再拟合包络拟合完截掉延拓部分省事的做法是直接裁掉 IMF 两端各 5% 的采样点再做后面的分析代价是损失一点数据。5.3 诊断类问题包络谱峰值对不上故障频率现象包络谱上峰很漂亮但位置和理论特征频率差了将近一倍或者差着一个固定比例。原因横轴换算用的是归一化频率而不是真实频率f (0:N-1)/N少了fs这个乘子。解决统一用f (0:N-1) * fs / N并且只画前一半到fs/2后半段是镜像。这个错误看起来低级但改参数时把fs覆盖掉的情况真的会发生。现象什么峰都看不出来包络谱一片平。原因选错了 IMF把低频趋势层拿去做解调了或者信号本身信噪比太低故障冲击根本没被激发出来。解决先人工把每一层 IMF 都画出来肉眼看哪一层有明显的周期性冲击再对照峭度值排序。如果所有层都很平问题在数据采集端不在算法换算法也救不回来。6. 二次开发把它改成能批量跑的诊断流水线单文件跑通只是起点。真实项目里你手上是几十上百条测点数据一条条改参数、存图、记结论做上十条就想砸键盘。我通常按下表这几个方向改投入产出比从高到低排改造点做法验证方式并行化每轮 EMD 塞进parfor结果存 cell循环外统一平均与串行结果逐元素比对误差应小于 1e-10自动选层峭度与相关系数归一化后加权打分取最高分那层人工挑层结果与自动结果比对看命中率批量处理外层遍历文件结果存 struct 数组含参数快照抽 3 条数据手工复算核对特征频率与判定结论算法替换把固定噪声换成 CEEMDAN 的自适应噪声同一数据上比较重构误差与 IMF 正交性指标参数快照每次运行把NE/Nstd/maxIMF与数据哈希写进结果隔一周重跑能完全复现上一次的图并行化那段的关键改动在于parfor不允许对imfSum做归约累加得改成下面这个形状imfCell cell(NE, 1); parfor e 1:NE w Nstd * sigma * randn(1, N); [imfE, ~] emd_single(x w, maxIMF, stopCri); imfCell{e} imfE; end % 循环外对齐层数再平均以最小层数为准避免维度错位 mMin min(cellfun((c) size(c,1), imfCell)); imf zeros(mMin, N); for e 1:NE imf imf imfCell{e}(1:mMin, :); end imf imf / NE;mMin那一行是这段代码能不能用的分水岭。并行之后各轮的层数差异比串行更明显直接按第一轮的长度切片后面层数少的那些轮次会下标越界。取所有轮次的最小层数代价是丢掉末尾一两层——那些层基本都是趋势项和虚假分量丢掉不影响诊断结论。再往上一层是算法替换。EEMD 的固定噪声幅度Nstd是个需要人工试的参数换成 CEEMDAN 之后噪声自适应重构误差更小代价是计算量翻倍、代码复杂度上一个台阶。我的建议是先拿 EEMD 把整套流水线跑顺、把特征频率判据和阈值标定好再考虑换算法反过来做算法一换所有阈值都得重标前面的标定工作全废。最后说一个血泪经验。这类信号处理代码最要命的不是写错而是「跑出来一张好看的图然后你信了」。我曾经因为漏了按层计数的cnt某一层幅值整体偏低包络谱上恰好少了一个谐波结论直接从「内圈故障」变成「疑似正常」。从那以后我每次改完核心循环都会先用一段构造好的仿真信号——两个已知频率的调幅分量叠加——跑一遍看分解出来的 IMF 频率和幅值是不是和设定值对得上对不上就不往下走。这道自检花不了五分钟但它是我唯一信的后悔药。希望帮到你。本文还有配套的精品资源点击获取