简介围绕EEMD集成经验模态分解的MATLAB程序包面向从事轴承故障诊断、机械振动信号分析的研究人员与工程技术人员也可供学习非线性、非平稳信号处理的学生参考。它针对传统EMD易出现模式混叠、噪声干扰等问题通过白噪声辅助与多次分解平均提升IMF提取稳定性便于从振动信号中分离异常频率成分并开展故障特征分析。压缩包共11个文件以9个m脚本为主体另有1个mat数据文件和1个asv备份整体约761KB体积轻便适合直接放入MATLAB环境调试与二次修改。已有2034人学习说明其在相关学习中具有一定参考价值。使用者可据此加载自己的振动数据完成EEMD分解、IMF分量筛选与结果可视化进而理解噪声辅助分解流程、参数设置和轴承故障特征提取思路为设备状态监测与早期故障识别提供可复用的代码基础。1. 从一段被噪声淹没的轴承振动信号说起轴承早期故障的冲击成分非常弱一段 12 kHz 采样、上万点的加速度信号直接丢进 FFT谱线上全是转频谐波和背景噪声故障特征频率根本冒不出来。EEMD集成经验模态分解处理的就是这种非平稳、低信噪比的场景往原始信号里多次叠加不同幅值的白噪声再做 EMD利用噪声在频谱上均匀分布的特性把模态混叠压下去最后对多次分解结果取集成平均得到一组按时间尺度排列的 IMF 分量。挑出承载故障冲击的高频 IMF 做 Hilbert 包络 demodulation包络谱上才会出现清晰的 BPFO、BPFI 及其倍频。这份 MATLAB 程序把这条链路从数据读入、EEMD 分解、IMF 筛选到包络谱出图串成了可运行的原型适合手里有台架或现场振动数据、想做轴承故障特征提取验证的开发者也适合刚接触经验模态分解、需要一份能跑起来再动手改的代码作参照。2. 程序结构与技术选型一个 rar 包里装的是什么2.1 目录骨架与函数分工rаr 解压之后正常的工程目录不会把所有东西塞进一个脚本文件里而是按“分解核心 / 辅助函数 / 数据 / 出图”拆开。我拆这类包最先看的是eemd.m和emd.m的调用关系因为它决定了整套代码是自研实现还是套用现成的 EMD 内核。常见的目录形态大概是这样eemd_bearing/ ├── main_eemd.m % 主入口读数据、调分解、画图 ├── eemd.m % 集成经验模态分解负责加噪与平均 ├── emd.m % 单次 EMD 筛分过程 ├── extrema.m % 极值点定位局部极大/极小 ├── envelope.m % 三次样条上下包络插值 ├── utils/ │ ├── read_signal.m % 读取 txt / csv / mat 振动数据 │ ├── select_imf.m % IMF 筛选相关系数或峭度准则 │ └── envelope_spectrum.m % Hilbert 包络 频谱 ├── data/ │ └── bearing_inner.mat % 示例数据内含信号与采样率 └── fig/ % 出图保存目录这个拆分方式对应的是标准 EEMD 的三层结构。最底层是emd.m它只做一件事给定一段信号找出所有局部极大值和极小值用三次样条插出上、下包络算出均值曲线原信号减均值得到一阶 IMF 候选然后对残差重复这一过程直到剩余分量单调或层数达到上限。中间层eemd.m不关心筛分细节它只负责循环加一组白噪声、调一次emd、把结果按层对齐并累加跑完NE次后除以NE。最上层main_eemd.m才是面向使用者的负责读数据、设参数、选分量、画包络谱。函数边界清不清楚直接决定二次开发时的修改成本。我见过不少把加噪循环和样条插值写在同一段脚本里的版本改个噪声幅值要翻两百行这种就不适合做实验对比。分层之后想把 EEMD 换成 CEEMDAN理论上只需要换掉eemd.m这一层main_eemd.m和后面的包络分析一行都不用动。2.2 为什么是 EEMD 而不是 EMD 或小波轴承振动信号是典型的非平稳信号冲击发生时频率成分在毫秒级内突变。小波要预先选定基函数和分解层数基函数选得不对故障频带就被切到两个尺度之间重构出来能量散掉。EMD 是数据驱动的分解基由信号自身的极值点决定不需要预设基函数这一点对冲击类信号很友好。但裸 EMD 有个绕不过去的问题模态混叠。当信号里同时存在高频冲击和低频转频时同一个 IMF 里会混进差异很大的时间尺度表现为某个分量一会儿是冲击包络、一会儿是转频正弦后续做包络谱时特征频率被调制得认不出来。EEMD 的补丁思路很直接——往信号里加白噪声让极值点分布更均匀噪声把原本会被“跳过”的小尺度极值点激活出来单次分解的混叠程度就下降独立重复NE次每次噪声不同最后平均噪声本身因为零均值被抵消大半而 IMF 的真实成分因为每次都能复现被保留下来。代价是计算量。EMD 一次分解已经要反复做样条插值EEMD 要跑NE次NE取 100 是常规值信号长度一万点、分解 10 层在普通笔记本上跑几分钟到十几分钟很正常。包里如果带示例数据第一次运行建议先把NE调小到 10 跑通流程确认数据路径和绘图没问题再改回 100 做正式分析。这是我自己踩过坑之后养成的习惯拿新代码先降参数跑通比一上来就全量跑然后卡在某个报错上强得多。3. 跑通第一遍从解压到画出 IMF 图3.1 环境准备与路径管理MATLAB 版本上没有太苛刻的要求R2016a 之后基本都能跑因为用到的hilbert、spline、fft、corr都是老函数没有依赖新的工具箱接口。信号处理工具箱建议装上hilbert在部分精简安装里会被拆出去。解压之后第一件事是把当前工作目录切到eemd_bearing/根目录或者用addpath(genpath(pwd))把子目录一次性加进来否则main_eemd.m里调utils/下的函数会报未定义。% 在 MATLAB 命令行执行把工程根目录及其子目录加入搜索路径 root D:\work\eemd_bearing; % 换成自己的解压路径 addpath(genpath(root)); cd(root); savepath; % 可选把路径存下来下次打开不用重复加genpath会递归包含所有子文件夹data/和fig/也会被加进去这不影响运行但如果工程里混进了.git或者备份目录加载会变慢路径也容易冲突。我一般会手动指定utils一个子目录而不是无脑genpath。savepath会把当前路径写进pathdef.m换项目之后记得清掉不然下次打开 MATLAB 会看到一堆黄色警告说路径不存在。3.2 数据读取与采样参数包里的示例数据通常是.mat格式里面存两个变量振动信号列向量和采样频率。有些版本会把采样率硬编码在main_eemd.m顶部有些放在数据文件的元信息里。不管哪种跑之前必须把fs和信号长度核对一遍因为后面包络谱的横轴完全依赖fs这个值填错故障特征频率的定位就会整体平移看上去像“算法不准”其实是参数没对上。% 读取示例数据兼容 .mat 与纯文本两种常见格式 function [x, fs] read_signal(filepath, fs_default) [~, ~, ext] fileparts(filepath); switch lower(ext) case .mat S load(filepath); % 变量名不固定时按字段名猜含 sig / x / data 的优先 fn fieldnames(S); x []; for i 1:numel(fn) v S.(fn{i}); if isnumeric(v) numel(v) 1000 isvector(v) x v(:); break; end end fs fs_default; case {.txt, .csv} x load(filepath); % 单列数值直接 load x x(:); fs fs_default; otherwise error(不支持的信号文件格式: %s, ext); end x detrend(x); % 去线性趋势避免低频漂移污染末阶 IMF end这段函数做三件事按扩展名分流读取、把信号强制成列向量、去趋势。detrend这一步容易被忽略但轴承信号里的直流偏置和缓慢温漂会在 EMD 里被当成一个超低频分量占掉最后一两阶 IMF导致真正的高频冲击分量层数被往后挤。fs_default作为兜底参数传进来是因为.txt文件里没有地方存采样率只能从主脚本传入。如果读者的数据是 DAQ 导出的带表头 CSVload会失败需要换成readmatrix并跳过前几行这个改动放在read_signal里就行不用动主流程。3.3 主脚本执行链路main_eemd.m从上到下的执行顺序是固定的清空环境、设参数、读数据、调 EEMD、画 IMF 时域图、选分量、算包络谱。跑通的标准不是“没报错”而是能看到两张图一张是 8 到 10 阶 IMF 的时域波形堆叠图一张是选定分量对应的包络谱谱上能看到一个明显高于周围底噪的峰。%% main_eemd.m 关键段落 clear; clc; close all; fs 12000; % 采样频率Hz必须和采集设置一致 Nstd 0.2; % 白噪声标准差与信号标准差的比值 NE 100; % 集成次数 maxIMF 10; % 最大分解层数 [x, fs] read_signal(fullfile(data,bearing_inner.mat), fs); % EEMD 分解imf 每列是一阶分量res 是残余趋势项 [imf, res] eemd(x, Nstd, NE, maxIMF); % 画 IMF 时域堆叠图 figure(Name,IMF 分量); nIMF size(imf, 2); for i 1:nIMF subplot(nIMF,1,i); plot((0:numel(x)-1)/fs, imf(:,i)); ylabel(sprintf(IMF%d, i)); if i nIMF, xlabel(时间 / s); end end % 选分量 包络谱 idx select_imf(x, imf, kurtosis); [f, P] envelope_spectrum(imf(:,idx), fs); figure(Name,包络谱); plot(f, P); xlabel(频率 / Hz); ylabel(幅值);第一次跑建议把NE改成 10、maxIMF改成 6几分钟内出结果确认data/路径对、read_signal返回的向量长度合理、eemd没有因为层数不足报维度不匹配。确认链路通了再恢复正式参数。select_imf的第三个参数是准则名包内如果支持corr和kurtosis两种先用kurtosis因为轴承故障冲击在峭度上的表现比相关系数更敏感尤其在内圈故障早期。4. 核心算法拆解加噪集成、IMF 筛选与包络谱4.1 EEMD 主循环白噪声怎么加、怎么平均eemd.m是整个包的心脏逻辑不复杂但几个参数的含义和边界处理容易理解偏。核心循环是每次迭代生成一组与信号等长的高斯白噪声按Nstd * std(x)缩放后叠加到原信号上调一次 EMD 得到该次的 IMF 集合把这组 IMF 按层累加最后除以集成次数。function [imf, res] eemd(x, Nstd, NE, maxIMF) % EEMD 集成经验模态分解 % 输入: % x 原始信号列向量 % Nstd 噪声标准差比例常用 0.1~0.3太大污染低频太小压不住混叠 % NE 集成次数常用 100~200越大越平滑耗时线性增长 % maxIMF 最大 IMF 层数 % 输出: % imf N x K 矩阵每列一阶 IMF按频率从高到低排列 % res 残余趋势项 x x(:); N numel(x); sig std(x); rng(20240101); % 固定随机种子保证结果可复现 imfAcc []; % 累加器 for k 1:NE noise Nstd * sig * randn(N, 1); imfK emd_core(x noise, maxIMF); % 单次 EMD K size(imfK, 2); % 首次迭代初始化累加矩阵后续层数不一致时补零对齐 if isempty(imfAcc) imfAcc zeros(N, K); elseif K size(imfAcc, 2) imfAcc [imfAcc, zeros(N, K - size(imfAcc,2))]; elseif K size(imfAcc, 2) imfK [imfK, zeros(N, size(imfAcc,2) - K)]; end imfAcc imfAcc imfK; end imf imfAcc / NE; res x - sum(imf, 2); end有三处细节值得单独说。第一是rng。EEMD 本身带随机噪声不固定种子的话同一段信号跑两次得到的 IMF 会有可见差异写报告或做对比实验时没法复现这是最容易被忽略的一处“玄学”。固定种子之后结果完全确定代价是失去了随机性带来的统计意义如果要做多次独立实验取均值就在外层再套一层改变种子的循环。第二是层数对齐。emd_core每次分解得到的层数不一定相同信号端点处极值点分布变化会导致某次多分出一层。直接矩阵拼接会报维度不一致所以要先判断当前累加矩阵和目标矩阵的列数少的补零。补零的位置在右侧也就是低频端因为层数差异几乎总是出现在最后几阶补零不会影响高频 IMF 的平均。第三是Nstd的取值。0.2 是文献里的常规值但对于冲击成分占比很低的早期故障信号噪声幅值太大会把微弱冲击直接淹掉反过来Nstd降到 0.05 以下噪声不足以激活均匀分布的极值点混叠又回来了。我的经验是先扫一遍 0.1、0.2、0.3 三档看包络谱上故障频率峰的信噪比选最高的那一档不要迷信固定值。4.2 分量挑选相关系数还是峭度EEMD 分解完得到一堆 IMF真正和轴承故障相关的是其中一两阶。选错分量后面的包络谱再漂亮也是噪声。包内select_imf.m一般支持两种准则逻辑差异很大。相关系数准则算的是每阶 IMF 与原信号的皮尔逊相关系数取最大的那一阶。这个准则的假设是“与原始信号最像的分量包含主要能量”对转频主导的信号有效但轴承故障冲击能量占比低最相关的往往是低频转频分量选出来做包络谱只能看到转频谐波。峭度准则算的是每阶 IMF 的峭度值取最大的那一阶依据是冲击成分会让分布尾部变厚、峭度升高。内圈故障早期冲击稀疏但幅值高峭度对这类信号更敏感。function idx select_imf(x, imf, mode) % 按准则挑选最可能承载故障冲击的 IMF % mode: corr 相关系数最大kurtosis 峭度最大 K size(imf, 2); switch lower(mode) case corr score zeros(1, K); for i 1:K c corrcoef(x, imf(:,i)); score(i) abs(c(1,2)); end case kurtosis score kurtosis(imf); % 每列一个峭度值 otherwise error(未知筛选准则: %s, mode); end [~, idx] max(score); end实际使用中更稳的做法是两级筛选先按峭度排序取前 3 阶再在这 3 阶里看谁的包络谱在理论故障频率附近有峰。只用单一指标容易翻车尤其当信号里混进了电磁干扰或传感器松动产生的瞬态冲击时峭度最高的那一阶可能对应的是干扰而不是轴承故障。这种情况在台架数据上不明显现场数据上非常常见所以选完分量最好把包络谱和理论频率表对照一眼再下结论。4.3 Hilbert 包络与故障特征频率定位选定 IMF 之后做包络解调。轴承故障冲击会激发系统的高频固有频率冲击以固有频率为载波、以故障特征频率为重复率出现所以对高频 IMF 做 Hilbert 变换取模得到的就是冲击重复率的包络再对包络做 FFT故障特征频率就落到低频段和转频谐波分开。function [f, P] envelope_spectrum(x, fs) % 对单分量信号做 Hilbert 包络谱 x x(:); env abs(hilbert(x)); % 解析信号取模 包络 env env - mean(env); % 去直流否则 0 Hz 处一根巨峰压掉一切 L numel(env); Y fft(env); P abs(Y(1:floor(L/2)1)) / L * 2; % 单边幅值谱 f (0:floor(L/2)) * fs / L; endenv - mean(env)这行不能省包络本身是正值序列均值不为零FFT 之后 0 Hz 处的幅值会比故障频率峰高一到两个数量级画图时纵轴被拉爆故障峰看起来像贴着地面。去均值之后横轴从 0 开始纵轴量级正常故障频率峰才看得清。理论特征频率得自己算包里通常不带轴承参数需要按型号查滚动体直径、节圆直径、接触角和滚动体个数function [BPFO, BPFI, BSF, FTF] bearing_freqs(d, D, n, fr, alpha) % d 滚动体直径, D 节圆直径, n 滚动体个数, fr 轴转频(Hz), alpha 接触角(rad) r d / D * cos(alpha); BPFO n / 2 * fr * (1 - r); % 外圈故障特征频率 BPFI n / 2 * fr * (1 r); % 内圈故障特征频率 BSF D / (2 * d) * fr * (1 - r^2); % 滚动体故障特征频率 FTF fr / 2 * (1 - r); % 保持架故障特征频率 end算出来之后在包络谱上找BPFO或BPFI的 1 倍、2 倍、3 倍频内圈故障还会出现以转频为间隔的边频带因为内圈随轴旋转故障点进出载荷区时冲击幅值被转频调制。看到“特征频率 转频间隔边带”这个组合基本可以确认是内圈故障。只看到单根峰、没有边带更可能是外圈或干扰。5. 常见问题排查跑不出结果的几个真实现场5.1 每次运行结果都不一样现象是同一份数据、同一组参数跑两次得到的 IMF 波形肉眼可见地不同包络谱峰值位置也跟着跳。原因是 EEMD 内部用了randn生成白噪声MATLAB 默认的随机数生成器在每次启动时按时间种子初始化不固定就不复现。解决办法是在调用eemd之前加一行rng(固定整数)或者把种子作为参数传进eemd.m内部。固定之后结果完全确定方便写对比实验和出图。如果做统计性分析需要多次独立实验就在外层循环里改变种子把每次的包络谱峰值记录下来取均值而不是让内部随机。5.2 包络谱上找不到故障特征频率现象是 IMF 图看着正常、包络谱也出了图但理论BPFI位置没有明显峰全是低幅值起伏。原因通常有三类fs填错导致横轴整体缩放这个最常见先核对采样率选的 IMF 不对峭度最大那阶可能对应传感器共振而不是故障冲击信号的故障频带不在所选 IMF 覆盖的范围内。解决办法是按顺序排查先用已知频率的正弦信号做一次自检确认fs链路没问题再打印每阶 IMF 的峭度和相关系数手动切换分量看包络谱变化最后考虑在 EEMD 之前加一级带通滤波把故障激发的固有频率频带先框出来。5.3 跑到一半内存溢出或耗时失控现象是NE100、信号长度十几万点时 MATLAB 卡死任务管理器里内存一路涨到几个 G或者跑了一晚上没出结果。原因是eemd.m里每次迭代都生成一个N x K的 IMF 矩阵K取 10 层、N取 20 万点单次就是 1600 万浮点数累加器再加一份NE循环里还有临时变量没及时释放。解决办法是分块处理长信号按每段 2 万到 5 万点切分、重叠 50% 做分解再拼接或者把maxIMF降到 6 到 8 层因为轴承故障分析真正有用的就是前几阶高频分量后面几阶基本是趋势项。NE也可以从 100 降到 50模态混叠的抑制效果差别不大耗时直接减半。5.4 IMF 两端出现发散振荡现象是每阶 IMF 的起始和结束几十个采样点幅值明显放大像被拉出去一样越靠后的分量越严重。原因是 EMD 用三次样条拟合上下包络时端点处缺少极值点约束样条会向外 extrapolate这就是常说的端点效应EEMD 加噪声只能在统计意义上缓解不能根除。解决办法是在分解前对信号做镜像延拓把首尾各一段数据翻转接到两端分解完再裁掉或者用包内如果自带的极值点延拓函数替换envelope.m里的边界处理。要求不高的话直接裁掉两端各 5% 的数据也能用代价是损失一点样本。5.5 换一份数据就报维度错误现象是示例数据跑得好好的换成自己的 CSV 之后eemd报矩阵维度不匹配。原因是read_signal里用load读单列文本时如果文件带表头或者有多列时间戳 幅值返回的是二维矩阵直接传进eemd后randn(N,1)和N x 1的加法对不上。解决办法是在read_signal里显式处理用readmatrix指定列范围跳过表头行取幅值那一列并(:)拉平。养成在读取之后打印size(x)和std(x)的习惯长度和量级不对后面全是白费功夫。6. 二次开发把单文件脚本改成批量诊断工具单个文件跑通只是起点真要做一组台架实验的对比分析得把这份代码改成能批量吃数据、自动出特征、结果落到表里的工具。我的改法是把main_eemd.m拆成两层一个run_one(filepath, cfg)负责单文件全流程并返回结构化结果一个batch_run(filelist, cfg)循环调用并汇总。function R run_one(filepath, cfg) % 单文件诊断返回选定 IMF、包络谱峰值及对应频率 [x, fs] read_signal(filepath, cfg.fs); [imf, ~] eemd(x, cfg.Nstd, cfg.NE, cfg.maxIMF); idx select_imf(x, imf, cfg.imfMode); [f, P] envelope_spectrum(imf(:,idx), fs); % 在理论特征频率 ±2% 邻域内找峰 band abs(f - cfg.fTarget) 0.02 * cfg.fTarget; [pk, loc] max(P(band)); fAxis f(band); R.file filepath; R.imfIdx idx; R.peakF fAxis(loc); R.peakA pk; R.snr pk / median(P); % 简易信噪比用于横向对比 end这里新增的snr字段是横向对比的关键单看峰值幅值没有意义不同文件之间增益和量程可能不同用峰值除以全谱中位数得到的相对值才有可比性。batch_run把R收集成 struct 数组再用struct2table转表导出 CSV一行一个样本按snr排序就能快速定位哪组数据的故障特征最明显。参数扫描是另一个高频需求。Nstd和NE到底取多少靠一次实验定不下来写一个双重循环把Nstd从 0.05 扫到 0.4、NE从 20 扫到 200每个组合记录snr画热图峰值所在区域就是这组数据的合理参数区间。这套流程跑一遍要几十分钟但比反复手动改参数试错省心得多。验证改动有没有效果我的习惯是保留一份原始示例数据的基准结果固定的rng种子、固定的参数、固定的输入文件每次改完代码先跑基准imf矩阵和peakF与上一次完全一致才算没有破坏原有逻辑。EMD 这类迭代算法对数值扰动很敏感改一行边界插值就可能让整条 IMF 序列偏移没有基准对照根本发现不了。从那以后我每次拿到这类分解算法的包都先把NE降下来跑一遍基准、把结果存成.mat锁死再动任何一行代码都跟基准对一遍——这条规矩让我少返工很多次。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/26 5:08:29 · 阅读 1.2k 基于MATLAB的EEMD轴承故障诊断与包络谱特征提取
简介围绕EEMD集成经验模态分解的MATLAB程序包面向从事轴承故障诊断、机械振动信号分析的研究人员与工程技术人员也可供学习非线性、非平稳信号处理的学生参考。它针对传统EMD易出现模式混叠、噪声干扰等问题通过白噪声辅助与多次分解平均提升IMF提取稳定性便于从振动信号中分离异常频率成分并开展故障特征分析。压缩包共11个文件以9个m脚本为主体另有1个mat数据文件和1个asv备份整体约761KB体积轻便适合直接放入MATLAB环境调试与二次修改。已有2034人学习说明其在相关学习中具有一定参考价值。使用者可据此加载自己的振动数据完成EEMD分解、IMF分量筛选与结果可视化进而理解噪声辅助分解流程、参数设置和轴承故障特征提取思路为设备状态监测与早期故障识别提供可复用的代码基础。1. 从一段被噪声淹没的轴承振动信号说起轴承早期故障的冲击成分非常弱一段 12 kHz 采样、上万点的加速度信号直接丢进 FFT谱线上全是转频谐波和背景噪声故障特征频率根本冒不出来。EEMD集成经验模态分解处理的就是这种非平稳、低信噪比的场景往原始信号里多次叠加不同幅值的白噪声再做 EMD利用噪声在频谱上均匀分布的特性把模态混叠压下去最后对多次分解结果取集成平均得到一组按时间尺度排列的 IMF 分量。挑出承载故障冲击的高频 IMF 做 Hilbert 包络 demodulation包络谱上才会出现清晰的 BPFO、BPFI 及其倍频。这份 MATLAB 程序把这条链路从数据读入、EEMD 分解、IMF 筛选到包络谱出图串成了可运行的原型适合手里有台架或现场振动数据、想做轴承故障特征提取验证的开发者也适合刚接触经验模态分解、需要一份能跑起来再动手改的代码作参照。2. 程序结构与技术选型一个 rar 包里装的是什么2.1 目录骨架与函数分工rаr 解压之后正常的工程目录不会把所有东西塞进一个脚本文件里而是按“分解核心 / 辅助函数 / 数据 / 出图”拆开。我拆这类包最先看的是eemd.m和emd.m的调用关系因为它决定了整套代码是自研实现还是套用现成的 EMD 内核。常见的目录形态大概是这样eemd_bearing/ ├── main_eemd.m % 主入口读数据、调分解、画图 ├── eemd.m % 集成经验模态分解负责加噪与平均 ├── emd.m % 单次 EMD 筛分过程 ├── extrema.m % 极值点定位局部极大/极小 ├── envelope.m % 三次样条上下包络插值 ├── utils/ │ ├── read_signal.m % 读取 txt / csv / mat 振动数据 │ ├── select_imf.m % IMF 筛选相关系数或峭度准则 │ └── envelope_spectrum.m % Hilbert 包络 频谱 ├── data/ │ └── bearing_inner.mat % 示例数据内含信号与采样率 └── fig/ % 出图保存目录这个拆分方式对应的是标准 EEMD 的三层结构。最底层是emd.m它只做一件事给定一段信号找出所有局部极大值和极小值用三次样条插出上、下包络算出均值曲线原信号减均值得到一阶 IMF 候选然后对残差重复这一过程直到剩余分量单调或层数达到上限。中间层eemd.m不关心筛分细节它只负责循环加一组白噪声、调一次emd、把结果按层对齐并累加跑完NE次后除以NE。最上层main_eemd.m才是面向使用者的负责读数据、设参数、选分量、画包络谱。函数边界清不清楚直接决定二次开发时的修改成本。我见过不少把加噪循环和样条插值写在同一段脚本里的版本改个噪声幅值要翻两百行这种就不适合做实验对比。分层之后想把 EEMD 换成 CEEMDAN理论上只需要换掉eemd.m这一层main_eemd.m和后面的包络分析一行都不用动。2.2 为什么是 EEMD 而不是 EMD 或小波轴承振动信号是典型的非平稳信号冲击发生时频率成分在毫秒级内突变。小波要预先选定基函数和分解层数基函数选得不对故障频带就被切到两个尺度之间重构出来能量散掉。EMD 是数据驱动的分解基由信号自身的极值点决定不需要预设基函数这一点对冲击类信号很友好。但裸 EMD 有个绕不过去的问题模态混叠。当信号里同时存在高频冲击和低频转频时同一个 IMF 里会混进差异很大的时间尺度表现为某个分量一会儿是冲击包络、一会儿是转频正弦后续做包络谱时特征频率被调制得认不出来。EEMD 的补丁思路很直接——往信号里加白噪声让极值点分布更均匀噪声把原本会被“跳过”的小尺度极值点激活出来单次分解的混叠程度就下降独立重复NE次每次噪声不同最后平均噪声本身因为零均值被抵消大半而 IMF 的真实成分因为每次都能复现被保留下来。代价是计算量。EMD 一次分解已经要反复做样条插值EEMD 要跑NE次NE取 100 是常规值信号长度一万点、分解 10 层在普通笔记本上跑几分钟到十几分钟很正常。包里如果带示例数据第一次运行建议先把NE调小到 10 跑通流程确认数据路径和绘图没问题再改回 100 做正式分析。这是我自己踩过坑之后养成的习惯拿新代码先降参数跑通比一上来就全量跑然后卡在某个报错上强得多。3. 跑通第一遍从解压到画出 IMF 图3.1 环境准备与路径管理MATLAB 版本上没有太苛刻的要求R2016a 之后基本都能跑因为用到的hilbert、spline、fft、corr都是老函数没有依赖新的工具箱接口。信号处理工具箱建议装上hilbert在部分精简安装里会被拆出去。解压之后第一件事是把当前工作目录切到eemd_bearing/根目录或者用addpath(genpath(pwd))把子目录一次性加进来否则main_eemd.m里调utils/下的函数会报未定义。% 在 MATLAB 命令行执行把工程根目录及其子目录加入搜索路径 root D:\work\eemd_bearing; % 换成自己的解压路径 addpath(genpath(root)); cd(root); savepath; % 可选把路径存下来下次打开不用重复加genpath会递归包含所有子文件夹data/和fig/也会被加进去这不影响运行但如果工程里混进了.git或者备份目录加载会变慢路径也容易冲突。我一般会手动指定utils一个子目录而不是无脑genpath。savepath会把当前路径写进pathdef.m换项目之后记得清掉不然下次打开 MATLAB 会看到一堆黄色警告说路径不存在。3.2 数据读取与采样参数包里的示例数据通常是.mat格式里面存两个变量振动信号列向量和采样频率。有些版本会把采样率硬编码在main_eemd.m顶部有些放在数据文件的元信息里。不管哪种跑之前必须把fs和信号长度核对一遍因为后面包络谱的横轴完全依赖fs这个值填错故障特征频率的定位就会整体平移看上去像“算法不准”其实是参数没对上。% 读取示例数据兼容 .mat 与纯文本两种常见格式 function [x, fs] read_signal(filepath, fs_default) [~, ~, ext] fileparts(filepath); switch lower(ext) case .mat S load(filepath); % 变量名不固定时按字段名猜含 sig / x / data 的优先 fn fieldnames(S); x []; for i 1:numel(fn) v S.(fn{i}); if isnumeric(v) numel(v) 1000 isvector(v) x v(:); break; end end fs fs_default; case {.txt, .csv} x load(filepath); % 单列数值直接 load x x(:); fs fs_default; otherwise error(不支持的信号文件格式: %s, ext); end x detrend(x); % 去线性趋势避免低频漂移污染末阶 IMF end这段函数做三件事按扩展名分流读取、把信号强制成列向量、去趋势。detrend这一步容易被忽略但轴承信号里的直流偏置和缓慢温漂会在 EMD 里被当成一个超低频分量占掉最后一两阶 IMF导致真正的高频冲击分量层数被往后挤。fs_default作为兜底参数传进来是因为.txt文件里没有地方存采样率只能从主脚本传入。如果读者的数据是 DAQ 导出的带表头 CSVload会失败需要换成readmatrix并跳过前几行这个改动放在read_signal里就行不用动主流程。3.3 主脚本执行链路main_eemd.m从上到下的执行顺序是固定的清空环境、设参数、读数据、调 EEMD、画 IMF 时域图、选分量、算包络谱。跑通的标准不是“没报错”而是能看到两张图一张是 8 到 10 阶 IMF 的时域波形堆叠图一张是选定分量对应的包络谱谱上能看到一个明显高于周围底噪的峰。%% main_eemd.m 关键段落 clear; clc; close all; fs 12000; % 采样频率Hz必须和采集设置一致 Nstd 0.2; % 白噪声标准差与信号标准差的比值 NE 100; % 集成次数 maxIMF 10; % 最大分解层数 [x, fs] read_signal(fullfile(data,bearing_inner.mat), fs); % EEMD 分解imf 每列是一阶分量res 是残余趋势项 [imf, res] eemd(x, Nstd, NE, maxIMF); % 画 IMF 时域堆叠图 figure(Name,IMF 分量); nIMF size(imf, 2); for i 1:nIMF subplot(nIMF,1,i); plot((0:numel(x)-1)/fs, imf(:,i)); ylabel(sprintf(IMF%d, i)); if i nIMF, xlabel(时间 / s); end end % 选分量 包络谱 idx select_imf(x, imf, kurtosis); [f, P] envelope_spectrum(imf(:,idx), fs); figure(Name,包络谱); plot(f, P); xlabel(频率 / Hz); ylabel(幅值);第一次跑建议把NE改成 10、maxIMF改成 6几分钟内出结果确认data/路径对、read_signal返回的向量长度合理、eemd没有因为层数不足报维度不匹配。确认链路通了再恢复正式参数。select_imf的第三个参数是准则名包内如果支持corr和kurtosis两种先用kurtosis因为轴承故障冲击在峭度上的表现比相关系数更敏感尤其在内圈故障早期。4. 核心算法拆解加噪集成、IMF 筛选与包络谱4.1 EEMD 主循环白噪声怎么加、怎么平均eemd.m是整个包的心脏逻辑不复杂但几个参数的含义和边界处理容易理解偏。核心循环是每次迭代生成一组与信号等长的高斯白噪声按Nstd * std(x)缩放后叠加到原信号上调一次 EMD 得到该次的 IMF 集合把这组 IMF 按层累加最后除以集成次数。function [imf, res] eemd(x, Nstd, NE, maxIMF) % EEMD 集成经验模态分解 % 输入: % x 原始信号列向量 % Nstd 噪声标准差比例常用 0.1~0.3太大污染低频太小压不住混叠 % NE 集成次数常用 100~200越大越平滑耗时线性增长 % maxIMF 最大 IMF 层数 % 输出: % imf N x K 矩阵每列一阶 IMF按频率从高到低排列 % res 残余趋势项 x x(:); N numel(x); sig std(x); rng(20240101); % 固定随机种子保证结果可复现 imfAcc []; % 累加器 for k 1:NE noise Nstd * sig * randn(N, 1); imfK emd_core(x noise, maxIMF); % 单次 EMD K size(imfK, 2); % 首次迭代初始化累加矩阵后续层数不一致时补零对齐 if isempty(imfAcc) imfAcc zeros(N, K); elseif K size(imfAcc, 2) imfAcc [imfAcc, zeros(N, K - size(imfAcc,2))]; elseif K size(imfAcc, 2) imfK [imfK, zeros(N, size(imfAcc,2) - K)]; end imfAcc imfAcc imfK; end imf imfAcc / NE; res x - sum(imf, 2); end有三处细节值得单独说。第一是rng。EEMD 本身带随机噪声不固定种子的话同一段信号跑两次得到的 IMF 会有可见差异写报告或做对比实验时没法复现这是最容易被忽略的一处“玄学”。固定种子之后结果完全确定代价是失去了随机性带来的统计意义如果要做多次独立实验取均值就在外层再套一层改变种子的循环。第二是层数对齐。emd_core每次分解得到的层数不一定相同信号端点处极值点分布变化会导致某次多分出一层。直接矩阵拼接会报维度不一致所以要先判断当前累加矩阵和目标矩阵的列数少的补零。补零的位置在右侧也就是低频端因为层数差异几乎总是出现在最后几阶补零不会影响高频 IMF 的平均。第三是Nstd的取值。0.2 是文献里的常规值但对于冲击成分占比很低的早期故障信号噪声幅值太大会把微弱冲击直接淹掉反过来Nstd降到 0.05 以下噪声不足以激活均匀分布的极值点混叠又回来了。我的经验是先扫一遍 0.1、0.2、0.3 三档看包络谱上故障频率峰的信噪比选最高的那一档不要迷信固定值。4.2 分量挑选相关系数还是峭度EEMD 分解完得到一堆 IMF真正和轴承故障相关的是其中一两阶。选错分量后面的包络谱再漂亮也是噪声。包内select_imf.m一般支持两种准则逻辑差异很大。相关系数准则算的是每阶 IMF 与原信号的皮尔逊相关系数取最大的那一阶。这个准则的假设是“与原始信号最像的分量包含主要能量”对转频主导的信号有效但轴承故障冲击能量占比低最相关的往往是低频转频分量选出来做包络谱只能看到转频谐波。峭度准则算的是每阶 IMF 的峭度值取最大的那一阶依据是冲击成分会让分布尾部变厚、峭度升高。内圈故障早期冲击稀疏但幅值高峭度对这类信号更敏感。function idx select_imf(x, imf, mode) % 按准则挑选最可能承载故障冲击的 IMF % mode: corr 相关系数最大kurtosis 峭度最大 K size(imf, 2); switch lower(mode) case corr score zeros(1, K); for i 1:K c corrcoef(x, imf(:,i)); score(i) abs(c(1,2)); end case kurtosis score kurtosis(imf); % 每列一个峭度值 otherwise error(未知筛选准则: %s, mode); end [~, idx] max(score); end实际使用中更稳的做法是两级筛选先按峭度排序取前 3 阶再在这 3 阶里看谁的包络谱在理论故障频率附近有峰。只用单一指标容易翻车尤其当信号里混进了电磁干扰或传感器松动产生的瞬态冲击时峭度最高的那一阶可能对应的是干扰而不是轴承故障。这种情况在台架数据上不明显现场数据上非常常见所以选完分量最好把包络谱和理论频率表对照一眼再下结论。4.3 Hilbert 包络与故障特征频率定位选定 IMF 之后做包络解调。轴承故障冲击会激发系统的高频固有频率冲击以固有频率为载波、以故障特征频率为重复率出现所以对高频 IMF 做 Hilbert 变换取模得到的就是冲击重复率的包络再对包络做 FFT故障特征频率就落到低频段和转频谐波分开。function [f, P] envelope_spectrum(x, fs) % 对单分量信号做 Hilbert 包络谱 x x(:); env abs(hilbert(x)); % 解析信号取模 包络 env env - mean(env); % 去直流否则 0 Hz 处一根巨峰压掉一切 L numel(env); Y fft(env); P abs(Y(1:floor(L/2)1)) / L * 2; % 单边幅值谱 f (0:floor(L/2)) * fs / L; endenv - mean(env)这行不能省包络本身是正值序列均值不为零FFT 之后 0 Hz 处的幅值会比故障频率峰高一到两个数量级画图时纵轴被拉爆故障峰看起来像贴着地面。去均值之后横轴从 0 开始纵轴量级正常故障频率峰才看得清。理论特征频率得自己算包里通常不带轴承参数需要按型号查滚动体直径、节圆直径、接触角和滚动体个数function [BPFO, BPFI, BSF, FTF] bearing_freqs(d, D, n, fr, alpha) % d 滚动体直径, D 节圆直径, n 滚动体个数, fr 轴转频(Hz), alpha 接触角(rad) r d / D * cos(alpha); BPFO n / 2 * fr * (1 - r); % 外圈故障特征频率 BPFI n / 2 * fr * (1 r); % 内圈故障特征频率 BSF D / (2 * d) * fr * (1 - r^2); % 滚动体故障特征频率 FTF fr / 2 * (1 - r); % 保持架故障特征频率 end算出来之后在包络谱上找BPFO或BPFI的 1 倍、2 倍、3 倍频内圈故障还会出现以转频为间隔的边频带因为内圈随轴旋转故障点进出载荷区时冲击幅值被转频调制。看到“特征频率 转频间隔边带”这个组合基本可以确认是内圈故障。只看到单根峰、没有边带更可能是外圈或干扰。5. 常见问题排查跑不出结果的几个真实现场5.1 每次运行结果都不一样现象是同一份数据、同一组参数跑两次得到的 IMF 波形肉眼可见地不同包络谱峰值位置也跟着跳。原因是 EEMD 内部用了randn生成白噪声MATLAB 默认的随机数生成器在每次启动时按时间种子初始化不固定就不复现。解决办法是在调用eemd之前加一行rng(固定整数)或者把种子作为参数传进eemd.m内部。固定之后结果完全确定方便写对比实验和出图。如果做统计性分析需要多次独立实验就在外层循环里改变种子把每次的包络谱峰值记录下来取均值而不是让内部随机。5.2 包络谱上找不到故障特征频率现象是 IMF 图看着正常、包络谱也出了图但理论BPFI位置没有明显峰全是低幅值起伏。原因通常有三类fs填错导致横轴整体缩放这个最常见先核对采样率选的 IMF 不对峭度最大那阶可能对应传感器共振而不是故障冲击信号的故障频带不在所选 IMF 覆盖的范围内。解决办法是按顺序排查先用已知频率的正弦信号做一次自检确认fs链路没问题再打印每阶 IMF 的峭度和相关系数手动切换分量看包络谱变化最后考虑在 EEMD 之前加一级带通滤波把故障激发的固有频率频带先框出来。5.3 跑到一半内存溢出或耗时失控现象是NE100、信号长度十几万点时 MATLAB 卡死任务管理器里内存一路涨到几个 G或者跑了一晚上没出结果。原因是eemd.m里每次迭代都生成一个N x K的 IMF 矩阵K取 10 层、N取 20 万点单次就是 1600 万浮点数累加器再加一份NE循环里还有临时变量没及时释放。解决办法是分块处理长信号按每段 2 万到 5 万点切分、重叠 50% 做分解再拼接或者把maxIMF降到 6 到 8 层因为轴承故障分析真正有用的就是前几阶高频分量后面几阶基本是趋势项。NE也可以从 100 降到 50模态混叠的抑制效果差别不大耗时直接减半。5.4 IMF 两端出现发散振荡现象是每阶 IMF 的起始和结束几十个采样点幅值明显放大像被拉出去一样越靠后的分量越严重。原因是 EMD 用三次样条拟合上下包络时端点处缺少极值点约束样条会向外 extrapolate这就是常说的端点效应EEMD 加噪声只能在统计意义上缓解不能根除。解决办法是在分解前对信号做镜像延拓把首尾各一段数据翻转接到两端分解完再裁掉或者用包内如果自带的极值点延拓函数替换envelope.m里的边界处理。要求不高的话直接裁掉两端各 5% 的数据也能用代价是损失一点样本。5.5 换一份数据就报维度错误现象是示例数据跑得好好的换成自己的 CSV 之后eemd报矩阵维度不匹配。原因是read_signal里用load读单列文本时如果文件带表头或者有多列时间戳 幅值返回的是二维矩阵直接传进eemd后randn(N,1)和N x 1的加法对不上。解决办法是在read_signal里显式处理用readmatrix指定列范围跳过表头行取幅值那一列并(:)拉平。养成在读取之后打印size(x)和std(x)的习惯长度和量级不对后面全是白费功夫。6. 二次开发把单文件脚本改成批量诊断工具单个文件跑通只是起点真要做一组台架实验的对比分析得把这份代码改成能批量吃数据、自动出特征、结果落到表里的工具。我的改法是把main_eemd.m拆成两层一个run_one(filepath, cfg)负责单文件全流程并返回结构化结果一个batch_run(filelist, cfg)循环调用并汇总。function R run_one(filepath, cfg) % 单文件诊断返回选定 IMF、包络谱峰值及对应频率 [x, fs] read_signal(filepath, cfg.fs); [imf, ~] eemd(x, cfg.Nstd, cfg.NE, cfg.maxIMF); idx select_imf(x, imf, cfg.imfMode); [f, P] envelope_spectrum(imf(:,idx), fs); % 在理论特征频率 ±2% 邻域内找峰 band abs(f - cfg.fTarget) 0.02 * cfg.fTarget; [pk, loc] max(P(band)); fAxis f(band); R.file filepath; R.imfIdx idx; R.peakF fAxis(loc); R.peakA pk; R.snr pk / median(P); % 简易信噪比用于横向对比 end这里新增的snr字段是横向对比的关键单看峰值幅值没有意义不同文件之间增益和量程可能不同用峰值除以全谱中位数得到的相对值才有可比性。batch_run把R收集成 struct 数组再用struct2table转表导出 CSV一行一个样本按snr排序就能快速定位哪组数据的故障特征最明显。参数扫描是另一个高频需求。Nstd和NE到底取多少靠一次实验定不下来写一个双重循环把Nstd从 0.05 扫到 0.4、NE从 20 扫到 200每个组合记录snr画热图峰值所在区域就是这组数据的合理参数区间。这套流程跑一遍要几十分钟但比反复手动改参数试错省心得多。验证改动有没有效果我的习惯是保留一份原始示例数据的基准结果固定的rng种子、固定的参数、固定的输入文件每次改完代码先跑基准imf矩阵和peakF与上一次完全一致才算没有破坏原有逻辑。EMD 这类迭代算法对数值扰动很敏感改一行边界插值就可能让整条 IMF 序列偏移没有基准对照根本发现不了。从那以后我每次拿到这类分解算法的包都先把NE降下来跑一遍基准、把结果存成.mat锁死再动任何一行代码都跟基准对一遍——这条规矩让我少返工很多次。希望帮到你。本文还有配套的精品资源点击获取