EEMD分解结合样本熵的IMF筛选与信号低中高频重构实战
2026/9/14 4:52:56 网站建设 项目流程

做信号处理这些年,最头疼的事情之一就是面对一堆混叠在一起的序列信号。你明明知道里面既有趋势性的缓变成分,也有周期性的中频脉动,还有高频噪声,但滤波器一刀切下去,边界糊成一片,相位还给你扭得乱七八糟。后来我开始用EEMD分解来做这件事,再搭配样本熵来筛选本征模态函数IMF,最后按熵值大小把信号重构为低频、中频、高频三段,整个流程才算真正跑通。这个方法尤其适合机械振动信号、供水管网噪声、生理电信号这类非平稳、非线性的序列数据。这篇就把整套思路和实操代码写出来,包括我踩过的坑和参数调优经验,给正在做信号频带分离的朋友一个可以直接抄作业的参考。

1. 整体思路:为什么要在EEMD之后还要用样本熵来挑IMF

1.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, endpoint=False) 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的个数,再统计m+1维的数量,两者比值的负对数就是样本熵。

我写了一个简洁版,直接用numpy实现,适用于中等长度序列。r一般取原始信号标准差的0.2倍,m取2,这是经过验证的常见配置。

def sample_entropy(signal, m=2, r=0.2): signal = np.asarray(signal, dtype=np.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:i+m] 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], axis=0) mid_signal = np.sum(imfs[mid_idx], axis=0) high_signal = np.sum(imfs[high_idx], axis=0)

这里有个细节:如果IMF数量不是3的整数倍,需要手动处理边界。我一般用np.array_split对排序后的IMF索引做切分,它会自动分配多出来的几个给前面的组,避免分组不均。

groups = np.array_split(order, 3) low_signal = np.sum(imfs[groups[0]], axis=0) mid_signal = np.sum(imfs[groups[1]], axis=0) high_signal = np.sum(imfs[groups[2]], axis=0)

我试过按熵值绝对阈值切分(比如熵值小于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,这是时间序列复杂度分析的标准选择。m=1时抗噪性差,m=3时对数据长度要求变高。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倍信号标准差,这套配置在振动、噪声、生理信号上表现都很稳定。最后再分享一个小技巧:如果只想分低中高三个频段,不一定要硬性按样本熵的绝对阈值去卡,可以先对熵值排序再做三等分,这样分组结果更贴合信号自身的复杂度和能量分布,而且几乎不用调参。后续如果遇到频段边界不清晰的情况,可以增加分组数,比如分成五段,再观察每段重构信号的频谱形态,找到最适合你研究对象的划分粒度。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询