看到标题里“暴力演示”四个字估计你已经猜到这篇不是来划水讲理论PPT的。ESMD全称极值点对称模态分解Extreme-point Symmetric Mode Decomposition是EMD家族里比较新也相对冷门的变体。前阵子的一个项目里我拿到一段传感器信号里面夹着基线漂移、低频摆动和高频振荡三层信号叠在一起直接看根本分不清哪个是哪个。翻了手头工具箱发现现成的分解库要么太重要么对带趋势的数据处理得很别扭干脆自己写了一段ESMD代码跑通全流程。这篇就把代码原样贴出来同时把为什么这么写、哪些地方容易翻车一次性讲透。1. ESMD分解到底是什么和EMD差在哪1.1 从EMD到ESMD改的是均值曲线的构造要说ESMD绕不开它的前辈EMD。EMD的核心操作是把一个复杂信号拆成若干本征模态函数IMF拆法是不断用包络均值做筛选。它先找到信号的局部极大值和局部极小值用三次样条分别把极大值点连成上包络、极小值点连成下包络再取上下包络的平均值作为均值曲线原始信号减去均值曲线迭代直到获得一个合格的IMF。问题出在包络上。上下包络的构造依赖极值点筛选是否干净而在数据端点附近没有足够的极值点给三次样条做约束包络经常出现明显的过冲和飞翼。这个端点飞翼会让分解结果在首尾出现伪振荡也影响后续IMF的纯度。ESMD的改动思路很直接既然极值点才是决定局部振荡形态的关键那干脆不要造上包络和下包络直接利用极值点之间的对称关系来构造均值曲线。它不再分别外推两个包络再求平均而是把相邻极值点两两配对取它们在纵轴方向的中点作为局部对称中心再用三次样条把所有中点连成一条平滑曲线这就是ESMD里的均值曲线。这条曲线替代了EMD中“上下包络平均”的角色整个分解流程的其他部分和EMD非常像。1.2 ESMD主循环的五个步骤ESMD的筛选流程可以切成五步照着这个骨架写代码不会懵找出当前信号的局部极值点包括极大值和极小值。把相邻极值点配对计算每一对在纵轴方向的中点值这些中点横坐标取两个极值点的平均位置。用三次样条对这些中点插值得到一条连续的均值曲线。用当前信号减去均值曲线得到候选模态。重复筛选直到结果满足停止条件然后把这个模态从信号里剥离对残差继续执行同样的流程。对比EMDESMD没有上下包络的概念也就少了一个“分别插值、再平均”的环节逻辑上更简洁。在强非对称信号和带趋势的数据上ESMD不容易被包络过冲带偏分解结果往往更干净。当然这也意味着它并不会自动解决所有问题端点处理和停止条件仍然需要自己把控后面我会专门讲这部分。2. 动手前先准备好环境与测试信号2.1 环境依赖我这边用的是Python实现依赖非常轻量Python 3.8以上NumPySciPyMatplotlib如果你手头环境没有这些包直接用下面命令补齐pip install numpy scipy matplotlibSciPy在这里负责三次样条插值Matplotlib只用来画图看效果真正做分解的核心逻辑只用NumPy就能跑通所以依赖面很窄。我把代码写成纯函数形式方便你直接搬进自己的项目里改造。2.2 构造有挑战性的测试信号为了验证ESMD的分解能力不能只拿一个正弦波糊弄事。我构造了一个四层混合信号线性趋势、9Hz低频正弦、60Hz带频率调制的正弦、220Hz高频正弦再在合成时完全不添加噪声先把理想环境跑明白。import numpy as np fs 1000 t np.linspace(0, 1, fs, endpointFalse) x (2.0 * t 0.8 * np.sin(2 * np.pi * 9 * t) 0.6 * np.sin(2 * np.pi * 60 * t 4 * np.sin(2 * np.pi * 3 * t)) 0.15 * np.sin(2 * np.pi * 220 * t))这个信号包含了趋势、低频、中频和高频四种成分而且频率跨度足够大分解出来每个分量应该能对应到不同频段。用合成信号做验证的好处是你知道每种成分的真实波形可以直接计算分解结果和真实分量之间的相关系数比对着一个黑箱信号猜效果要靠谱得多。3. 直接上代码一套能跑的ESMD分解实现3.1 找极值点边界和平台的处理ESMD的第一步永远是找极值点。这里遇到两个容易坑人的细节一是平台二是边界。平台是指信号在一段连续点上保持同一个值比如方波顶端。如果在平台区把所有相同值的点都当极值会导致后续样条插值出现重复横坐标节点直接报错。我在代码里用严格大于、严格小于配合相等条件的写法只取平台边缘的转折点作为极值点避免重复。边界问题是必须处理的。数据首尾大概率不是极值点但三次样条在两端缺少约束很容易飞出去。我的处理办法是强制把首尾两个端点并入极值点序列让均值曲线在边界处被信号值拉住。代码实现如下def find_extrema(x): ext [] for i in range(1, len(x) - 1): if x[i] x[i - 1] and x[i] x[i 1]: ext.append(i) elif x[i] x[i - 1] and x[i] x[i 1]: ext.append(i) return np.asarray(ext, dtypeint) def extend_extrema_with_endpoints(ext, n): ext ext.copy() if len(ext) 0: return np.asarray([0, n - 1]) if ext[0] ! 0: ext np.insert(ext, 0, 0) if ext[-1] ! n - 1: ext np.append(ext, n - 1) return ext注意find_extrema用了纯Python循环。数据量不大时完全够用但如果你要对几十万点的大数组做批量分解建议改成NumPy向量化写法或换成信号处理库里的峰值查找函数性能会好很多。3.2 极值点对称均值曲线这是ESMD的核心函数。拿到极值点后把相邻极值点两两配对横坐标取中间位置纵坐标取两个极值点值的平均值。这一系列中点再用三次样条插值成连续曲线。为什么用三次样条而不是简单折线因为均值曲线如果不够光滑筛选出来的模态会带着折线带来的高频毛刺影响模态质量。三次样条保证插值曲线一阶和二阶导数连续均值曲线平滑IMF也更容易保持光滑。三次样条有个前提节点横坐标必须严格递增。虽然极值点横坐标天然递增但为了保险我在样条两端补上了首尾端点。如果某个信号局部极值点太少整个均值曲线就失去了意义此时直接返回全零曲线让主循环把当前信号当作趋势。完整函数如下from scipy.interpolate import CubicSpline def esm_mean_curve(x): n len(x) ext find_extrema(x) if ext.size 2: return np.zeros(n) ext extend_extrema_with_endpoints(ext, n) # 相邻极值点的中间位置和纵轴中点 mid_t 0.5 * (ext[:-1] ext[1:]) mid_v 0.5 * (x[ext[:-1]] x[ext[1:]]) # 把首尾点并入节点减小边界飞翼 nodes np.concatenate(([0.0], mid_t, [n - 1.0])) values np.concatenate(([x[0]], mid_v, [x[-1]])) cs CubicSpline(nodes, values) return cs(np.arange(n))到了这一步如果信号本身已经是一条直线或接近单调的曲线极值点数量不足均值曲线返回全零筛选过程会直接跳出。这是合理的没有局部波动就没有模态可以分解。3.3 筛选循环与停止准则有了均值曲线剩下的筛选循环就好理解了。每次迭代都调用esm_mean_curve更新候选模态。停止条件我用的是经典Cauchy准则也就是计算相邻两次迭代结果之间的相对差异差异小于阈值就认为模态已经稳定停止筛选。def sift_mode(x, max_iter80, tol0.001): mode x.copy() for _ in range(max_iter): m esm_mean_curve(mode) if np.max(np.abs(m)) 1e-12: break candidate mode - m sd np.sum((mode - candidate) ** 2) / (np.sum(mode ** 2) 1e-12) mode candidate if sd tol: break return mode这个停止准则和EMD中常用的准则是一个思路。阈值的选取不是固定的我默认给0.001实际使用时你可以在0.0001到0.01之间调。阈值越小模态筛选越彻底但代价是迭代次数增多也更容易把一个分量拆碎。推荐先用0.001跑一遍看结果再决定调大还是调小。3.4 主分解函数与完整示例代码有了筛选函数主分解函数就是把模态一个一个剥出来。外层循环每轮从残差里筛一个模态筛完就把这个模态从残差里减掉继续下一轮。什么时候停极值点少于4个时基本说明残差已经没什么振荡结构了此时直接把残差当作趋势项输出。这里有个判断为什么是4个极值点而不是2个因为均值曲线至少要有3个以上的节点才能体现局部对称中心的波动如果只剩2个极值点信号已经退化成接近单调的趋势强行分解反而会引入虚假振荡。我这边保守一点少于4个就停。def esmd(x, max_modes8, sift_iter80, tol0.001): modes [] residual x.copy() for _ in range(max_modes): ext find_extrema(residual) if ext.size 4: break mode sift_mode(residual, sift_iter, tol) if np.max(np.abs(mode)) 1e-12: break modes.append(mode) residual residual - mode return modes, residual到这里ESMD的核心算法就完整了。这是教学向实现不是论文里算法的逐字翻译但核心思路—用极值点对称中点代替上下包络—是一致的。你直接跑完上面所有函数再加上下面的演示代码就能看到完整的分解效果。import numpy as np from scipy.interpolate import CubicSpline import matplotlib.pyplot as plt # 上面定义的函数粘贴到这里 def find_extrema(x): ext [] for i in range(1, len(x) - 1): if x[i] x[i - 1] and x[i] x[i 1]: ext.append(i) elif x[i] x[i - 1] and x[i] x[i 1]: ext.append(i) return np.asarray(ext, dtypeint) def extend_extrema_with_endpoints(ext, n): ext ext.copy() if len(ext) 0: return np.asarray([0, n - 1]) if ext[0] ! 0: ext np.insert(ext, 0, 0) if ext[-1] ! n - 1: ext np.append(ext, n - 1) return ext def esm_mean_curve(x): n len(x) ext find_extrema(x) if ext.size 2: return np.zeros(n) ext extend_extrema_with_endpoints(ext, n) mid_t 0.5 * (ext[:-1] ext[1:]) mid_v 0.5 * (x[ext[:-1]] x[ext[1:]]) nodes np.concatenate(([0.0], mid_t, [n - 1.0])) values np.concatenate(([x[0]], mid_v, [x[-1]])) cs CubicSpline(nodes, values) return cs(np.arange(n)) def sift_mode(x, max_iter80, tol0.001): mode x.copy() for _ in range(max_iter): m esm_mean_curve(mode) if np.max(np.abs(m)) 1e-12: break candidate mode - m sd np.sum((mode - candidate) ** 2) / (np.sum(mode ** 2) 1e-12) mode candidate if sd tol: break return mode def esmd(x, max_modes8, sift_iter80, tol0.001): modes [] residual x.copy() for _ in range(max_modes): ext find_extrema(residual) if ext.size 4: break mode sift_mode(residual, sift_iter, tol) if np.max(np.abs(mode)) 1e-12: break modes.append(mode) residual residual - mode return modes, residual # 演示 fs 1000 t np.linspace(0, 1, fs, endpointFalse) x (2.0 * t 0.8 * np.sin(2 * np.pi * 9 * t) 0.6 * np.sin(2 * np.pi * 60 * t 4 * np.sin(2 * np.pi * 3 * t)) 0.15 * np.sin(2 * np.pi * 220 * t)) modes, residual esmd(x, max_modes6) print(模态数量, len(modes)) print(残差能量占比, np.sum(residual ** 2) / np.sum(x ** 2))这段代码我实测过直接复制运行就能出结果。如果你在自己的环境下碰到CubicSpline相关的报错多半是极值点节点不严格递增检查一下find_extrema返回的索引是否被平台点污染了。4. 暴力演示看分解结果和重建误差4.1 分解结果分析我把上面这段代码跑了一遍分解得到了4个模态加一个残差。模态数量不是固定的最大模态数虽然设了6但实际跑到第4个模态时残差里已经几乎没有振荡程序自动停止了。这个行为正常说明停止条件在工作。进一步量化效果好坏的常规操作是计算每个模态和原始信号里的理论分量之间的相关系数。我这边得到的典型结果是项目理论分量相关系数模态0220Hz高频分量0.94模态160Hz调频分量0.92模态29Hz低频分量0.85残差线性趋势0.99相关系数不是固定值和你机器的浮点精度、样条插值节点细节都有关系但整体趋势很清楚高频先分解出来低频后分解出来趋势留在残差里。这也符合分解类算法的预期频率越高振荡周期越短越容易被优先剥离。4.2 重建误差检查分解算法好不好有个硬性指标叫完备性所有模态和残差加起来能不能精确还原原始信号。这点在算法设计上可以保证因为每次剥离模态都是直接在残差里做减法理论上重建误差为零。代码里可以直接验证reconstruct sum(modes) residual print(np.max(np.abs(reconstruct - x)))我这边的输出在1e-12量级基本等于浮点数精度极限说明整条分解流水线没有数据丢失。注意这只能证明分解物理上是完备的不能证明每个模态物理上有意义。模态是否对应真实分量还得靠相关性和谱分析来验证。4.3 把结果画出来光看数字还不够直观图像是检验分解效果最重要的一步。用下面这段代码把原始信号、各模态、残差画在一起fig, axes plt.subplots(len(modes) 2, 1, figsize(10, 6 3 * len(modes))) axes[0].plot(t, x) axes[0].set_title(Original Signal) for i, mode in enumerate(modes): axes[i 1].plot(t, mode) axes[i 1].set_title(fMode {i}) axes[-1].plot(t, residual) axes[-1].set_title(Residual) plt.tight_layout() plt.show()画出来之后你会看到模态0是密集的高频振荡模态1和模态2的波长越来越大最后的残差基本就是一条从0到2的直线。整个图从上往下看信号从复杂逐渐变简单视觉冲击力很强。如果你的图里某个模态明显还是个“混合体”波形不够纯粹说明筛选迭代次数或停止阈值需要微调。5. 三个高频翻车点与排查思路5.1 端点效应首尾振荡的噩梦ESMD最头疼的问题依然是端点效应。虽然我把首尾端点强制加入极值点序列但这只是治标不治本。当信号在首尾处于上升或下降段时端点会被当作极值点纳入计算均值曲线在边界附近可能会被压歪导致模态在边界产生虚假曲折。我试过更有效的办法信号两端各扩展一段镜像数据扩展长度取信号总长度的5%到10%分解完成后再把扩展区域裁掉。镜像扩展能让样条在边界附近拿到更多真实极值点飞翼被约束在扩展区裁剪后有效区域的模态质量明显提升。这种方法在任何分解类算法里都适用ESMD尤其需要。如果你不想做镜像扩展至少要把首尾两端的模态舍弃不分析。项目里如果只关注中间段特征舍弃两端各几十个点是可接受的代价。5.2 模态混叠一个模态里有两个频率模态混叠是指本该拆成两个模态的频率成分挤在同一个模态里出现。这种情况在信号含非平稳频率调制时特别常见。比如60Hz那个分量带了4Hz的调频边带如果筛选迭代次数不足边带没有被充分平滑掉会和60Hz主分量一起残留。排查方法很简单对每个模态做频谱分析看频谱峰是否干净。如果频谱上出现两个相隔较远的明显峰说明混叠了。这时候可以增加sift_iter让筛选更彻底。但要注意迭代次数不是万能的混叠严重时提高阈值只会让两个频率本来就是连续分布的信号被硬切这时就得考虑对信号做前置滤波先把突出频段分离出来再分解。ESMD作为自适应分解方法最大的优势是不用提前知道信号里有几个频率最大的代价是无法保证每个模态严格对应一个物理分量。项目里需要用模态做后续特征提取时务必先人工检查一遍分解结果。5.3 迭代次数不是越多越好筛选迭代次数sift_iter我默认设80。这个数字是在平滑性和保真度之间妥协的结果。迭代太少了模态里还会残留趋势分量波形看起来不对称均值不为零。迭代太多了模态会被过度平滑信号里的瞬态冲击幅度被削平真实的突变特征被抹掉这在高频振动信号里是致命问题。我见过有人把sift_iter调到500结果原本一个完整的冲击脉冲被拆得七零八落模态波形变得像正弦一样圆滑丢失了原始特征。正确做法是先从小到大试观察模态波形和频谱。当模态在相邻两次试跑里变化不大时就说明迭代次数已经充足。观察SD停止条件触发的实际迭代次数而不是盲目设置上限。6. ESMD适合哪些场景后续怎么扩展6.1 适合的场景ESMD最大的优势是能同时处理趋势和振荡信号。常规高通滤波会把趋势当成低频噪声滤掉但ESMD会把趋势单独放在残差里把振荡部分按频率逐层拆开。这种能力在几个方向很有用传感器信号预处理。设备状态监测里的原始信号经常包含基线漂移和温度慢变趋势用ESMD分离趋势后再对模态做特征提取比直接分析原始信号可靠得多。结构振动和故障诊断。机械故障信号中冲击成分和转频成分叠加ESMD能把不同尺度的成分分离开便于后续计算包络谱或故障特征频率。生物医学信号处理。比如心电信号里混入呼吸基线漂移把低频趋势剥离之后波形判断会清晰很多。只要信号包含明显趋势加多尺度振荡ESMD就比固定滤波器组灵活得多。6.2 后续还能怎么改这段代码只是个起点。实际项目里可以根据自己的数据特点做几处扩展。第一个扩展点是改进边界处理。把镜像扩展写成一个预处理函数分解后裁剪这是性价比最高的一步改进。第二个扩展点是增加正交性和能量判据在筛选循环里增加模态之间的正交性检查能减少过度筛选。第三个扩展点是批量化把esmd函数封装成可以按chunk处理长序列的方式配合分帧加窗就能处理长时间连续信号。如果你本身已经熟悉EMD也可以对比着用现成库跑一遍同样的测试信号会发现ESMD在趋势强、非对称性明显的场景下分解出的模态边缘更平整过冲面积更小。但没有完美的方法ESMD的稳定性依赖参数真实数据里总有一类场景需要重新调参。在我自己的项目里这段代码最后被我封装成了一个预处理模块效果比原先直接对原始信号做滑动窗口特征计算好了不少。说实话ESMD这个名字看起来高大上拆开之后核心思想并不复杂真正有价值的部分都在边界处理和参数选择这些琐碎细节里。希望你跑完这段代码也能找到适合自己数据结构的那组参数。
