简介EEMD集成经验模态分解MATLAB程序面向从事信号处理与机械故障诊断的研究人员和工程师尤其适合需要分析非线性、非平稳振动信号的轴承故障诊断场景。资源包共11个文件以9个m脚本为核心涵盖EEMD主流程、极值点提取、包络计算、IMF分量筛选与显著性检验等模块另含1个mat振动数据样本和1个asv备份文件整体约761KB结构紧凑便于直接运行与二次开发。已有2034人学习下载说明其在相关领域具有较高的参考价值。读者可借助该程序加载自有振动信号完成EEMD分解、IMF分量选取与频谱可视化进而提取轴承故障特征频率判断滚珠磨损、内圈裂纹等异常模式是理解并实践EEMD方法的实用工具。1. 从一段抖动信号说起EEMD 集成经验模态分解 MATLAB 程序能帮你做什么如果你手上有轴承振动信号想做故障诊断大概率绕不开一个经典问题原始信号里混着噪声、调制成分和冲击成分直接做 FFT 只能看到一堆频率峰根本分不清哪个是故障特征。经验模态分解EMD本来是个好工具能把信号自适应地拆成若干本征模态函数IMF但 EMD 有个玄学毛病——模态混叠同一个 IMF 里混进了时间尺度差异很大的成分分解结果不稳定换个采样点结果就变了。EEMDEnsemble Empirical Mode Decomposition集成经验模态分解就是来治这个病的。它的思路很直接往原始信号里反复加入不同幅值的高斯白噪声做多次 EMD再把多次分解得到的 IMF 按对应阶数取平均。白噪声把原本可能混叠的成分“打散”到不同尺度上平均之后噪声互相抵消剩下的就是更干净的 IMF。这套 MATLAB 程序就是干这件事的适合做旋转机械故障诊断、信号处理课程设计、或者需要把 EEMD 当预处理步骤的工程开发者。拿到手能直接跑也能拆开改参数、换数据、接自己的诊断流程。2. 程序结构与核心模块拆解从 main 脚本到 EEMD 主循环2.1 目录结构与文件职责拿到这个 .rar 解压后常见的目录组织是这样的不同版本可能略有差异但核心文件跑不掉EEMD_Program/ ├── main.m % 入口脚本加载信号、调 EEMD、画图 ├── eemd.m % EEMD 核心函数 ├── emd.m % 单次 EMD 分解或调用工具箱 ├── get_imf.m % IMF 提取与筛选辅助 ├── plot_result.m % 分解结果可视化 ├── data/ │ └── bearing_signal.mat % 示例轴承振动数据 └── utils/ ├── extrema.m % 极值点检测 └── spline_envelope.m % 包络线拟合main.m是入口负责读数据、设参数、调eemd.m、出图。eemd.m是核心里面是集成循环。emd.m是单次分解负责筛分过程sifting。utils里是底层工具函数极值点检测和包络拟合是 EMD 的命根子这两个函数写得好不好直接决定分解质量。提示如果你机器上已经装了 MATLAB 的 EMD 工具箱emd.m可能只是个 wrapper如果没有程序里通常自带一个简化版实现精度够用但速度一般。2.2 EEMD 主循环加噪、分解、平均这是整个程序最核心的一段逻辑我把它简化成可读的伪代码形式方便你看清每一步在干什么function [imf, res] eemd(signal, Nstd, NE, maxIMF) % signal : 输入信号列向量 % Nstd : 添加白噪声的标准差比例相对信号标准差 % NE : 集成次数ensemble number % maxIMF : 最大 IMF 阶数 N length(signal); allIMF zeros(N, maxIMF, NE); % 存每次分解的 IMF for k 1:NE % 1. 生成白噪声并加到原信号上 noise Nstd * std(signal) * randn(N, 1); noisySignal signal noise; % 2. 对加噪信号做一次标准 EMD imf_k emd(noisySignal, maxIMF); % 3. 对齐阶数后存入三维数组 nIMF size(imf_k, 2); allIMF(:, 1:nIMF, k) imf_k; end % 4. 沿集成维度取平均得到最终 IMF imf mean(allIMF, 3); res signal - sum(imf, 2); % 残余分量 end逻辑说明第 1 步的Nstd控制加噪强度典型值 0.2 到 0.5太小起不到打散作用太大残留噪声多。第 2 步调emd做筛分每次分解出的 IMF 数量可能不一样所以第 3 步要做对齐不够的补零。第 4 步的mean(allIMF, 3)是 EEMD 的灵魂——沿集成维度平均噪声项因为零均值被压下去真实成分被保留。参数说明NE一般取 50 到 200取小了平均效果不够取大了计算时间线性增长。maxIMF通常设 8 到 12轴承信号一般前几阶就包含故障特征。Nstd的经验值是 0.2这是 EEMD 原始文献推荐的默认值。2.3 单次 EMD 的筛分过程emd.m里是经典的筛分循环核心就三步找极值、拟合包络、减均值。function imf emd(signal, maxIMF) imf []; r signal; % 残余信号初始为原信号 for i 1:maxIMF h r; for sift 1:10 % 筛分迭代通常 10 次以内收敛 % 找极大值和极小值点 [maxPos, maxVal] findpeaks(h); [minPos, minVal] findpeaks(-h); minVal -minVal; % 三次样条拟合上下包络 upperEnv spline(maxPos, maxVal, 1:length(h)); lowerEnv spline(minPos, minVal, 1:length(h)); % 包络均值 meanEnv (upperEnv lowerEnv) / 2; % 减去均值得到新的 h h_new h - meanEnv; % 收敛判断连续两次 h 的差异足够小 if norm(h_new - h) / norm(h) 0.01 break; end h h_new; end imf [imf, h]; % 保存这一阶 IMF r r - h; % 更新残余 % 残余如果单调或极值点太少停止分解 if length(findpeaks(r)) 2 break; end end end逻辑说明筛分的目的是把信号里最快变化的成分最高频先剥离出来作为第一阶 IMF。每次迭代用上下包络的均值去逼近信号的局部趋势减掉之后剩下的就是振荡成分。收敛条件用相对误差判断0.01 是个常用阈值。参数说明maxIMF控制最多分解几阶实际阶数由残余信号的极值点数量决定。筛分迭代上限设 10 是防止死循环正常信号 5 到 8 次就收敛了。spline用三次样条边界处理是 EMD 的经典难点程序里通常用简单端点延拓效果一般但够用。2.4 结果可视化与 IMF 筛选分解完之后plot_result.m会把原始信号和各阶 IMF 画成子图方便你肉眼判断哪一阶包含故障特征。轴承故障诊断里通常看前 3 到 5 阶 IMF计算它们的包络谱或者峭度值峭度大的那阶往往对应冲击成分。% 计算各阶 IMF 的峭度辅助筛选 for i 1:size(imf, 2) kurtVal(i) kurtosis(imf(:, i)); end [~, idx] max(kurtVal); % 峭度最大的 IMF 阶数 fprintf(峭度最大的 IMF 阶数: %d\n, idx);这段代码不是必须的但做故障诊断时很实用。峭度对冲击成分敏感轴承外圈或内圈故障会在特定 IMF 上表现为高峭度。你可以把idx对应的 IMF 拿去做包络谱找故障特征频率。3. 跑起来环境准备、参数设置与第一次运行3.1 MATLAB 版本与依赖检查这个程序对 MATLAB 版本要求不高R2016a 以上基本都能跑。用到的函数都是基础信号处理函数findpeaks、spline、randn、kurtosis不需要额外工具箱。如果你用的是 R2020b 之后的版本findpeaks的行为有细微变化极值点检测的阈值参数可能需要微调。检查依赖的方法很简单在命令行敲which findpeaks which spline which kurtosis三个都有路径返回就说明环境没问题。如果findpeaks报未定义说明你装的是精简版 MATLAB需要补装 Signal Processing Toolbox。3.2 加载数据与参数初始化main.m开头通常是这样的%% 1. 加载信号 load(data/bearing_signal.mat); % 变量名通常是 signal 或 x fs 12000; % 采样频率按实际数据改 t (0:length(signal)-1) / fs; %% 2. EEMD 参数 Nstd 0.2; % 噪声标准差比例 NE 100; % 集成次数 maxIMF 8; % 最大 IMF 阶数 %% 3. 执行 EEMD [imf, res] eemd(signal, Nstd, NE, maxIMF); %% 4. 绘图 plot_result(signal, imf, res, fs);参数怎么改fs必须和你数据的实际采样率一致否则包络谱的频率轴全错。NE从 100 起步如果结果不稳定就加到 200。Nstd先试 0.2分解出来的 IMF 如果模态混叠还严重加到 0.3 或 0.4。maxIMF设 8 够用轴承信号一般 5 阶以内就有明显特征。3.3 运行时间与性能预期EEMD 的计算量是NE乘以单次 EMD 的时间。单次 EMD 对 10000 点信号大概 0.1 到 0.3 秒NE100就是 10 到 30 秒。如果信号更长或者NE更大跑几分钟很正常。我一般会先把NE设成 20 快速跑一遍看结果趋势确认参数合理后再加到 100 出正式结果。注意MATLAB 的spline在数据点很多时是瓶颈如果嫌慢可以把emd.m里的样条换成interp1的pchip方法速度快一倍但包络平滑度略降对最终结果影响不大。3.4 结果怎么看跑完之后你会看到一张图最上面是原始信号下面依次是 IMF1 到 IMFn最后是残余分量。判断分解是否合理的标准IMF1 频率最高往后频率逐渐降低残余分量应该是单调趋势或者极缓变化各阶 IMF 之间没有明显的频率重叠。如果 IMF1 里能看到明显的低频调制说明模态混叠没压住把Nstd调大或者NE调大再跑。如果各阶 IMF 看起来都差不多可能是maxIMF设太大了后面几阶已经是噪声减到 5 或 6 再试。4. 避坑与排查EEMD 调参和运行中的 5 个常见问题4.1 分解结果每次运行都不一样现象同样的数据跑两次eemd得到的 IMF 数值有差异。原因EEMD 内部用randn生成白噪声MATLAB 默认的随机种子每次启动都不同导致每次加噪序列不一样平均后的结果自然有微小差异。解决在eemd.m开头加一行rng(42);固定随机种子。42 是随便选的你也可以用其他整数。固定之后每次运行结果完全一致方便复现和对比。4.2 IMF 阶数太多后面几阶全是噪声现象maxIMF8但 IMF6 到 IMF8 看起来就是随机抖动没有物理意义。原因EEMD 会把噪声也分解成 IMF集成平均只能压掉一部分残余噪声会落到高阶 IMF 上。另外轴承信号本身的有效成分可能只有 3 到 4 阶。解决不要盲目设大maxIMF。先跑一次看峭度分布峭度低于 3 的高阶 IMF 直接丢弃。或者在eemd.m里加个判断当残余信号的极值点少于 2 个就提前退出不强行凑满maxIMF。4.3 包络谱里找不到故障特征频率现象分解看起来正常但拿 IMF 做包络谱故障特征频率比如外圈故障的 BPFO对不上。原因可能是 IMF 选错了故障特征不在你选的那一阶也可能是采样率fs设错了频率轴整体偏移还有一种可能是包络谱的解析信号构造有问题。解决先确认fs和数据采集时一致。然后把前 5 阶 IMF 的包络谱都画出来逐个对照理论故障频率。轴承故障频率有公式可算BPFO (n/2) * fr * (1 - (d/D)cosθ)其中 n 是滚珠数fr 是转频d 是滚珠直径D 是节径θ 是接触角。算出来之后在包络谱上找对应峰。4.4 程序跑一半报错“矩阵维度不一致”现象eemd.m里allIMF(:, 1:nIMF, k) imf_k;这行报错。原因不同次 EMD 分解出的 IMF 数量不一样某次可能只分解出 5 阶另一次分解出 7 阶nIMF超过maxIMF或者对齐逻辑有 bug。解决在赋值前加判断nIMF min(nIMF, maxIMF);确保不超过预分配维度。另外emd.m返回的imf如果是行向量而不是列向量也会导致维度不匹配统一用imf imf(:);转成列向量。4.5 大数据量下内存爆掉现象信号长度 100000 点NE200maxIMF10跑的时候 MATLAB 卡死或者报 Out of Memory。原因allIMF三维数组的大小是N × maxIMF × NE100000 × 10 × 200 2 亿个 double约 1.6 GB加上中间变量很容易超内存。解决两个办法。一是减小NE50 到 100 通常够用。二是分段处理把长信号切成 10000 点的段每段单独做 EEMD最后拼接。分段处理要注意段与段之间的边界效应重叠 10% 再取中间部分能缓解。5. 二次开发从能跑到好用几个改造方向5.1 把 EEMD 封装成函数句柄接入自己的诊断流程原始程序是个脚本跑完出图就结束了。实际项目里你多半想把 EEMD 当成一个预处理步骤后面接特征提取和分类器。改造方法是把main.m里的核心逻辑抽成一个函数function [featureVec] extract_eemd_feature(signal, fs) % 输入原始信号、采样率 % 输出特征向量各阶 IMF 的峭度 包络谱峰值 Nstd 0.2; NE 100; maxIMF 6; [imf, ~] eemd(signal, Nstd, NE, maxIMF); featureVec []; for i 1:size(imf, 2) % 峭度特征 k kurtosis(imf(:, i)); % 包络谱峰值特征 env abs(hilbert(imf(:, i))); envSpec abs(fft(env - mean(env))); [~, peakIdx] max(envSpec(1:floor(end/2))); peakFreq (peakIdx - 1) * fs / length(env); featureVec [featureVec, k, peakFreq]; end end这样你就能在循环里批量处理多个样本把featureVec喂给 SVM 或神经网络做分类。参数maxIMF设 6 是因为特征太多反而干扰分类器前 6 阶的峭度和峰值频率已经能区分大部分轴承故障类型。5.2 用并行计算加速集成循环NE100的时候串行跑要几十秒如果你有 Parallel Computing Toolbox把eemd.m里的for k 1:NE改成parfor k 1:NE速度能提升 4 到 8 倍取决于核数。注意parfor里不能有依赖上一次迭代的变量allIMF的切片赋值是合法的但rng的种子设置要小心每个 worker 的随机流是独立的固定种子后结果仍然可复现。5.3 替换 EMD 内核从样条包络到更快的实现如果你嫌spline慢可以把emd.m里的包络拟合换成envelope函数R2015b 之后自带它内部用希尔伯特变换速度快很多但包络定义和极值点拟合不完全等价分解结果会有差异。我的经验是做故障诊断时包络的精确形状不如极值点位置重要用envelope替代后故障特征频率仍然能对上但 IMF 的数值会变如果你有后续的阈值判断逻辑需要重新标定。5.4 验证改造是否成功一个简单的对照实验改完之后怎么知道没改坏拿一段已知故障类型的轴承信号分别用原始程序和改造后的版本跑对比三件事IMF 阶数是否一致、峭度最大的阶数是否一致、包络谱峰值频率是否一致。三个都对上说明改造没引入偏差。如果只有数值有微小差异但特征一致可以接受如果特征频率对不上说明改造动了核心逻辑得回去查。我自己的习惯是每次改完eemd.m或emd.m先跑一遍main.m确认出图和改之前一样再跑对照实验。这个习惯帮我省过好几次“以为改的是可视化结果动了分解逻辑”的后悔药。EEMD 这种多步迭代的算法中间任何一步的改动都可能被后续放大验证必须走完整流程。希望帮到你。本文还有配套的精品资源点击获取
