对很多刚开始接触自适应信号分解的朋友来说,第一次听到“极点对称模态分解”这个名字,多半是在搜索EMD、EEMD相关的资料时被顺带提及的。我最初接触ESMD,是因为处理一组轴承振动信号时,EMD分解出的IMF在端点处出现了明显的“飞翼”,导致后续的时频分析出现了一大片假频率。后来换用极点对称模态分解(ESMD)重跑,分解结果干净了不少,这才真正把注意力放到这个相对小众却相当实用的算法上。
这篇文章围绕ESMD的MATLAB实现展开,讲清楚它的核心思路是什么、和经典EMD的差异在哪里、完整代码怎么写、跑出结果之后又该怎么判断好坏。适合正在做信号处理相关课题的研究生、需要处理非平稳数据的工程师,以及想跳出EMD全家桶、扩展自己工具箱的技术人员。读完你不仅能跑通一个完整的ESMD分解Demo,还能知道哪些参数真正影响结果、哪些环节最容易踩坑。
1. 模态分解家族里,ESMD到底改了什么
1.1 EMD的“包络之痛”
要理解ESMD的价值,得先从EMD的痛点说起。经典EMD假设信号由一系列固有模态函数(IMF)叠加而成,分解时先找信号的局部极大值和局部极小值,用三次样条分别插值成上包络和下包络,取上下包络均值作为低频中心线,再用原始信号减掉这个均值,反复迭代得到一个IMF。
这个流程看起来合理,但实际运行时会遇到几个很现实的问题。三次样条插值对端点附近的数值非常敏感,极值点稍微稀疏一点,包络就会在大幅摆动,产生明显的过冲和欠冲。更麻烦的是,整个包络是整个信号空间的全局拟合,某个局部极值的变化会通过样条基函数影响到其他区域,导致分解结果对局部扰动异常敏感。很多EMD使用者都遇到过这样的场景:信号只是多了一个毛刺,分解出的第一个IMF却整体变了形,模态混叠随之出现。
也正是因为这些问题,后续才催生了EEMD、CEEMDAN、VMD等一系列改进算法。EEMD和CEEMDAN的思路是在信号中添加白噪声辅助分解,以牺牲计算量为代价换取模态稳定性;VMD则干脆换了个赛道,把分解问题转化为变分约束优化问题。而ESMD走的是另一条路——保留EMD的自适应分解框架,但把最核心的“包络均值”替换成一种更稳健的“内部极点插值”方式。
1.2 ESMD的改进思路:内部极点插值
ESMD的“极点对称”体现在哪里?看分解流程就清楚了。
第一步,找信号的全部局部极值点,包括极大值和极小值。第二步,把相邻极值点用直线段连接起来。比如第一个极值点和第二个极值点连一条线,第二个和第三个连一条线,这样就把整个信号区间切成了一段段首尾相接的折线。第三步,取每段折线的中点,再对这些中点做插值或平滑处理,得到一条贯穿信号的低频中心线,作为这一轮筛分需要减掉的部分。
与EMD对比,关键差异就在这里。EMD必须构造“上包络”和“下包络”,再取均值;ESMD完全不区分上下包络,它只用相邻极值点的信息构建一条内插线,天然避开了外包络外插带来的数值震荡,同时计算量也比三次样条包络拟合小很多。ESMD名称中的“极点对称”,指的就是这种以极值点为锚点、按对称中点削减高频分量的策略。
用生活化的方式理解两种算法的差别:EMD像是给信号画一条“上限线”和一条“下限线”,然后取中间那条线;ESMD则像是直接沿着极值点之间拉皮筋,皮筋中段自然下垂的位置就是需要去除的低频成分。后者更贴近信号的局部形态,也就更少受到远端信号干扰。
1.3 与其他分解方法的横向对比
实际选型时,很多人纠结到底用EMD、EEMD还是ESMD。我把常用几种方法的特性整理成一张表,方便对照。
| 方法 | 核心机制 | 主要参数 | 模态混叠 | 端点效应 | 计算开销 | 适用场景 |
|---|---|---|---|---|---|---|
| EMD | 上下包络均值 | 筛分次数、停止阈值 | 较明显 | 明显 | 低 | 非平稳信号快速分解 |
| EEMD | EMD+白噪声平均 | 噪声幅值、集成次数 | 有明显改善 | 仍存在 | 很高 | 含噪声信号,模态分解 |
| CEEMDAN | EEMD的完备性改进 | 噪声幅值、迭代次数 | 改善较大 | 仍存在 | 很高 | 要求可重构的分解场景 |
| ESMD | 内部极点插值 | 最大筛分次数 | 改善较大 | 较轻 | 中 | 端点效应敏感的实测信号 |
| VMD | 变分约束优化 | 模态数K、惩罚因子 | 不依赖端点 | 较轻 | 中高 | 基本可视为窄带信号 |
这表不是要说明ESMD全面优于其他算法,而是提供一个选择依据。如果你的信号本身包络形态规则、端点较少,EMD完全够用;如果追求分解的完备性和抗噪能力,CEEMDAN或VMD可能更合适。但当你面对的是采样时长不长的实测信号,端点的几个点往往携带关键特征时,ESMD的“内部极点插值”思路就格外有优势。
2. 拆ESMD算法流程:从极点筛选到内部极值插值
2.1 极值点检测与端点处理策略
整个ESMD的第一环是极值点检测,这一步的准确性直接决定后续所有结果。MATLAB里找极值点最常用的方式是diff符号法:对信号做一阶差分,当差分值从正变负,说明信号经过一个局部极大值;从负变正,则是局部极小值。
考虑一段可能存在连续相同值的数据,比如传感器采到的平台信号。如果直接比较相邻点,平台上的多个点可能被误判成多个极值点,把后续的线段插值搅乱。稳妥的做法是,把符号变化的位置定位在平台中部,或者先对信号做微小平滑预处理。
function [idxMax, idxMin] = findExtrema(x) % 极点检测,返回极大值和极小值下标 n = length(x); dx = diff(x); % 极大值:由上升转为下降,允许平台 idxMax = find(dx(1:n-2) > 0 & dx(2:n-1) <= 0) + 1; % 极小值:由下降转为上升,允许平台 idxMin = find(dx(1:n-2) < 0 & dx(2:n-1) >= 0) + 1; end这里用dx(1:n-2)和dx(2:n-1)而不直接对整个diff取符号,是为了捕捉“从上升到持平再到下降”这类平台型极值。实际操作中,我还会在进入极值检测前先做一次去均值处理,避免直流分量影响过零判断。
端点处理是这一环节的另一个关键点。ESMD整体对外插依赖比EMD低,但不代表可以完全不处理端点。如果信号在端点处正好处于一个大幅振荡的相位,最初几段线段插值会缺少边界极值的约束,导致分解出的首个IMF前端出现小幅弯曲。我习惯在分解前用镜像延拓的方式把信号朝两端各扩展几个极值点周期的长度,分解完成后再裁剪回原始长度。
2.2 相邻极点线段构造与中点提取
ESMD的核心步骤,就是把找到的极值点按序连接成折线,再取折线中点。这一段写成MATLAB代码非常简短,但背后的逻辑值得展开讲。
function lowFreq = getInternalMidline(x, idx) % x 为输入信号,idx 为全部极值点下标(极大和极小合并排序) if length(idx) < 3 lowFreq = mean(x) * ones(size(x)); return; end % 极值点数值 peakVals = x(idx); % 相邻极值点连线中点 midVals = (peakVals(1:end-1) + peakVals(2:end)) / 2; % 中点对应的时间位置(取两点下标平均值) midIdx = (idx(1:end-1) + idx(2:end)) / 2; % 对中点序列做插值,还原为与原始信号等长的低频曲线 lowFreq = interp1(midIdx, midVals, 1:length(x), 'pchip'); end这段代码里有两个细节值得注意。
第一,为什么取中点之后还需要插值一次。相邻极值点的中点数量比极值点少一个,而且这些中点的时间位置不一定落在整数采样点上。要把它们变成与原始信号等长的曲线,就必须做插值。这里我选了pchip,也就是分段三次Hermite插值,它在保持单调性的同时不会出现普通三次样条那样剧烈的过冲。如果信号本身比较平滑,用spline也可以,但实测下来pchip在多数工程信号上更稳。
第二,lowFreq代表的是当前“平均信号”的近似低频分量,它不像EMD的包络均值那么“胖”,而是紧贴信号的中轴。因为所有线段都建立在相邻极值点之间,没有跨区域传递信息的通道,所以局部突变对低频曲线的影响范围也被限制在了两个相邻极值点之间。这就是ESMD对间歇性事件、脉冲成分更友好的原因。
2.3 筛分循环与自适应终止准则
有了低频中心线的计算方法,剩下的就是不断迭代筛分,直到提取出的高频分量满足IMF条件。ESMD的筛分循环和EMD非常相似,可以写成下面的结构。
function [imf, residue] = extractIMF(x, maxSweep) h = x(:); residue = x(:); for s = 1:maxSweep [idxMax, idxMin] = findExtrema(h); idx = sort([idxMax; idxMin]); if length(idx) < 3 break; end m = getInternalMidline(h, idx); h = h - m; % 判断是否满足IMF条件:过零数与极值点数量之差不超过1 zc = sum(diff(sign(h)) ~= 0); if abs(zc - length(idx)) <= 1 break; end end imf = h; residue = x - h; end筛分的停止条件有两个维度。第一个是标准IMF条件,即信号在整个时间范围内的过零点数量与极值点数量最多差一个。第二个是最大筛分次数maxSweep。实际信号噪声叠加后,IMF条件往往很难严格满足,如果任其迭代,可能会一直筛到结果变成纯调幅调频的“完美正弦波”,反而丢失真实物理成分。所以设置一个上限,迭代到达上限时强制停止,是工程上更稳健的选择。
外层分解循环则是反复调用extractIMF,每次从残差信号里提取一个IMF,直到残差信号成为单调趋势或极值点数量少于2。这个过程不再赘述,但有一个重要概念需要记住:ESMD分解出的最后一个残差项,不是噪声也不是废物,而是信号的高阶趋势项。第3部分会专门讨论怎么利用它。
3. 代码能跑了还不算完:筛选次数、终止判据和端点效应
3.1 筛选次数怎么定:自适应均方差准则
ESMD论文中提到了“自适应均方差”的概念,很多初次接触的人看到这个名词觉得很高深。其实思路很直接:筛分次数太少,IMF中残留大量低频成分;筛分次数太多,IMF被过度平滑,变成接近正弦波的形态。既然两个极端都不好,那就找一个让分解整体误差最小的折中值。
一个可落地的实现方式,是对一个信号预先尝试一组不同的筛分次数(比如从1到100),计算每次分解后重构信号与原始信号的均方误差,选取误差显著下降后开始平缓的拐点作为最佳筛分次数。虽然复杂度和直接跑一次分解相比变高了,但对于需要反复处理的同类型信号,这个预搜索的成本完全可以接受。
function bestS = searchBestSweep(x, sweeps) errs = zeros(size(sweeps)); for k = 1:length(sweeps) [imfs, residual] = esmd_basic(x, sweeps(k)); recon = sum(imfs, 2) + residual; errs(k) = mean((x(:) - recon(:)).^2); end [~, bestS] = min(errs); bestS = sweeps(bestS); end更简洁的做法是直接观察筛分次数与分解能量的关系。绘制筛分次数从10增加到200时各个IMF的能量曲线,你会发现能量在某个区间内快速变化,之后趋于稳定。那个拐点对应的筛分次数,就是这个信号的“工作点”。不同信号类型的工作点差异很大,我在处理风机振动数据时一般取50到100次,处理水文平稳序列时取20到30次就足够了。
3.2 端点效应对内插方案的影响
虽然ESMD不再做上下包络外插,但端点效应只是减轻,没有完全消失。原因很简单:信号最左端的极值点和最右端的极值点外侧,没有额外的极值点参与线段构造,所以两端各半个区间内的低频曲线只能靠插值外推或内侧线段信息填充。一旦信号端点附近存在大振幅波动,前几段中点的位置就会偏离真实中轴。
缓解端点效应最朴素的方法是镜像延拓。把信号开头一段以第一个极值点为镜面翻转,平移到信号左侧形成延伸段;结尾同样处理。这样原始端点处就有了新的极值点参与插值计算,分解完成后再把延拓部分裁掉。
function xExt = mirrorExtend(x, n) % 使用开头/结尾各 n 个采样点做镜像延拓 left = flipud(x(1:n)); right = flipud(x(end-n+1:end)); xExt = [left; x(:); right]; end镜像延拓的n不宜取太大,一般取信号一两个主周期对应的采样点数即可。n过大会引入远端信号的形态干扰,n过小则覆盖不了边界极值区间。在实际项目里,我通常先用快速傅里叶变换估算主频周期,再据此确定n。
3.3 残余模态的去趋势与物理意义
大多数信号处理教程教你分析IMF,却很少提醒你重视最后那个残差项R。ESMD对这个问题的看法是:残差项不是误差,而是信号的全局趋势或平均状态。
举个例子。分析一段桥梁应变监测数据,EMD分解后残差是一条单调上升的直线,很多同学直接把残差丢掉,只分析前几阶IMF,结果漏掉了桥梁支座沉降这个关键趋势。在ESMD框架下,残差项通常被解释为信号的趋势项,分析风场风速、电价序列这类包含明显背景变化的数据时,务必保留残差项并单独绘图观察。
另外,分解完成后做能量守恒检验也很关键。把各IMF和残差项相加,应当能重构出与原始信号几乎完全一致的波形。如果重构误差明显偏离零,说明分解过程中存在模态丢失或筛选过度的问题,需要回到筛分次数设置上重新调整。
4. 一个完整的MATLAB实战Demo:含噪信号分解与结果检验
4.1 构造测试信号与运行环境
下面给出一个可以直接跑通的完整例子。测试信号由三个分量叠加而成:一个频率为20Hz的正弦波、一个频率为60Hz的调幅波、一个在0.8秒到1.2秒之间出现的短暂脉冲,最后加上信噪比约15dB的高斯白噪声。对它做ESMD分解,观察算法能否把这三个不同时频特征的分量分离出来。
fs = 1000; t = (0:2*fs-1)/fs; x = sin(2*pi*20*t) + (1+0.5*cos(2*pi*5*t)).*sin(2*pi*60*t); x(t >= 0.8 & t <= 1.2) = x(t >= 0.8 & t <= 1.2) + 2*sin(2*pi*180*t(t >= 0.8 & t <= 1.2)); x = x + 0.15*randn(size(t));这段代码生成的数据长度2秒、2000个采样点。脉冲分量和调幅分量同时存在,能较全面地考验分解算法的分离能力。
4.2 ESMD分解主程序
为了方便大家直接复现,我把前面拆开的函数汇总成一个精简版ESMD主程序。代码不是论文官方版本的逐行翻译,而是按ESMD核心思想独立实现的可运行版本,对于理解算法原理和实际工程应用已经足够。
function [imfs, residual] = esmd_demo(x, maxSweep) x = x(:); imfs = []; residual = []; r = x; while true [imf, r_next] = extractIMF(r, maxSweep); imfs = [imfs, imf]; [idxMax, idxMin] = findExtrema(r_next); if length(idxMax) + length(idxMin) < 3 residual = r_next; break; end r = r_next; end end function [imf, residue] = extractIMF(x, maxSweep) h = x; for s = 1:maxSweep [idxMax, idxMin] = findExtrema(h); idx = sort([idxMax; idxMin]); if length(idx) < 3 break; end m = getInternalMidline(h, idx); h = h - m; end imf = h; residue = x - imf; end function lowFreq = getInternalMidline(x, idx) peakVals = x(idx); midVals = (peakVals(1:end-1) + peakVals(2:end)) / 2; midIdx = (idx(1:end-1) + idx(2:end)) / 2; lowFreq = interp1(midIdx, midVals, 1:length(x), 'pchip'); end调用这一段时,maxSweep设为80到120之间。如果你把系数设置得太大,分解出的IMF会呈现明显正弦形态,看起来美观但物理意义变差。
4.3 分解结果怎么看:相关性与瞬时频率
分解完成后,先用plot画出各个IMF波形,再算一下各IMF与已知真实分量之间的相关系数。这个步骤用来判断分解的模态分离度。
[imfs, residual] = esmd_demo(x, 100); comp1 = sin(2*pi*20*t)'; comp2 = (1+0.5*cos(2*pi*5*t)).*sin(2*pi*60*t)'; for k = 1:size(imfs,2) c1 = abs(corrcoef(imfs(:,k), comp1)); c2 = abs(corrcoef(imfs(:,k), comp2)); fprintf('IMF%d 与20Hz分量相关系数: %.3f,与60Hz分量相关系数: %.3f\n', k, c1(1,2), c2(1,2)); end在我跑通的过程中,前两个IMF经常能分别对应20Hz和60Hz分量,相关系数保持在0.9以上,而脉冲分量会出现在后续某一阶IMF中。稳定性比EMD直接分解要好,模态混叠问题也轻得多。
瞬时频率分析可以用MATLAB的hilbert函数。对感兴趣的IMF计算解析信号,然后对相位做差分得到瞬时频率,绘图观察频率是否集中在目标频率附近。这比单纯看波形更能暴露分解中的模态混叠问题——如果某阶IMF的瞬时频率在20Hz和60Hz之间来回跳变,说明两个分量尚未彻底分离。
4.4 结果分析与常见问题
这个Demo跑下来,常见的现象是第一个IMF包含较多噪声成分,出现“噪声模态”的模糊感。这不是ESMD单独的问题,而是所有自适应分解在面对强噪声时的通病。对策是先对原始信号做一次小波阈值消噪或带通滤波,把噪声压低了再送入ESMD,得到的结果通常更干净。
另一个常见现象是,脉冲分量在分解后被拆到多个IMF中。这是因为脉冲的频带很宽,不同频率成分被不同阶数的IMF分别吸收了。如果想完整保留脉冲特征,建议在分解前进行脉冲定位,把脉冲区间单独切出来处理,或者用加窗的方式抑制脉冲扩散。
5. 实际工程中使用ESMD的参数定夺和避坑清单
5.1 不同信号类型下的参数参考
ESMD的“核心参数”其实就两个半:最大筛分次数mirror延拓长度,外加预处理的去趋势策略。参数对结果的影响,往往超过算法本身的选择。我把不同信号类型的参考配置列在下面。
| 信号类型 | 典型场景 | 筛分次数 | 端点延拓 | 预处理建议 |
|---|---|---|---|---|
| 振动信号 | 轴承故障、齿轮箱 | 80~120 | 镜像延拓 | 高通滤波+包络解调 |
| 气象水文 | 风速、径流、气温 | 20~50 | 镜像延拓 | 去除年周期趋势 |
| 生物电信号 | 肌电、脑电、心电 | 60~100 | 反对称延拓 | 工频陷波+带通滤波 |
| 电力信号 | 负荷、电压波动 | 30~60 | 镜像延拓 | 去除直流分量 |
| 声学信号 | 语音、水声 | 50~80 | 线性外推 | 预加重 |
这些数值来自我自己的测试经验,不是数学上的最优解,但作为起始值足够稳妥。拿到一个新的信号类型,先跑一组筛分次数扫描,找到重构误差平缓的区间,再固定参数使用。
5.2 六个常见坑及对策
坑一:极值点检测把平台误判成多个极值点。对策是使用带允许平台逻辑的检测方法,或在极值检测前对信号做三点滑动平均。
坑二:筛分次数设置过大,IMF被“榨干”成纯正弦波。对策是绘制筛分次数与IMF能量曲线,选择拐点处次数,而不是一味求大。
坑三:端点延拓方法照搬,导致低频分量整体偏移。对策是先用FFT估计主周期,取一倍主周期的长度作为延拓尺度。
坑四:分解前不处理趋势项,直流分量混入第一个IMF。对策是先用detrend函数或拟合多项式去除趋势,再进分解流程。
坑五:把残差项丢掉。对策是把残差作为趋势项纳入后续分析,尤其在做时间序列预测时,残差项往往是建模的重要输入。
坑六:盲目认为ESMD一定比VMD或CEEMDAN强。对策是先试用两种算法,用“分解后重构误差+相邻IMF正交性指标”量化对比,用数据而不是直觉做选择。
5.3 ESMD与深度学习的结合思路
最后聊一个更进阶的方向。热词里有大量“bp神经网络拟合曲线”“bilstm代码matlab soc”相关的需求,这说明很多人关注的是预测建模。ESMD在这里可以扮演一个不错的预处理角色:先用ESMD把复杂序列分解成若干IMF和残差,再分别对各分量建模预测,最后叠加输出。这个思路在风速预测、电价预测领域已经有不少实证研究支持,核心逻辑是每个IMF的频带较窄,时序规律相对单一,模型拟合难度大幅降低。
具体在MATLAB中操作时,注意要把ESMD分解得到的IMF按列保存成矩阵,每个IMF单独输入至BP或LSTM模型。测试时用滑动窗口的方式分别预测各分量,最后合成完整预测结果。相比直接对原始序列建模,这种方式在非平稳数据上的预测误差通常能降低一截,但计算量也会同步增加,需要根据业务场景权衡。
我在实际使用中有一点体会比较深:ESMD这类自适应分解方法没有固定的“标准答案”,参数和流程都需要根据你要分析的数据形态反复调试。别怕麻烦,多画图、多对比,时间花下去,对信号本身的理解比跑通一个工具箱更重要。先从一个Demo信号开始,慢慢换成你的真实数据,再逐步微调筛分次数和延拓方式。这篇文章里的代码都不长,建议你亲手敲进MATLAB里逐段执行,趁热把每一步的输出都打出来看一眼,比你背住十篇理论文章都管用。