简介面向非线性非平稳信号分析需求这套MATLAB程序实现了集合经验模态分解EEMD通过向原始信号叠加均匀白噪声进行多次分解并取平均有效抑制EMD模态混叠问题适用于振动信号、生物医学信号等场景的时频特征提取。与经典EMD相比EEMD以噪声辅助分析改善分解稳定性是许多科研与工程项目的常用起点。压缩包共46个文件以42个m脚本为核心另含csv数据文件、txt说明文档及p文件整体仅189KB。代码覆盖EEMD主流程、IMF筛选、Hilbert谱与边际谱分析、显著性检验、绘图与后处理等环节模块划分清晰既便于初学者对照运行理解算法原理也方便研究者修改参数进行对比实验。目前已有5188人学习下载小巧体积与完整功能使其成为信号处理课程设计、科研预研的实用工具。1. EEMD集合经验模态分解的MATLAB程序包一套到手即用的非平稳信号分解工具EEMD集合经验模态分解的MATLAB程序代码是我处理振动信号时用得最频繁的一套工具。最早我用纯EMD分解轴承故障信号第一阶IMF里同时混着50Hz工频和冲击成分第三阶又把主频拆成两段典型的模态混叠。换成EEMD之后白噪声把不同尺度的成分撑开混叠基本消失。这套代码包封装了完整的集合经验模态分解流程输入一段一维信号输出从高频到低频的IMF分量和残差可直接用于故障诊断、信号去噪、时频分析和Hilbert谱绘制。适合做机械故障诊断和非平稳信号分析的工程师、研究生也适合想快速跑通一套可用代码的新手。2. EEMD原理与参数选型为什么加噪声反而更准NE和Nstd怎么定2.1 从EMD到EEMD模态混叠怎么破先回顾EMD的基本流程。经验模态分解的核心假设是任何复杂信号都能拆成若干个固有模态函数IMF加一个残差。IMF要满足两个条件极值点个数与过零点个数相等或最多相差一个上下包络关于时间轴局部对称。算法上就是反复用三次样条构造上下包络、取包络均值、从原信号里减掉均值筛分到满足停止条件后得到第一阶IMF再对残差重复这个过程直到残差近似单调。这里的停止条件是一个经常被忽略的细节。早期实现用相邻两次筛分结果的归一化标准差判断阈值常取0.2到0.3后来的实现改用包络均值振幅小于局部幅值一定百分比来停止。这套代码包内部具体用哪种准则跑出来的IMF数量会略有差别但不影响整体使用。如果你发现分解出的IMF数量异常多比如几十阶大概率是停止阈值设得太严筛分过度迭代了。EMD最让人头疼的问题叫模态混叠表现是同一尺度的成分被拆到多个IMF里或者一个IMF里混进差异很大的时间尺度。根源在于EMD依赖信号局部极值的分布一旦信号里有间歇性成分、脉冲冲击或噪声极值点分布不均匀包络就被带偏筛分结果自然乱。核心筛分循环长这样各家实现大同小异% 单次EMD筛分的核心循环示意 % r: 当前残差信号列向量 h r; SD 1; while SD 0.3 % 筛分停止阈值0.2~0.3常见 [upper, lower] envelope_win(h); % 三次样条构造上下包络 m 0.5 * (upper lower); % 包络均值 h_new h - m; % 减去均值向IMF靠近 SD sum((h_new - h).^2) / sum(h.^2); % 相邻两轮差异的归一化指标 h h_new; end imf_k h; % 这轮筛分得到的第 k 阶 IMF这段代码里envelope_win内部完成三次样条插值SD是无量纲的收敛指标两轮筛分结果几乎不变时循环退出。用while SD 0.3而不是固定迭代次数是为了让不同尺度的成分都能按需收敛如果改成固定迭代次数低频成分往往还没筛透循环就结束了分解出的低频IMF会带着明显的高频残留。EEMD解决模态混叠的思路很直接既然极值分布不均匀是根源就人为把极值分布抹匀。做法是在原始信号上叠加白噪声利用白噪声在频带内均匀分布的特性迫使信号在各尺度上的极值点被填充到近似均匀对加噪信号做一次标准EMD重复若干次再把对应阶次的IMF做集合平均。白噪声均值为零集合次数足够多时噪声相互抵消保留下来的就是真实的结构成分。这就是“集合经验模态分解”里“集合”两个字的来由。2.2 核心参数NE集合次数与Nstd噪声幅值EEMD能压住模态混叠背后有一个定量关系多次平均后残留在IMF里的噪声幅值大致满足 err Nstd / sqrt(NE)。Nstd是加入白噪声的标准差与原始信号标准差之比NE是集合次数。这个关系是整套参数选取的根基几乎所有调参建议都能从它推出来。参数含义常用范围设太小的后果设太大的后果Nstd噪声标准差 / 信号标准差0.05~0.3通常0.1~0.2模态混叠抑制不住伪IMF增多有效成分被淹没NE集合平均次数100~500噪声抵消不干净结果抖动计算时间线性增长从 err Nstd / sqrt(NE) 能读出三层实操含义。第一Nstd越大单次加噪抹匀极值分布的效果越强但残留噪声也越大需要在“抑制混叠”和“残留少”之间折中。第二NE越大残留噪声越小但计算量线性增长500次集合平均对长信号来说相当耗时。第三Nstd从0.1调到0.2要保持同样的残留噪声水平NE需要变成原来的4倍因为sqrt(NE)是平方根关系噪声幅值翻倍抵消次数就要翻四倍。我一般这样起步先用Nstd0.2、NE100跑一遍看IMF数量和波形如果某些IMF明显抖动把NE加到200~300如果模态混叠仍然明显本该单一频率的IMF出现大幅调幅包络再把Nstd加大到0.2~0.3。对仿真的高信噪比信号Nstd取0.05~0.1就够对实测振动信号这种信噪比低、间歇性冲击多的场景Nstd取0.1~0.2更稳。还有一个判断NE够不够的经验方法分别用NE100和NE300跑两遍把对应阶次的IMF相减看差异幅值。如果差异的均方根小于该阶IMF幅值均方根的1%说明NE已经够如果差异明显继续加NE。不要一上来就设1000先花一分钟做这个对比实验能省下几十分钟的无效计算。提示Nstd是相对值不是绝对值。它表示白噪声标准差与信号标准差的比值信号幅值不影响这个参数的量级不同量纲的信号可以直接用同一套推荐值。2.3 代码包文件构成与主函数接口常见的MATLAB EEMD实现文件组织大致是一个主入口函数通常叫eemd.m内部调用单次EMD筛分函数和包络计算函数外加一个示例脚本演示调用。这套代码包的典型结构如下文件作用eemd.m主函数完成加噪、集合平均全流程emd_once.m单次EMD筛分包含停止准则和包络插值envelope_win.m用三次样条构造上下包络example.m示例脚本生成仿真信号并演示调用主函数的接口在实现中大同小异最通用的一种是这样% 主函数入口通用接口 % modes eemd(y, Nstd, NE) % y : 输入信号N x 1 列向量 % Nstd : 白噪声幅值系数相对 y 的标准差 % NE : 集合次数 % modes: N x M 矩阵前 M-1 列是 IMF最后一列是残差这个接口有两点要注意。第一输入必须是列向量如果从文本文件或Excel读进来是行向量先转置再调用否则内部矩阵切分时很容易维度报错。第二残差和IMF是放在同一个矩阵里的用的时候要么对前M-1列单独操作要么把最后一列单独取出别直接拿整个矩阵去做相关分析或频谱分析否则残差的低频趋势会把结果带偏。另外输入信号在进EEMD之前最好做两步预处理一是去均值让信号围绕零轴波动否则残差里会多出一条直流分量二是检查是否有明显趋势项有的话先做一次多项式去趋势。很多初学者跳过预处理结果分解出的第一阶IMF和残差都带着奇怪的缓变形状还以为是代码问题其实是信号没收拾干净。3. 在MATLAB里跑通EEMD调用、参数设置与结果解读3.1 最小调用示例输入输出到底长什么样先搭一个最小可复现的调用环境。假设有一段1秒、采样率1000Hz的仿真信号由30Hz和80Hz两个正弦叠加再混入一点随机噪声。这是验证EEMD最干净的实验台因为真实成分已知分解完可以直接对答案。% 仿真信号30Hz 80Hz 正弦叠加 fs 1000; % 采样率 1000 Hz t (0:1/fs:1-1/fs); % 1 秒时长列向量 x sin(2*pi*30*t) 0.5*sin(2*pi*80*t) 0.05*randn(size(t)); % 调用EEMD Nstd 0.2; % 白噪声幅值系数 NE 200; % 集合次数 [modes, its] eemd(x, Nstd, NE); % 查看输出规模 fprintf(输入信号长度: %d\n, length(x)); fprintf(modes矩阵: %d x %d\n, size(modes,1), size(modes,2));这段代码里有两个细节值得说。第一0.05*randn(size(t))加的噪声很小目的是贴近实测信号里那种微小扰动——如果一点噪声都不加分解太理想反而掩盖参数的影响。第二[modes, its]双输出写法不是所有版本都有its是每阶IMF完成筛分用的迭代次数调试时才看如果报输出参数太多改成modes eemd(x, Nstd, NE)即可。跑完打印出的modes矩阵行数等于信号长度列数等于IMF数量加1。这段仿真信号正常会分解出3~4阶IMF加1列残差。有时候同样代码跑两次IMF列数会差一列——不是bug是白噪声随机性导致的第4章专门讲怎么处理。its输出也值得瞄一眼。通常前几阶高频IMF的筛分次数在几十次量级越到后面的低频IMF需要的筛分次数越少。如果its出现异常大或异常小的值多半是停止准则和信号尺度不匹配这时候优先排查代码内部的包络插值逻辑比调Nstd更有效。3.2 完整示例把仿真信号分解成IMF并画图拿到modes矩阵后第一步永远是画图。很多人上来就做Hilbert谱或算能量占比分解本身是乱的后面全是白做。先把每阶IMF按时间轴画出来基本形态对不对一眼就能看出来。% 把每个IMF和残差画出来 figure(Color,w,Position,[100 100 560 800]); n_modes size(modes,2); for k 1:n_modes subplot(n_modes,1,k); plot(t, modes(:,k), LineWidth, 0.8); grid on; if k n_modes title([残差 R, num2str(k-1)]); else title([IMF , num2str(k)]); end ylim([-1.5 1.5]); end xlabel(时间 (s));画完重点看三件事。第一IMF1应当基本是80Hz的高频分量、IMF2基本是30Hz的低频分量两者是近似等幅的正弦如果出现明显调幅包络说明模态混叠没压住Nstd要加大。第二各阶IMF在同一时刻不应同时出现大振幅且频率接近的成分否则就是分解边界不清。第三残差应当是缓慢变化的趋势线对这段信号残差振幅应接近0如果残差还有明显的正弦形状说明筛分没筛干净需要检查停止阈值或提前做去趋势。如果第一阶IMF上带着低幅值随机毛刺属于正常现象。EEMD的集合平均只能压制白噪声残留不能完全消除高频IMF带点毛刺是物理上必然的。不要试图把Nstd设成0来“让结果更干净”那会直接退回EMD的老路模态混叠马上回来。还可以把原信号和全分量重构结果叠加对比验证分解是否完备% 全分量重构所有IMF加残差 x_sum sum(modes, 2); figure(Color,w); plot(t, x, b, LineWidth, 1.2); hold on; plot(t, x_sum, r--, LineWidth, 1.0); legend(原始信号, IMF残差重构); xlabel(时间 (s)); ylabel(幅值); title(分解完备性验证);两条曲线应当几乎重合说明分解没有丢能量。如果明显分叉优先怀疑代码内部有没有做端点裁切或降采样这类隐藏操作会在边界处丢数据。3.3 后处理Hilbert谱与边际谱的可视化EEMD产出IMF后最常见的下游操作是时频分析。每个IMF是窄带信号对它做Hilbert变换求瞬时频率和瞬时幅值拼起来就是Hilbert时频谱。MATLAB从2018a起在信号处理工具箱里直接提供hht函数不用自己写Hilbert变换非常省事。% 对IMF做Hilbert-Huang变换 % 语法: [hs, f, t] hht(imf, fs) imf modes(:, 1:end-1); % 去掉残差 [H, f_axis, t_axis] hht(imf, fs, FrequencyLimits, [0 300]); % Hilbert谱 figure(Color,w); imagesc(t_axis, f_axis, H); axis xy; colormap(parula); colorbar; xlabel(时间 (s)); ylabel(频率 (Hz)); title(Hilbert-Huang 时频谱); % 边际谱对时间维求和 marginal sum(H, 2); figure(Color,w); plot(f_axis, marginal, LineWidth, 1.2); xlabel(频率 (Hz)); ylabel(幅值); title(边际谱);两个参数要解释。一个是FrequencyLimits, [0 300]采样率1000Hz时Nyquist上限是500Hz但这段信号只有30Hz和80Hz把频率轴限制到300Hz以内能拉开展示如果信号有高频冲击成分上限按需放宽。另一个是imf modes(:, 1:end-1)残差不满足IMF窄带条件直接丢给hht会产生大量无意义的瞬时频率跳动必须去掉。边际谱上应出现两个明显峰一个在30Hz附近、一个在80Hz附近。峰值位置和原始频率对得上说明EEMD分解可信如果峰偏了或冒出一堆杂峰回查Nstd和NE大概率参数没收敛。如果信号本身包含频率随时间快速变化或突变成分Hilbert谱上会出现从高频扫到低频的曲线或竖直亮线这是时频定位的正常表现不要误判成干扰。另有一点要记住hht输出的瞬时频率是逐点计算的对噪声敏感谱图里零星的高能量点大概率是噪声残留把频率轴聚合到几个主峰再看结论会稳定很多。4. EEMD使用中的常见问题与避坑指南4.1 同样的代码跑两次IMF结果不一样现象输入信号、参数完全一致分两次运行IMF波形和列数有细微差异有人以为程序有随机性错误。原因EEMD每次迭代都生成新白噪声序列MATLAB的randn默认随机种子每次启动不同噪声序列不同集合平均结果自然有微小差异。这是算法的统计特性不是bug。解决调用前固定随机种子让噪声序列可复现rng(0); % 固定随机种子保证结果可复现 modes eemd(x, 0.2, 200);论文图表或项目报告用的结果一定要固定种子否则重跑一遍图就对不上非常被动。算法调试阶段不固定也没关系差异通常在1%以内。4.2 信号两端出现大幅摆动端点效应现象IMF在首尾两端突然大幅振荡尤其第一阶IMF在起点和终点有明显包络尖峰中间段反而正常。原因三次样条包络依赖信号极值点端点处没有数据支撑样条在端点的斜率外推不受约束包络在两端发散筛分过程把发散带进了IMF。解决分解前做端点延拓常见做法是镜像延拓两头各延长一段镜像数据分解后剪掉% 镜像延拓前后各延长300个点 ext_len 300; x_ext [flipud(x(2:ext_len1)); x; flipud(x(end-ext_len:end-1))]; % 对延拓信号做EEMD modes_ext eemd(x_ext, 0.2, 200); % 剪掉延拓部分只保留原始区间 modes_cut modes_ext(ext_len1 : end-ext_len, :);延拓长度取信号长度的10%左右比较稳妥太长增加计算量太短压不住端点发散。EEMD集合平均对端点效应本身有一定平滑作用但如果信号只有几百个点端点效应占比很大镜像延拓几乎是必须的。4.3 Nstd设太大分解出一堆伪IMF现象Nstd设成0.5甚至1分解结果多出好几阶IMF每阶波形杂乱无规律边际谱上大量没有物理意义的峰。原因噪声幅值过大时单次筛分过程被噪声主导产生很多完全由噪声引导的极值点平均后保留下来就成了伪分量。解决把Nstd降回0.05~0.2区间同时提高NE补偿残留噪声。判断伪IMF的常用手段是计算每阶IMF与原始信号的相关系数% 各IMF与原始信号的相关系数 corr_vals zeros(size(modes,2), 1); for k 1:size(modes,2) cc corrcoef(x, modes(:,k)); corr_vals(k) cc(1,2); end % 经验阈值最大相关系数的十分之一 threshold max(corr_vals) / 10; valid_idx find(corr_vals threshold); fprintf(有效IMF索引: %s\n, mat2str(valid_idx));相关系数高的IMF与原始信号结构强相关大概率是真实成分很低的基本都是噪声残留。阈值是经验值不同场景要微调但作为快速筛选很实用。4.4 NE设太大计算时间长到怀疑人生现象NE设1000、信号几万点跑几十分钟没结果电脑风扇狂转。原因EEMD要对加噪信号重复做完整EMDNE次串行分解每次都涉及包络构造和筛分迭代计算量随NE线性增长。解决先算账。NE从100加到400残留噪声只减一半计算时间翻4倍。我的做法是先用NE50快速跑通流程看IMF形态确认参数方向和预处理没问题再提高到200~500做正式结果。如果装了MATLAB并行计算工具箱可以把NE次循环改成parfor并行原理上是把集合平均的每个独立EMD分配到不同worker% 并行加速集合平均示意具体实现随代码版本而异 % 初始化累加器后 % parfor i 1:NE % noise_std Nstd * std(y); % y_noisy y noise_std * randn(size(y)); % modes_i emd_once(y_noisy); % 单次EMD结果 % acc acc modes_i; % 累加最后统一除以NE % end % modes acc / NE;注意parfor里累加、最后统一除以NE的顺序不能乱否则和串行版结果对不上。另外对超长信号先用降采样验证算法参数再回到原始采样率跑正式结果是省时间的血泪经验。4.5 输入是行向量程序直接报错现象调用eemd时报维度不匹配或下标超界检查代码逻辑发现没问题。原因很多数据读取方式比如CSV导入得到1xN行向量主函数内部按Nx1列向量矩阵化尺寸对不上。解决调用前强制转列向量x x(:); % 无论原来什么形状统一转成列向量用x(:)比用x稳因为x已经是列向量时x(:)无副作用而x会把行向量和列向量互换反而容易埋雷。5. 进阶用相关系数筛选IMF并重构去噪信号EEMD的最终目的不只是画几张IMF图而是把分解结果用于实际任务最常见的是去噪和故障特征提取。去噪思路很直接把噪声主导的高频IMF剔除用有意义IMF重构信号。我一般分三步。第一步计算每阶IMF与原始信号的相关系数按4.3的方法筛出有效IMF第二步把不满足条件的IMF尤其是第1阶纯噪声IMF置零第三步把剩余IMF和残差相加得到重构信号。好处是既保留信号主要能量又滤掉随机噪声比全量重构更有工程意义。% 筛选有效IMF并重构 corr_vals zeros(size(modes,2), 1); for k 1:size(modes,2) cc corrcoef(x, modes(:,k)); corr_vals(k) cc(1,2); end threshold max(corr_vals) / 10; valid_idx find(corr_vals threshold); % 重构只累加有效IMF x_recon sum(modes(:, valid_idx), 2); fprintf(参与重构的IMF编号: %s\n, mat2str(valid_idx)); % 重构误差 err x - x_recon; fprintf(重构相对误差: %.4f%%\n, 100 * norm(err) / norm(x));重构误差反过来能验证参数选得对不对。误差在5%以内说明有效IMF抓住了主要成分超过20%多半是阈值定太高丢掉了真实成分或Nstd偏大破坏了有效IMF。我有个习惯每次调完参数都重构一次对答案误差突然变大就意味着前面某个环节出了岔子。这套筛选重构的思路也可以直接用于故障特征提取把重构后的信号做包络解调或频谱分析故障特征频率比在原始信号上清晰得多因为EEMD替你滤掉了一部分宽带噪声。从那以后我每次用EEMD都强制走一遍固定流程先固定随机种子再镜像延拓然后Nstd0.2、NE100快速试跑画IMF图确认形态最后相关系数筛选加重构误差验证全部通过才进入正式分析。这套习惯帮我省掉了无数个调参的加班夜。希望帮到你。本文还有配套的精品资源点击获取