MEEMD程序详解:从EEMD到排列熵的MATLAB实现与参数调优 简介MEEMD改进集合经验模态分解与EEMD的MATLAB源码包面向信号处理、故障诊断等领域的科研人员与工程师用于解决非线性、非平稳信号的分解与特征提取问题适用于课程设计、论文复现和工程预研。资源采用RAR压缩共33个文件其中27个.m源码和6个.mat数据文件压缩包仅113KB涵盖meemd_ZHY.m、ceemd.m、emd.m、extrema.m等核心函数以及meemd_ZHY_jiedu.m注释版脚本s1.mat、s2.mat、ecg.mat等数据文件可直接作为测试信号。压缩包内按MD4_PEfenxi、MD5_MPEfenxi、MD6_MEEMDandPEfenxi等专题组织分别对应排列熵、多尺度排列熵、MEEMD与PE联合分析等实验便于读者按模块逐步深入。已有118人学习下载。通过源码可掌握噪声辅助分解、IMF平均、残余处理等关键步骤了解EEMD到MEEMD的改进逻辑并能将代码迁移至振动分析、电力系统故障诊断等实际场景是学习HHT方法及其改进算法的高性价比实用工具。1. MEEMD 程序到手后先搞清楚它在分解什么“附件2_MEEMD程序”这类压缩包在信号处理圈子里流通很广文件名里同时出现 MEEMD、EEMD 甚至手滑写错的“EEME”其实指向的是同一件事用 MATLAB 实现改进的集合经验模态分解。网上代码质量参差不齐直接运行经常报“未定义函数 emd”或者分解结果里模态混叠依旧严重因为 MEEMD 不是一个固定算法它在不同论文里代表不同的改进策略最常见的是把排列熵和 EEMD 结合用熵值识别异常 IMF再针对性处理。这篇内容把下载到的 MEEMD 程序拆开讲告诉你哪些代码是必须保留的哪些参数不调就白跑以及拿到源码后怎么在 MATLAB 环境里快速验证它分解得对不对。适合正在做振动分析、故障诊断、股票或气象时间序列分解的工程师。2. MEEMD 的算法基础从 EEMD 到排列熵改进2.1 EEMD 为什么需要“集合”经验模态分解EMD把信号按局部特征时间尺度分解成若干本征模态函数 IMF但模态混叠问题始终存在同一次分解中相近频率成分可能被分到不同 IMF一个 IMF 内部也可能混入多个频率。EEMD 的思路是在原始信号上叠加高斯白噪声利用白噪声在时域上均匀分布的统计特性让不同尺度的信号自动映射到合适的参考尺度上再重复多次取平均。每次加噪副本执行 EMD 后把对应序号 IMF 做集总平均白噪声在平均中互相抵消真实信号成分被保留。[ x_i(t) x(t) \epsilon w_i(t) ]这里的 \epsilon 是噪声幅值系数w_i 是零均值单位方差的高斯白噪声。EEMD 把“一次分解”变成“一群分解的均值”代价是计算量成倍增长同时引入两个新麻烦噪声幅值怎么选才既压制混叠又不污染信号少数 IMF 在特定时段可能残留较多噪声集总平均无法完全抹掉。MEEMD 的两个目标就是解决这两点它比 EEMD 多出来的改进通常表现在对异常 IMF 的识别和处理上。2.2 MEEMD 针对 EEMD 改进的三个位置排列熵的引入是 MEEMD 最常见的技术路径。我一般遇到的 MEEMD 源码主要在三个位置改进第一处在集总完成后计算每个 IMF 的排列熵将熵值高于阈值的 IMF 标为异常分量第二处不直接输出这些 IMF而是先对原信号剔除异常分量后的剩余信号继续做 EMD 或者再次集总第三处改进 IMF 的筛分停止条件避免过度筛分。需要说明的是不同版本 MEEMD 对“异常 IMF”的处理不一样部分程序会把排列熵高的异常 IMF 直接丢弃另一部分会将其叠加小噪声后再分解一次。拿到别人代码先看清楚它走的是哪条分支本文讨论的“附件2_MEEMD程序”走的是排列熵识别异常、从原信号中剔除后再分解一次的标准路径。对照来看EEMD 的问题在于“加噪后平均”对噪声残留的容忍度较高而 MEEMD 在 EEMD 的结果之上加了一道质检哪个 IMF 熵值过高就认为它还含有未抵消的随机成分或异常事件不让它直接进入最终结果。这个设计对冲击信号、突变信号处理效果提升明显但对纯周期信号排列熵识别反而可能把有效分量误删所以阈值选择是整个 MEEMD 代码里最重要的参数没有之一。2.3 排列熵用符号化方式测复杂度排列熵的优点是计算快、抗噪强、对数据长度要求不像样本熵那么苛刻。对长度为 N 的时间序列先重构为 m 维延迟向量[ X_i [x(i), x(i\tau), \dots, x(i(m-1)\tau)] ]对每个向量里的元素排序把排序后得到的序号组合作为排列模式统计每种模式出现频率 p_j归一化排列熵的计算公式为[ PE -\frac{\sum p_j \ln p_j}{\ln(m!)} ]周期信号的排列模式非常有限熵值低白噪声或随机冲击的排列模式高度复杂熵值接近 1。在 MEEMD 程序里m 取 57、τ 取 1 是最常用配置阈值 th 取 0.6 或 0.7。下面这张表列出频率混叠场景下参数的一般选法其中 s 是原始信号的标准差。参数常见范围选取依据噪声幅值系数 k0.1s 0.4s与信号标准差挂钩过小模态混叠抑制不住过大引入伪分量集总次数 N100 500越大集总平均越干净运行时间线性增长嵌入维数 m5 7样本点少于 200 时取 3防止模式种类不够延迟 τ1高采样率信号可适当增大到 23排列熵阈值 th0.6 0.8低于 0.5 会把正常 IMF 误判为异常高于 0.8 则漏判需要提醒的是m! 种排列模式需要足够多的重构向量去统计如果信号长度 N 明显小于 m! 的数值直方图会出现大量零频次归一化熵偏低排列熵失去区分能力。短数据场景下宁可把 m 调小也不要硬套论文里的推荐值。3. 读 MEEMD 的 MATLAB 源码函数结构与核心代码逐段拆解3.1 文件包里常见的文件组织方式这类源码文件夹一般有以下几类文件主函数 meemd.m辅助的 EEMD 或 EMD 函数排列熵函数 permutation_entropy.m有时还有边界延拓 mirror_extension.m 和画图脚本 test_demo.m。主函数负责信号输入、参数传递、调用集总循环EMD 部分通常基于 Rilling 的经典实现或作者自带的简化版排列熵函数单独成一个文件方便单独调试。直接双击附件里的 test_demo.m 跑不通时先看 MATLAB 当前路径有没有把子文件夹加入搜索路径再看主函数内部调用的函数名与实际文件名是否一致。手动在命令窗执行 which meemd 和 which emd 各试一次输出结果是空白或者“未找到”就说明路径或函数名不匹配。3.2 MEEMD 主函数的工作流下面这段代码是常见 MEEMD 主函数结构的浓缩版够在 R2016b 及之后版本运行注释里标明了每步在做什么。这里用占位函数 emd、pad2len 表示经典 EMD 内核和对齐工具实际运行时要替换成你自己源码里的对应实现。function [imfs, res] meemd(x, k, N, m, tau, th) % 输入x 信号列向量k 噪声幅值系数N 集总次数 % m/tau 排列熵维数与延迟th 异常判据阈值 x x(:); n length(x); raw zeros(n, N); % 存储第一次集总的中间结果 % 步骤1EEMD 集总平均 imf_layer {}; for i 1:N xn x k * std(x) * randn(n, 1); imfs_i emd(xn); % 每列一个 IMF最后一列为残差 raw raw pad2len(imfs_i, n); end raw raw / N; % 步骤2对每个 IMF 计算排列熵 imf_num size(raw, 2); pe zeros(imf_num, 1); for j 1:imf_num pe(j) permutation_entropy(raw(:, j), m, tau); end % 步骤3熵值高于阈值的 IMF 视为异常从原信号中剔除 valid pe th; reject ~valid; x_clean x - sum(raw(:, reject), 2); % 步骤4对净化后的信号再做一次 EMD imfs_final emd(x_clean); res imfs_final(:, end); imfs [raw(:, valid), imfs_final(:, 1:end-1)]; end逻辑说明步骤 1 完成经典的 EEMD所有加噪副本的 IMF 逐列对齐相加后取平均得到初步分解步骤 2 用排列熵逐个检查 IMF熵值高说明该分量还带有较强随机性步骤 3 把高熵分量从原始信号中减去实现“异常剔除”步骤 4 对净化后的信号再做一次普通 EMD最终输出由保留下来的低熵 IMF 和新分解出的 IMF 拼接而成。参数说明k 控制噪声强度std(x) 让噪声能量跟随信号尺度自动调整N 决定平均次数N 越大结果越稳定但耗时越长m 和 tau 直接影响排列熵的区分度th 决定“异常”的判定松紧度。这套流程里面最容易出问题的是 pad2len 这一步因为每次加噪后 EMD 分解出的 IMF 数量可能不一致常见做法是把缺少的行补零或者统一在最短 IMF 数量处截断。3.3 排列熵计算的实现与边界情况排列熵函数本身不长容易出错的是重复值的排序顺序。下面代码对相同值统一用其出现次序作区分避免 sort 的默认行为不稳定。function pe permutation_entropy(x, m, tau) N length(x); M N - (m - 1) * tau; patterns zeros(1, M); for i 1:M v x(i:tau:i (m - 1) * tau); [~, idx] sort(v, stable); patterns(i) polyval(idx, m); % 将序号转成唯一整数 end counts histcounts(patterns, 0:max(patterns)1); p counts(counts 0) / M; pe -sum(p .* log(p)) / log(factorial(m)); end逻辑说明每得到一个排序索引序列就用 polyval 把它编码成整数键histcounts 统计各模式出现次数除以重构向量总数 M 得到概率估计最后除以 log(m!) 把熵值归一化到 01 之间。参数说明tau 取 1 时等价于连续样本碰到采样率特别高、相邻点相关性过强的情况适当增大 tau 能让排列模式更丰富。这个函数还有一个隐含问题对直流分量敏感。输入 x 里若有明显均值大量重构向量的排序结果会完全相同PE 被严重低估在 MEEMD 流程里表现就是趋势项被误当成低熵有效分量保留下来。因此实际调用前最好先对信号做一次 detrend 或减去均值再做排列熵计算。4. 用 MEEMD 跑通一段仿真信号参数设置、输出校验与常见误用4.1 先造一个能评价分解质量的仿真信号真实数据不知道原始成分很难判断分解结果到底对不对。先用已知频率成分的信号验证算法本身再上真实数据。下面代码生成 1 秒、采样率 1000 Hz、包含 50 Hz 和 220 Hz 两个正弦分量以及白噪声的仿真信号fs 1000; t (0:999) / fs; x 0.8 * sin(2 * pi * 50 * t) 0.4 * sin(2 * pi * 220 * t); x x 0.15 * randn(size(t));参数说明两个频率 50 Hz 和 220 Hz 之间相差不到 2.5 倍频程普通 EMD 在这个频率比下容易出现模态混叠适合用来验证 MEEMD 的改进效果。噪声标准差取 0.15约为 50 Hz 分量幅值的五分之一信噪比不算低但已经足够让单一 EMD 分解产生端点飞翼和模式混合。对这个测试信号调用 MEEMDk 0.25; N 200; m 6; tau 1; th 0.7; [imfs, res] meemd(x., k, N, m, tau, th); figure; for i 1:size(imfs, 2) subplot(size(imfs, 2), 1, i); plot(t, imfs(:, i)); ylabel([IMF, num2str(i)]); end建议先把每个 IMF 的频谱画出来确认 50 Hz 和 220 Hz 分别落在哪个分量里。如果两个频率出现在同一个 IMF 中说明 k 太小或 N 不够如果某个 IMF 频率成分散乱优先调 k。4.2 参数调整的三条核心经验判断 MEEMD 结果好坏主要看三点分解出来的 IMF 是否对应原始频率是否存在一个 IMF 里同时出现两个主导频率原信号减去所有 IMF 加残差后的重建误差是否远小于信号本身。基于这三点参数调整有一个先后顺序先调 k把 k 从 0.1 慢慢升到 0.4观察模态混叠消失的临界点再调 NN 从 100 提到 300看 IMF 波形是否还抖动抖动明显就继续加大最后调 thth 决定哪些分量被当作异常剔除。th 过小正常周期分量因熵值偏高被误删th 过大异常冲击混入最终结果排列熵检测形同虚设。打印每个 IMF 的排列熵值分布能帮助定阈值for i 1:size(imfs, 2) fprintf(IMF%d PE%.3f\n, i, permutation_entropy(imfs(:, i), 6, 1)); end0.7 这个阈值不是万能值。如果打印结果显示所有 IMF 的 PE 都低于 0.5说明信号本身很干净th 可以适当降低到 0.6 以增强异常检测灵敏度如果大部分 IMF 的 PE 都在 0.8 以上说明噪声强度远高于预期应该先增大 N 而不是继续调 th否则整体都会被误删。4.3 重建误差与正交性校验分解是否可信可以用重建误差和正交性两个指标验证。重建误差衡量的是“分解再相加能否还原原信号”代码如下recon sum(imfs, 2) res; err recon - x; fprintf(最大重建误差%.3e\n, max(abs(err)));正常情况最大重建误差应该在 10^-14 到 10^-12 量级这个量级只受浮点精度影响。如果误差达到 10^-2 量级说明中间某一步把信号尺度改了最常见的错误是排列熵识别后把有效分量也从原信号里整体减掉却又没有在新分解结果中补回来。正交性指标按以下近似方式计算idx size(imfs, 2); orth zeros(idx, 1); for i 1:idx orth(i) sum(imfs(:, i) .* recon - imfs(:, i).^2) / sum(imfs(:, i).^2); end逻辑说明这个指标把每个 IMF 与重构信号的内积减去 IMF 自身能量再除以 IMF 自身能量衡量该 IMF 与其他分量的混叠程度。理想情况下两个不同频率的 IMF 正交指标接近 0指标绝对值大于 0.1 就要怀疑存在模态混叠或虚假分量。将 EEMD 和 MEEMD 的分解结果放到同一张图上对比功率谱可以直观看出 MEEMD 的改进是否有效这一步在故障诊断场景里几乎是必须的。4.4 误用把 MEEMD 当滤波器用最常见的误用是以为 MEEMD 可以替代带通滤波器直接拿分解结果中感兴趣的 IMF 做后续分析而不检查其物理意义。MEEMD 分出来的分量不一定对应某个确定的物理振源特别是噪声较强时一个 IMF 可能由多个频率成分拼凑出来只是排列熵恰好不高而已。正确做法是先做频谱分析确认每个 IMF 的主频率与实际工况一致再决定保留还是剔除。另一个常见误用是短数据配合高嵌入维数比如 0.2 秒、采样率 500 Hz、一共 100 个样本点却设置 m7此时重构向量个数太少排列熵接近随机阈值的区分能力几乎为零代码能跑但结果没有任何统计意义。5. 最后一道关边界效应、停止条件与 MEEMD 结果验证技巧5.1 边界效应先解决EMD 类算法在信号两端天然容易发散端点处极值点缺失包络线向外延伸时产生大幅振荡。常见做法是分解前做镜像延拓把信号两端反射扩展若干个极值周期再分解分解后截掉延拓部分。源码里如果没有这步运行结果的前几个点和后几个点往往异常大解决办法是在调用 meemd 前手动扩展和截断n_ext round(0.1 * length(x)); x_ext [flipud(x(1:n_ext)); x; flipud(x(end-n_ext1:end))]; [imfs_ext, res_ext] meemd(x_ext, k, N, m, tau, th); imfs imfs_ext(n_ext1:end-n_ext, :); res res_ext(n_ext1:end-n_ext);这段代码把开头和结尾各向外复制 10% 样本形成镜像分解完成后把扩展部分直接裁掉。边界效应在 IMF 高频分量中最明显裁掉后检查首尾两个周期是否平滑如果仍有明显跳变把 n_ext 提高到 20% 再试一次。5.2 停止条件与 IMF 筛选次数EMD 内部筛分过程有停止条件经典实现用相邻筛分结果的标准差 SD 作判据[ SD \frac{\sum_{t0}^{T}|h_{i-1}(t) - h_i(t)|^2}{\sum_{t0}^{T} h_{i-1}(t)^2} ]当 SD 低于设定阈值常见 0.2 或 0.25时筛分停止。代码中如果找不到这个阈值搜索 sd、sift 或 stop 相关变量。阈值过大导致欠筛分IMF 不满足局部对称性阈值过小导致过度筛分把幅值调制的信号变成频率调制产生没有物理意义的振荡。MEEMD 对筛分次数的敏感度比 EEMD 低一些集总平均会中和部分误差但完全不管它同样会得到过度平滑的波形处理非平稳信号时尤其明显。5.3 用希尔伯特包络验证 MEEMD 是否保留有效成分验证 MEEMD 分解质量最直观的进阶手段是检查单分量 IMF 的希尔伯特包络是否是缓变曲线。如果某个 IMF 是纯净的单频分量其解析信号的幅值包络应该是缓慢变化的如果包络线持续快速跳动说明该 IMF 内部还有残余噪声或混叠分量。验证代码如下h hilbert(imfs(:, 1)); amp abs(h); env smooth(amp, 20); figure; plot(t, imfs(:, 1)); hold on; plot(t, env, r-, LineWidth, 1.5);不要只盯着第一个 IMF。MEEMD 输出的第一个 IMF 往往包含最高频噪声成分去噪应用里保留 IMF1 有时反而会引入高频毛刺。建议把每个 IMF 的包络都画一遍包络跳变明显的分量优先回炉重调。最后补一个批处理技巧把多个 k 和 th 组合的结果同时算出排列熵分布横向对比后选出熵值分布最集中的一组参数比单次试错更有判断力也能省掉反复读图的精力。本文还有配套的精品资源点击获取