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 取 5～7、τ 取 1 是最常用配置，阈值 th 取 0.6 或 0.7。下面这张表列出频率混叠场景下参数的一般选法，其中 s 是原始信号的标准差。

需要提醒的是，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); % 存储第一次集总的中间结果 % 步骤1：EEMD 集总平均 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!) 把熵值归一化到 0～1 之间。参数说明：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 分解产生端点飞翼和模式混合。对这个测试信号调用 MEEMD：

k = 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，观察模态混叠消失的临界点；再调 N，N 从 100 提到 300，看 IMF 波形是否还抖动，抖动明显就继续加大；最后调 th，th 决定哪些分量被当作异常剔除。th 过小，正常周期分量因熵值偏高被误删；th 过大，异常冲击混入最终结果，排列熵检测形同虚设。打印每个 IMF 的排列熵值分布能帮助定阈值：

for i = 1:size(imfs, 2) fprintf('IMF%d PE=%.3f\n', i, permutation_entropy(imfs(:, i), 6, 1)); end

0.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 个样本点，却设置 m=7，此时重构向量个数太少，排列熵接近随机，阈值的区分能力几乎为零，代码能跑，但结果没有任何统计意义。

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_ext+1:end))]; [imfs_ext, res_ext] = meemd(x_ext, k, N, m, tau, th); imfs = imfs_ext(n_ext+1:end-n_ext, :); res = res_ext(n_ext+1:end-n_ext);

这段代码把开头和结尾各向外复制 10% 样本形成镜像，分解完成后把扩展部分直接裁掉。边界效应在 IMF 高频分量中最明显，裁掉后检查首尾两个周期是否平滑，如果仍有明显跳变，把 n_ext 提高到 20% 再试一次。

5.2 停止条件与 IMF 筛选次数

EMD 内部筛分过程有停止条件，经典实现用相邻筛分结果的标准差 SD 作判据：

[ SD = \frac{\sum_{t=0}^{T}|h_{i-1}(t) - h_i(t)|^2}{\sum_{t=0}^{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 组合的结果同时算出排列熵分布，横向对比后选出熵值分布最集中的一组参数，比单次试错更有判断力，也能省掉反复读图的精力。

本文还有配套的精品资源，点击获取