做信号处理这些年最头疼的事情之一就是面对一堆混叠在一起的序列信号。你明明知道里面既有趋势性的缓变成分也有周期性的中频脉动还有高频噪声但滤波器一刀切下去边界糊成一片相位还给你扭得乱七八糟。后来我开始用EEMD分解来做这件事再搭配样本熵来筛选本征模态函数IMF最后按熵值大小把信号重构为低频、中频、高频三段整个流程才算真正跑通。这个方法尤其适合机械振动信号、供水管网噪声、生理电信号这类非平稳、非线性的序列数据。这篇就把整套思路和实操代码写出来包括我踩过的坑和参数调优经验给正在做信号频带分离的朋友一个可以直接抄作业的参考。1. 整体思路为什么要在EEMD之后还要用样本熵来挑IMF1.1 EEMD解决了EMD的什么问题先说说EMD。经验模态分解会把一个复杂信号按时间尺度逐级拆成若干个IMF每个IMF要求局部对称、过零点数和极值点数最多差一个。听起来很完美实际用起来却有个让人抓狂的毛病模式混叠。假如信号里同时存在一个连续的高频小振幅成分和一个间歇出现的大振幅同频成分EMD会把这俩拆到同一个IMF里或者把一个完整模态撕裂到相邻两个IMF里导致后续分析全是错的。EEMD的核心思路就是给原始信号加多次白噪声利用噪声在集成平均中互相抵消的特性把不同尺度的信号强制分离到不同IMF中去。每次加入不同的噪声序列做一次EMD然后把所有分解结果按IMF序号做平均。这个“加噪-分解-平均”的过程看起来简单但效果立竿见影。我在实际处理供水管网噪声数据时原始EMD的IMF2和IMF3总会出现一段频率跳跃换成EEMD之后IMF边界干净了很多后续做频带划分终于不用再靠肉眼硬猜。1.2 样本熵为什么适合做IMF筛选分解完拿到十几个IMF问题来了哪些是高频噪声哪些是有意义的振荡哪些是趋势项一般做法是看频谱、算能量、算相关系数但这些方法要么需要人为选阈值要么对非平稳信号不敏感。样本熵的优势在于它反映的是时间序列在模式维度上出现新信息的概率值越大说明序列越复杂、越随机值越小说明越有规律、越接近确定性成分。放到IMF筛选这个场景里高频IMF通常是杂乱噪声样本熵很高低频趋势项变化平缓样本熵较低中频有用成分介于两者之间。我试过用方差和能量排序效果远不如样本熵直观。原因很简单方差大不代表复杂度高一个大幅值的正弦波方差可能比小幅值噪声高得多但它的熵并不高。样本熵是从“可预测性”角度去度量信号天然适合判断一个IMF里到底是“有组织”的成分还是“无组织”的残差。1.3 信号重构为低中高频的判定逻辑有了每个IMF的样本熵剩下的工作就是分组重构。我的做法是先计算每个IMF的样本熵值然后按熵值从小到大排序把序列分成三段熵值最小的几个IMF相加作为低频重构信号中间熵值的作为中频重构信号熵值最大的几个作为高频重构信号。这种“按熵排序后再三等分”的方式比单纯设定绝对阈值更稳健因为不同信号的熵值分布范围千差万别没有一套固定的阈值能通吃所有场景。用熵值大小来对应频率高低初看有点反直觉但实际效果很好。因为EEMD分解出的IMF本身就有从高到低的频率分布趋势样本熵和频率有一定的相关性。高频分量随机性强熵高低频分量规律性强熵低。重构出的三频段信号再叠加回去和原始信号的误差能控制在很小的范围内说明这个分组策略没有丢失关键信息。2. 工具选型与数据准备2.1 Python环境与关键库这一整套流程我用Python实现核心库三个numpy负责数组运算scipy负责信号处理和积分PyEMD负责EEMD分解。PyEMD不是标准库需要单独安装pip install EMD-signal就行。另外还要用到matplotlib画图做验证。我建议用Anaconda管理环境Python版本3.9到3.11之间都能跑不要追求最新版PyEMD对新版Python的兼容性偶尔会出问题。我机器上是Python 3.10跑得很顺。scipy版本用1.10以上里面信号处理的接口更稳定。2.2 仿真信号的构造为了验证整个链路是否可靠我习惯先构造一个已知成分的仿真信号。这个信号包含三部分一个2Hz的低频正弦波一个50Hz的中频正弦波还有一个200Hz的高频衰减振荡加随机噪声。采样率设为1000Hz时长1秒一共1000个点。这样设计能让三个频段在频谱图上分得很开方便观察重构效果。import numpy as np fs 1000 # 采样率 1000 Hz t np.linspace(0, 1, fs, endpointFalse) low_freq_sig 1.5 * np.sin(2 * np.pi * 2 * t) # 2 Hz 低频 mid_freq_sig 1.0 * np.sin(2 * np.pi * 50 * t) # 50 Hz 中频 high_freq_sig 0.6 * np.sin(2 * np.pi * 200 * t) * np.exp(-t * 20) # 200 Hz 衰减高频 noise 0.3 * np.random.randn(len(t)) # 随机噪声 signal low_freq_sig mid_freq_sig high_freq_sig noise为什么要加噪声因为实际信号不可能那么干净。EEMD本身就是为了抵抗噪声影响如果仿真信号不带噪EEMD的优势就体现不出来。另外噪声的存在会让样本熵的计算结果更有区分度便于观察不同IMF之间的熵值跨度。2.3 采样率与频谱的关系很多新手会问“不知道采样率怎么求频率频谱”。这个问题必须先搞清楚否则后面做频谱验证全是乱猜。采样率决定了奈奎斯特频率也就是频谱能分析到的最高频率是采样率的一半。如果只知道数据点数和时间长度可以用点数 / 总时长倒推采样率。比如一段信号有1000个点记录了0.5秒那采样率就是2000Hz频谱能分析的最高频率是1000Hz。实际操作里我见过有人做频谱分析忘了带采样率参数结果scipy.signal.periodogram给出的频率轴完全是错的。所以先用fs len(data) / duration确认采样率再往下走。这个坑很小但一旦踩了后面的所有分析全部作废。3. 实操过程从EEMD分解到样本熵计算再到重构3.1 EEMD分解的完整代码PyEMD库的EEMD接口使用起来很直观。关键参数有三个n_std是添加噪声的标准差相对于信号的比例n_splits是集成次数S_number是筛迭代次数。我的初始参数是噪声标准差0.1倍信号标准差集成次数100次实测效果比较稳。from PyEMD import EEMD eemd EEMD() eemd.trials 100 # 集成次数 eemd.noise_width 0.1 # 噪声幅值比例相对信号标准差 imfs eemd.eemd(signal)分解完成后imfs是一个二维数组每一行是一个IMF最后一行通常是残差项。代码里没有直接展示聚类或者分组这些要我们自己算。先看分解出来的IMF数量print(imfs.shape)一般会得到八九个IMF残差是单调趋势。如果数量太少说明EEMD参数没调好噪声加少了或者集成次数不够。3.2 每个IMF的样本熵怎么算样本熵的计算不复杂但很多库没有现成函数需要自己写。核心逻辑是给定时间序列长度为N设置嵌入维数m和相似容限r统计m维向量中任意两向量距离小于r的个数再统计m1维的数量两者比值的负对数就是样本熵。我写了一个简洁版直接用numpy实现适用于中等长度序列。r一般取原始信号标准差的0.2倍m取2这是经过验证的常见配置。def sample_entropy(signal, m2, r0.2): signal np.asarray(signal, dtypenp.float64) N len(signal) r r * np.std(signal) def _max_dist(x, y): return np.max(np.abs(x - y)) # 统计m维匹配对数 def _count_matches(m): count 0 templates np.array([signal[i:im] for i in range(N - m 1)]) for i in range(len(templates)): for j in range(i 1, len(templates)): if _max_dist(templates[i], templates[j]) r: count 1 return count A _count_matches(m 1) B _count_matches(m) if A 0 or B 0: return 0.0 return -np.log(A / B)这个实现是O(N²)复杂度序列特别长时会很慢。实际处理振动信号时每个IMF通常有几万个点我一般先用滑窗平均降采样到2000点以内再算熵结果影响很小。对速度有要求的话可以用KD树加速但这不是本文重点先跑通再说。计算每个IMF的样本熵entropies [] for imf in imfs: entropies.append(sample_entropy(imf))3.3 按样本熵排序并重构三频段信号拿到所有IMF的熵值后我习惯打出来看一眼。正常情况下高频IMF的熵值会明显高于低频IMF。如果出现中间某阶IMF的熵值异常说明EEMD分解可能有混叠这时需要调整参数重新分解。分组的代码逻辑如下order np.argsort(entropies) # 从小到大排列的IMF索引 n_imfs len(imfs) # 三等分 low_idx order[:n_imfs // 3] mid_idx order[n_imfs // 3 : 2 * n_imfs // 3] high_idx order[2 * n_imfs // 3:] low_signal np.sum(imfs[low_idx], axis0) mid_signal np.sum(imfs[mid_idx], axis0) high_signal np.sum(imfs[high_idx], axis0)这里有个细节如果IMF数量不是3的整数倍需要手动处理边界。我一般用np.array_split对排序后的IMF索引做切分它会自动分配多出来的几个给前面的组避免分组不均。groups np.array_split(order, 3) low_signal np.sum(imfs[groups[0]], axis0) mid_signal np.sum(imfs[groups[1]], axis0) high_signal np.sum(imfs[groups[2]], axis0)我试过按熵值绝对阈值切分比如熵值小于0.5归低频但不同批数据的熵值分布很不固定阈值难调。排序后三等分最大程度保留了信号本身的能量分布适应性更强。3.4 用频谱验证重构效果重构完了不能直接扔上去必须做验证。最直观的验证方式就是频谱图。对三个重构信号分别做FFT看它们的主频是否落在预期频段。from scipy.signal import periodogram f_low, p_low periodogram(low_signal, fs) f_mid, p_mid periodogram(mid_signal, fs) f_high, p_high periodogram(high_signal, fs)画图时我会把三个频谱放在同一个图里用不同颜色区分观察峰值频率是否分别集中在低频段、中频段、高频段。如果是说明分组成功如果低频重构信号里出现了高频尖峰说明该IMF被错误分到了低频组需要检查样本熵有没有计算错。这里顺带提一下“分段频谱图”的做法。有时原始信号长达几分钟直接整段FFT会把不同时段的频率特征平均掉。我会把信号按固定时窗切成段对每段做FFT再把频谱按时间拼接成二维图。这个方法在供水管网噪声记录仪的频带划分上特别有用因为管网泄漏噪声在某个频带上的能量变化是随时间波动的整段频谱看不出特征分段频谱图能呈现动态变化过程。4. 踩坑记录与参数调优4.1 EEMD参数噪声幅值与集成次数怎么设这是EEMD里影响最大的两个参数。噪声幅值noise_width太小抑制模式混叠的作用就不明显太大又会产生虚假IMF甚至把真实信号淹没。我做过一组对照实验同一段信号noise_width从0.01逐步提到0.5当超过0.3时分解结果的残差项出现明显振荡说明过拟合了。个人经验是0.1到0.15之间最稳尤其是信噪比不高的实测信号。集成次数trials理论上越大越好但计算时间成倍增加。100次大概是精度和速度的平衡点。如果信号本身比较干净可以降到50次如果信号里有强烈的间歇性成分建议加到200次。我处理供水管网脉冲噪声时用到过300次效果确实更好但单个文件要跑十几分钟前期调试时很不划算。所以建议先用100次跑通最后批量处理时再按需加次数。4.2 样本熵参数m、r怎么选择样本熵的m一般取2这是时间序列复杂度分析的标准选择。m1时抗噪性差m3时对数据长度要求变高。r的选取更微妙它决定了两个序列段算不算“相似”。r太小几乎没有匹配对熵值趋近无穷大或无定义r太大所有序列段都相似熵值趋近0区分度消失。我测试过不同r值对IMF排序结果的影响最后固定在0.2倍信号标准差。这个值在大多数情况下能拉开不同IMF的熵值差距。有个小技巧如果发现所有熵值都特别接近先看是不是r取大了如果只有个别IMF熵值是0或无穷大说明数据长度不够或者序列太规则。序列长度最好大于500个点少于200个点算出来的熵基本没有参考价值。4.3 真实采集信号的注意事项真实采集信号和仿真信号完全是两回事。首先EEMD要求输入是等间隔采样的一维序列如果采集设备有丢包或者时间戳抖动必须先重采样到均匀时间轴。我遇到过供水管网记录仪时间戳偶尔跳变导致EEMD分解出的IMF前几阶全是毛刺最后用三次样条插值重采解决了。其次真实信号往往有明显趋势和直流分量。EEMD的第一阶IMF可能包含很大的趋势项直接计算样本熵会把整个分组带偏。我建议先做去趋势预处理signal signal - np.mean(signal)再用高通滤波器滤掉0.5Hz以下成分。低频趋势对样本熵影响极大因为趋势项会让序列看起来很有规律熵值虚低干扰分组。还有一个容易忽略的点EEMD分解得到的IMF数量越多计算样本熵的时间越长。如果信号长达几十万点建议先做分段处理或者降采样。比如100kHz采样率的机械振动信号我通常会降到10kHz再做分析保留主要频带信息速度提升明显。5. 嵌入式移植的可行性参考5.1 STM32F4做FFT与EEMD的差距最近热搜里能看到“基于STM32F4的嵌入式FFT频谱分析系统设计”很多人想把信号频带分离这套算法嫁接到单片机上。我明确说EEMD加样本熵整套流程放在STM32F4上实时跑现实意义不大。STM32F4主频最高168MHz做1024点FFT是绰绰有余但EEMD需要多次迭代每次迭代都要做一次EMD筛再加上100次集成平均运算量是FFT的几千倍实时计算根本不现实。更合理的嵌入式方案是在STM32F4上只做数据采集和FFT把频谱数据通过串口或CAN发给上位机由上位机完成EEMD分解和样本熵计算。如果一定要在MCU端做频带重构建议用“分段频谱图”的思路代替EEMD先对信号分段做FFT再按频带划分能量实现低复杂度的高中低频带分离。这种做法虽然不如EEMD自适应但在资源受限的平台上更加可行。5.2 分段频谱图思路在重构信号中的应用分段频谱图可以视为一种“面向工程实现”的频带划分手段。把信号每256点或512点一段对每段做FFT得到频谱序列。然后按照目标频带划分规则比如低频10Hz到100Hz、中频100Hz到1kHz、高频1kHz到5kHz分别累加各段频谱幅度得到每段信号的频带能量。再对这些能量序列做平滑就能绘制出“频率-时间-能量”二维图。我把这个方法用在供水管网噪声记录仪上用来判断管道泄漏的频带特征。那台记录仪就是基于STM32F4的方案采样率20kHz每段512点FFT基本能实时输出9个频带的能量分布。后来在和EEMD重构结果对比时发现如果只关心低频、中频、高频三个趋势分段FFT的频带能量累积曲线和EEMD重构后的信号包络非常接近。区别在于EEMD得到的三个重构信号可以直接做时域分析比如算有效值、峰峰值、互相关系数而分段FFT只能给频带能量。5.3 不知道采样率时怎么反推频率最后补一个实操细节也是很多人问过的采集到的数据文件里没有保存采样率怎么求频谱频率最原始的思路是用设备固件里的采样配置反推但一旦配置文件丢失就抓瞎了。我常用的办法是找信号里的已知特征频率。比如电网相关的信号里肯定有50Hz工频机械振动信号里可以通过转轴转速计算转频然后以这个已知频率为基准去校正频谱坐标。具体做法是先对信号做FFT找到工频峰值对应的索引k那么真实频率对应的分辨率是真实频率 / k也就是df 50 / k采样率就等于df * N。这个方法在信噪比尚可的情况下误差很小。如果信号里连已知频率都没有就只能靠数据记录时间戳换算。整个文件从首尾时间戳算出总时长T点数N已知平均采样率就是N / T。需要注意的是这个平均采样率只能用于整段FFT频谱分辨率是1/T如果中间有重采样或者时间戳跳变这个方法也会失真。我个人做信号频带分离这几年最深的体会是算法本身都不复杂难的是每个环节的参数适配和异常判断。EEMD加样本熵这套组合胜在自适应性强不需要频繁调整滤波器系数但你必须对每一步的原理有清晰认知才敢在参数异常时做判断。我自己的经验是噪声幅值设在0.1倍信号标准差、集成次数100次、样本熵的容限取0.2倍信号标准差这套配置在振动、噪声、生理信号上表现都很稳定。最后再分享一个小技巧如果只想分低中高三个频段不一定要硬性按样本熵的绝对阈值去卡可以先对熵值排序再做三等分这样分组结果更贴合信号自身的复杂度和能量分布而且几乎不用调参。后续如果遇到频段边界不清晰的情况可以增加分组数比如分成五段再观察每段重构信号的频谱形态找到最适合你研究对象的划分粒度。