先说个背景。经验模态分解（EMD）这类自适应方法，在单通道非平稳信号分析里已经不算新鲜了，它能按频率从高到低把信号拆成若干个固有模态函数（IMF）。但真实工程信号很少有单通道的，脑电是几十个导联，振动监测是多个测点，气象数据也是多站点同时采样。逐通道单做EMD有个很头疼的问题：每个通道分解出来的IMF个数不一样、中心频率对不上、相位也七零八落，后续想做通道间的相关性或者统计比较，根本无从下手。多元经验模态分解（MEMD）就是为了解决这个问题出现的。它把多通道数据当作一个整体，在多元空间里联合筛选IMF，保证各通道分解出的模态天然对齐。这篇文章就用Python把MEMD的核心逻辑拆开、手写一遍，顺便给出能直接跑起来的完整代码。

这套算法适合谁看？如果你在用多通道数据做时频分析、想避免逐通道EMD的模态错位问题，或者单纯好奇“多元分解”内部到底怎么运作，那这篇应该对你有用。我会尽量把原理讲得通俗，把代码控制在核心可用的程度，同时把那些文档里不会写的细节坑点也一并交代。

1. MEMD是什么：单通道EMD的局限与多元分解思路

1.1 从单通道EMD说起

先回顾一下单通道EMD在做什么。对一条信号x(t)，EMD的筛选过程是：找到信号所有局部极大值和局部极小值，分别用样条插值得到上包络和下包络，取两者均值m(t)，然后用x(t)减去均值得到候选分量d(t)。这个“减去均值包络”的操作会反复迭代，直到d(t)满足IMF的两个条件：一是极值点数量与过零点数量相差不超过1，二是上下包络均值在局部上接近零。分离出一个IMF后，从原信号里减去它，对剩余信号继续重复，直到残余分量变成单调趋势或幅度小到可以忽略。

这套逻辑本身是自适应的，不需要先验的基函数，对非线性、非平稳信号尤其友好。但问题也很明显：它是单通道算法。多通道数据如果每个通道独立做一遍，得到的IMF集合完全没法保证一一对应。刚做过这类工作的朋友应该都见过这种场面——通道1筛出了5个IMF，通道2只筛出4个，明明同一个3Hz的振荡成分，通道1落在IMF2，通道2却落在IMF3，相位和幅值也都对不上。想比较通道间的关系，要么强行对齐，要么手动挑选成分，既费劲又容易出错。

1.2 什么叫“联合分解”

MEMD的核心改动是：不再把每个通道当成孤立信号处理，而是把全部通道拼成一个多元信号X(t)=[x1(t), x2(t), ..., xn(t)]，在n维空间里做联合筛选。每次迭代时，算法会在n维单位超球面上取很多方向向量，把多元信号向这些方向投影，再根据投影的极值信息构造一个“多元包络均值”，所有通道共享同一个均值包络，一起作差更新。这样一来，每个IMF都是通道间对齐的——同一阶IMF在时间维上天然同步，频率成分也趋于一致，便于做跨通道的时频分析、相干分析和统计建模。

这里可以做个简单类比。单通道EMD就像让一堆人各自独立走路，步调和路线全凭个人发挥，最后你想让他们排成一排几乎不可能；MEMD等于给他们栓了一根共同的绳子，每步都朝同一个方向、按同一个节拍走，虽然每个人姿态还有差异，但整体步调是对齐的。后面要讲的多变量包络均值，就是这根绳子的工程实现。

需要明确，MEMD不是简单的“多通道信号拼一起做EMD”。它真正的技术点有三个：方向向量的采样、多元包络均值的计算、筛选停止准则的设计。下面一节就来拆这三个点。

2. 算法核心拆解：方向向量、多元包络均值与筛选停止

2.1 方向向量：为什么必须均匀采样

MEMD第一步是在n维单位超球面上生成K个方向向量v_k，然后计算投影：

p_k(t) = v_k^T X(t)

这个投影把多元信号压成一维标量。之所以要投影，是因为高维空间里没法直接定义“极值点和包络”，只能先沿特定方向看信号形态，再用这些方向上的信息反推多元包络。

方向向量的质量直接影响分解结果。如果方向分布不均匀，很多方向密集在一起，某些方向的信号形态会被过度强调，另一些方向则完全没被观测到，包络均值就会有偏。实际中不建议直接用随机均匀分布，因为随机采样在有限样本量下很容易在球面上聚团，每跑一次结果还不同。更稳妥的办法是用低差异序列在球面上做确定性采样，其中最常用的是Hammersley序列。

Hammersley序列本质上是多维的拟随机序列，它能把采样点均匀地填充到单位超立方体里。以二维为例，第一维取i/N，第二维通过以2为底的radical inverse生成。把它推广到多维：第一维仍然是i/N，后面每一维用不同的素数做radical inverse。得到超立方体上的点后，再利用球坐标变换映射到超球面：

第一个角度取φ0 = 2π·u0，作为经度方向的角；其余角度取φi = arccos(2·ui - 1)，对应纬度方向的角。转换到笛卡尔坐标后，得到的单位向量就近似均匀分布在n维球面上。方向数K到底取多少合适，我在后面调参部分会详细说，先记住经验值：多通道场景至少要64个方向。

2.2 多元包络均值：MEMD最关键的一个细节

方向向量有了，下一步是计算多元包络。这是和单通道EMD差异最大的地方，也是最容易理解错的地方。

对每个方向v_k，先算投影p_k(t)，然后找p_k(t)的全部局部极大值点，记录这些点的时间位置。搞清楚一件事：找的是局部极大值，不是把极大和极小合在一起。为什么只找极大值？因为方向向量遍布整个球面，既包含正方向也包含负方向。在正方向上的投影极大值对应信号在某方向上的“上包络触碰点”，在负方向上投影的极大值对应原方向上的“下包络触碰点”。正负方向合在一起，上下包络的信息就都包含进来了。

找到极值位置后，在这些时间点上对原始多元信号的每个通道分别做三次样条插值，得到该方向下的多元包络曲线e_k(t)。对K个方向重复这个过程，把所有方向得到的包络逐个通道取平均：

m(t) = (1/K) Σ_k e_k(t)

m(t)是一个n维向量时间序列，就是这次筛选要用的均值包络。整个流程用伪代码描述就是：

1. 生成K个方向向量v_k；
2. 计算每个方向的投影p_k(t) = v_k^T X(t)；
3. 找出p_k(t)的局部极大值点时间索引；
4. 在这些索引位置，对X(t)的每个通道分别样条插值，得到e_k(t)；
5. 所有方向得到的e_k(t)取平均，得到m(t)；
6. X(t)减去m(t)，判断是否满足多元IMF停止条件。

这里有个容易出现原则性错误的地方：插值的对象是原始多元信号X(t)的各通道，而不是投影p_k(t)。如果你对投影本身插值再投影回去，得到的包络会偏离真实多元包络。我在初版代码里就踩过这个坑，后面会用仿真实验说明区别。

还有一个细节值得注意。经典EMD里上包络和下包络分开插值、再取平均；MEMD则是按方向生成包络、再跨方向平均。当通道数为1时，方向向量只有+1和-1两个，正方向包络穿过局部极大值点，负方向包络穿过局部极小值点，两者平均恰好退化成经典EMD的上下包络均值逻辑。这说明MEMD在单通道场景下是可以自然兼容EMD思路的。

2.3 筛选停止准则：包络均值的“足够小”怎么定义

单通道EMD判断筛选是否结束，看的是IMF的两个条件。多元场景下，结束条件变成：所有通道上，均值包络m(t)在任意时刻都“足够接近零”。数值实现上当然不可能要求严格为零，一般用能量比值作为停止阈值：

criterion = Σ_t ‖m(t)‖² / Σ_t ‖X(t)‖²

分子是当前均值包络的能量，分母是待筛选信号的能量。当这个比值小于某个阈值（比如0.001或0.01）时，认为包络均值已经小到可以忽略，筛选结束，当前的d(t)就是一个多元IMF。

这个准则直观理解很舒服：多次迭代后均值包络逐渐逼近零，说明信号和它的对称中心越来越重合，就可以收手了。阈值设太小会造成过筛分，把噪声一点点扣出来；设太大又会提前停止，模态分解不彻底。我在实际使用中一般先试0.01，如果发现模态混叠明显再往小调。

另外要限制最大筛选次数。因为高维包络插值在边界处可能不稳定，偶尔会出现包络均值不小、但信号已经快被磨平的尴尬情况，没有最大迭代次数保护就会陷入死循环。经验上单次筛选不要超过100轮，多数信号在20到60轮之间就能收敛。

当残余信号在所有通道上的极值点总数太少，连包络插值都做不出来时，整个分解过程就该停了。这个停止条件和单通道EMD“残余信号单调即停止”是类似的。

3. 从零写出可运行的Python实现

下面给出我整理过的精简版实现。为了便于阅读，我把函数拆成了几个模块：方向向量生成、极值点查找、包络均值、筛选主循环。完整代码可以在本地直接复制运行，依赖只需要numpy和scipy两个库。

3.1 工具函数：Hammersley采样与球面方向向量

import numpy as np from scipy.interpolate import CubicSpline PRIMES = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71] def _radical_inverse(idx, base): """以base为底数的radical inverse，把整数映射到[0,1)""" value = 0.0 f = 1.0 / base i = idx while i > 0: value += f * (i % base) i = i // base f /= base return value def hammersley_points(n_points, dim): """在[0,1]^dim内生成n_points个Hammersley低差异序列点""" seq = np.empty((n_points, dim)) for i in range(n_points): seq[i, 0] = (i + 0.5) / n_points for d in range(1, dim): seq[i, d] = _radical_inverse(i, PRIMES[d - 1]) return seq def unit_vectors_on_hypersphere(n_vecs, dim): """在dim维单位超球面上生成n_vecs个方向向量""" if dim == 1: return np.array([[1.0], [-1.0]]) u = hammersley_points(n_vecs, dim - 1) angles = np.empty((n_vecs, dim - 1)) angles[:, 0] = 2 * np.pi * u[:, 0] if dim > 2: angles[:, 1:] = np.arccos(2 * u[:, 1:] - 1) vecs = np.ones((n_vecs, dim)) for i in range(dim - 1): p = np.ones(n_vecs) for j in range(i): p *= np.sin(angles[:, j]) p *= np.cos(angles[:, i]) vecs[:, i] = p p = np.ones(n_vecs) for j in range(dim - 1): p *= np.sin(angles[:, j]) vecs[:, -1] = p return vecs

这里Hammersley序列的作用是把方向向量尽量均匀铺在球面上。如果没有这个工具，直接用np.random.randn再加归一化也能得到方向向量，但效果和稳定性都不如低差异序列。你可以在实际项目里对比一下，随机采样的方向向量在不同运行之间会造成微小但可见的分解差异，而Hammersley序列每次跑出来的结果是完全一致的。

3.2 极值点查找与多元包络插值

def _local_maxima(x): """返回一维信号的局部极大值索引（不处理平顶，教学简化版）""" if len(x) < 3: return np.array([], dtype=int) dx = np.diff(x) prev = dx[:-1] curr = dx[1:] maxima = np.where((prev > 0) & (curr < 0))[0] + 1 return maxima def _mean_envelope(X, directions): """ 计算多元信号的均值包络 X: shape = (n_channels, n_samples) directions: shape = (n_dirs, n_channels) """ n_ch, n_samples = X.shape t = np.arange(n_samples) env_sum = np.zeros_like(X) valid_dirs = 0 for k in range(directions.shape[0]): proj = directions[k] @ X # 方向向量与多元信号的内积 maxima = _local_maxima(proj) if len(maxima) < 3: continue env_k = np.empty((n_ch, n_samples)) for c in range(n_ch): cs = CubicSpline(t[maxima], X[c, maxima], bc_type='natural') env_k[c] = cs(t) env_sum += env_k valid_dirs += 1 if valid_dirs == 0: return np.zeros_like(X) return env_sum / valid_dirs

这段代码就是上面第二节说的包络均值计算的直接落地。每个方向投影找局部极大值，然后在极值点的位置上对每个原始通道做三次样条插值。注意CubicSpline用的边界条件是bc_type='natural'，自然样条在边界处二阶导为零，对时间序列的端点效应会比默认的“非扭结”条件温和一些。

valid_dirs变量做了个保护，万一某条方向投影的极值点太少，无法构造样条，这个方向就直接跳过。这在实际数据里很常见，尤其是高阶IMF分离之后，信号越来越平缓，极值点越来越少。

3.3 单次筛选与完整MEMD主循环

def _sift_one_imf(X, directions, max_sift=80, sd_thresh=0.01): """筛出一个多元IMF，返回筛选后的细节分量""" n_ch, n_samples = X.shape residue = X.copy() t = np.arange(n_samples) for _ in range(max_sift): mean_env = _mean_envelope(residue, directions) if np.all(np.abs(mean_env) < 1e-12): break detail = residue - mean_env num = np.sum(mean_env ** 2) den = np.sum(residue ** 2) + 1e-12 criterion = num / den residue = detail if criterion < sd_thresh: break return detail def memd(X, n_dirs=None, max_imfs=10, sd_thresh=0.01, max_sift=100): """ 多元经验模态分解 参数 ---- X : shape = (n_channels, n_samples) 输入多元信号 n_dirs : int 方向向量数量，None时按通道数自动估算 max_imfs : int 最多分解出的IMF数量 sd_thresh : float 筛选停止的能量比值阈值 max_sift : int 单次筛选最大迭代次数 返回 ---- imfs : shape = (n_imfs, n_channels, n_samples) 分解出的全部IMF residual : shape = (n_channels, n_samples) 剩余趋势分量 directions : shape = (n_dirs, n_channels) 本次使用的方向向量 """ X = np.atleast_2d(X).astype(np.float64) n_ch, n_samples = X.shape if n_dirs is None: n_dirs = max(64, 16 * n_ch * (n_ch + 1)) directions = unit_vectors_on_hypersphere(n_dirs, n_ch) imfs = [] residual = X.copy() for _ in range(max_imfs): imf = _sift_one_imf(residual, directions, sd_thresh=sd_thresh, max_sift=max_sift) residual = residual - imf imfs.append(imf) total_extrema = 0 for c in range(n_ch): total_extrema += len(_local_maxima(residual[c])) if total_extrema < n_ch * 2: break return np.array(imfs), residual, directions

主流程比单通道EMD还要简洁，因为关键逻辑都收进_mean_envelope里了。每次筛出一个IMF后，从残差中减去，然后检查剩余信号是否还有足够多的极值点来判断要不要继续。注意这里判断残余信号是否还能继续分解，用的是每个通道的局部极大值个数之和。当极值点太少，样条包络已经没有意义，强制停止。

默认的n_dirs估算规则是我自己的经验值：通道数低时至少64个方向，通道数高时按16 * n_ch * (n_ch + 1)放大，比如3通道大概192个方向，8通道就是1152个方向。这个规则不一定最优，但方向数太少会影响包络质量，方向数太多又太耗时，算是一个比较平衡的起点。

3.4 先做单通道退化验证

写算法有个习惯，先用最简单的场景验证实现没有结构性错误。单通道时MEMD应当接近经典EMD的思路，我们用一条合成信号来测：

t = np.linspace(0, 1, 500) x = np.sin(2 * np.pi * 5 * t) + 0.4 * np.sin(2 * np.pi * 40 * t) x = x.reshape(1, -1) imfs, residual, _ = memd(x, n_dirs=2) print("IMF数量:", imfs.shape[0]) print("每个IMF的近似频率幅度比:", [ np.round(np.max(np.abs(np.fft.rfft(imfs[i][0]))), 3) for i in range(imfs.shape[0])])

单通道时我把方向数改成2，因为球面上只有两个方向：+1和-1。正方向包络过极大值，负方向包络过极小值，两者平均效果和经典EMD一致。跑下来应该能得到两个IMF，第一个能量集中在40Hz附近，第二个集中在5Hz附近，残余是一小段接近零的趋势。如果你把这里的IMF波形和某个成熟库的单通道EMD结果对比，会发现主要振荡成分基本一致，微小差异来自样条插值边界条件和极值点查找细节，属于正常偏差。

4. 三通道仿真实验：验证模态对齐效果

算法写完了，得拿一个有说服力的例子测一下。我构造了一个三通道仿真信号，每个通道共享同一个3Hz低频振荡，但各自叠加了不同频率、不同幅度的高频成分：

fs = 200 T = 2.0 N = int(fs * T) t = np.arange(N) / fs np.random.seed(42) x1 = 1.0 * np.sin(2 * np.pi * 3 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.05 * np.random.randn(N) x2 = 1.2 * np.sin(2 * np.pi * 3 * t + 0.8) + 0.4 * np.sin(2 * np.pi * 40 * t) + 0.05 * np.random.randn(N) x3 = 0.8 * np.sin(2 * np.pi * 3 * t + 1.2) + 0.3 * np.sin(2 * np.pi * 25 * t + 0.5) + 0.05 * np.random.randn(N) X = np.vstack([x1, x2, x3]) imfs, residual, directions = memd(X, n_dirs=128, max_imfs=8, sd_thresh=0.005) print("IMF shape:", imfs.shape)

这个信号的设计意图很明确：三个通道有共同的3Hz成分，这是验证“对齐”的关键；每个通道又各有一个几十赫兹的高频成分，用来考验高频成分能否被正确分离到各自通道的IMF里，而不是互相混叠。

实际分解结果我描述一下，方便你验证。以n_dirs=128、sd_thresh=0.005跑，通常得到3到4个IMF加一个残余。第一个IMF主要包含各通道的高频成分，但通道1的50Hz、通道2的40Hz、通道3的25Hz会在各自的通道分量里保持独立，互不串扰；第二个或第三个IMF会集中承载3Hz附近的低频振荡，三个通道的IMF波形在时间上高度同步，相位差保持在初始构造时的0.8和1.2弧度附近；最后的残差基本是一条直流分量或者微弱趋势线。

要量化“对齐”效果，一个直接办法是看三个通道同一阶IMF之间的峰值位置是否在同一时刻出现。在3Hz这个IMF里，三个通道的极大值时间索引应该非常接近，偏差在1个采样点以内。而如果你拿每个通道单独做EMD，大概率会出现某个通道的3Hz成分落在IMF2、另一个落在IMF3的尴尬局面，这就是MEMD相对逐通道处理最直观的优势。

绘图代码可以这样写，方便对照波形：

import matplotlib.pyplot as plt plt.figure(figsize=(10, 8)) for i in range(min(3, imfs.shape[0])): plt.subplot(5, 1, i + 1) plt.plot(t, imfs[i][0], label='ch1') plt.plot(t, imfs[i][1], label='ch2') plt.plot(t, imfs[i][2], label='ch3') plt.legend(loc='upper right') plt.xlim([0, 1.5]) plt.subplot(5, 1, 4) plt.plot(t, residual[0], label='residual ch1') plt.plot(t, residual[1], label='residual ch2') plt.plot(t, residual[2], label='residual ch3') plt.legend(loc='upper right') plt.tight_layout() plt.show()

跑完看到的波形就像我上面说的那样：高频IMF里三个通道各有自己的波形，低频IMF里三条曲线紧密贴在一起，只是幅度和相位不同。这个“低频振荡紧密同步”的效果，正是MEMD后续做跨通道时频分析的基础。

5. 调参、踩坑与不同场景的取舍经验

5.1 方向数K到底怎么选

方向数K是MEMD最主要的超参数。它影响包络质量，也直接决定计算时间。用n_dirs=64和三通道组合时，实测多数合成信号都能分解得不错；但遇到通道间幅度差异特别大的数据，比如一个通道振幅是另一个的10倍以上，64个方向可能不够，包络均值会被强通道带偏，弱通道的模态就容易混叠。这时候把方向数提到128或256，效果会明显改善。

维度很高的数据，比如几十个通道同时处理，方向数需求会急剧上升。因为高维球面上的均匀覆盖需要更多采样点，几百个方向都是常态。但方向数增大带来的是K倍的计算量，所以基本策略是：先用较小的K快速验证逻辑，确认模态结构稳定后再逐步增大K。如果某次增大K后分解结果变化很小，说明当前方向数已经够用。

5.2 样条插值的过冲问题

三次样条插值在信号变化剧烈的地方容易出现过冲，比如尖峰附近的插值曲线会甩出一个超过实际极值幅度的“大波浪”，这种包络过冲会直接污染IMF。低阶IMF（高频振荡）比较敏感，因为包络在几十赫兹上快速摆动，样条更容易过头。

缓解办法有两个。一是换用PCHIP插值，分段三次Hermite插值多项式能保证插值结果在节点间单调，不会产生新的极值点。scipy里对应的是PchipInterpolator。在实际项目中处理生理信号或振动信号时，我经常用PCHIP替代CubicSpline，代价是包络曲线没那么光滑，但换来了稳定。代码改动就一行：

from scipy.interpolate import PchipInterpolator cs = PchipInterpolator(t[maxima], X[c, maxima])

二是把插值极值点的边界条件处理得更保守一些。自然样条已经比非扭结条件稳，但如果你发现边界IMF变形特别厉害，可以考虑直接在包络计算前给信号两端加一段镜像数据，用完再裁掉。端点效应这个坑下面单独说。

5.3 停止阈值和最大筛分次数

sd_thresh是筛选停止阈值。我在合成信号上对比过不同取值的差异：0.05时低频IMF里还能明显看到高频残留，0.001时高频IMF被反复磨削，噪声成分被当作信号抠出来。多数情况下0.005到0.01是一个比较合理的区间。如果你关注的是低频慢波，可以稍微放宽阈值，因为低频成分对高频噪声不那么敏感；如果你要分析高频细节，阈值就得收紧，同时把max_sift适当调大，比如150。

max_sift设太大还有一个隐患：筛分次数太多会让IMF变成纯正弦式的窄带信号，丢失原始波形的非线性特征。EMD类算法本来的优势之一就是保留信号的形态特征，筛分过头反而弄巧成拙。我通常限制在100次以内，看到收敛不了就改阈值而不是穷举迭代。

5.4 边界效应与数据长度

端点效应是所有基于包络插值算法的通病。信号两端没有足够的极值点，样条包络在边界会发散或剧烈摆动，导致刚开始的十几二十个采样点上的IMF不可信。MEMD的高维投影和平均过程会让边界问题更复杂，因为每个方向投影的端点情况都不同。

最实用的处理方式不是试图把边界完全修好，而是接受失真并留出余量。分解前把信号两头各延拓一段，比如镜像延拓或线性预测延拓，分解后再把延拓部分裁掉；或者分析时直接忽略每个IMF首尾5%的数据点。尤其在做谱分析时，不建议把边界那一小段算进结果，后面还得花时间清理边缘伪影。

5.5 模态混叠与噪声辅助思路

MEMD能缓解逐通道EMD的模态错位，但并不能根治模态混叠。如果信号里混有间歇性强振荡，比如时有时无的50Hz脉冲，分解时低频模态里还是可能嵌进高频片段。这类问题的常规解法是往MEMD里引入噪声辅助机制，思路类似于单通道领域的集成EMD思路：给多元信号加多次白噪声，重复MEMD分解，再把各次结果取平均。这样能把间歇成分“扯”到合适的IMF里，代价是计算量成倍上涨。本次给出的代码没有包含这层，因为核心逻辑已经足够复杂，先跑通基础版再接噪声辅助也不迟。

5.6 单通道场景不要滥用MEMD

最后提醒一句：如果数据本质上只有单通道，直接用经典EMD就好。MEMD在单通道时虽然能退化，但毕竟是把问题绕了一圈，边际效益不高。它的主战场一定是多通道联合分析，尤其是G通道数较多、通道间存在相关性、需要统一分解框架再往下做时频统计的场景。方向数、阈值这些参数的纠错成本，也只有在多通道数据上才算花得值。

我在实际使用中最深的体会是：MEMD的价值不只在“分解”本身，更在“对齐”带来的后处理便利。只要IMF对齐了，后面做希尔伯特谱、通道间相干性、事件相关振荡分析，都会顺畅很多。现在手头有多通道信号处理需求的朋友，不妨拿这份代码先跑起来，再根据自己的数据特性调方向数和阈值。过程里如果看到某阶IMF波形突然发飘，多半是边界题或样条过冲，往回翻第五小节就能找到对策。