做信号处理这些年我有个越来越强烈的感受傅里叶变换就像一副度数固定的眼镜戴上它看平稳信号很清晰一碰到非平稳、非线性信号就“花眼”。后来接触了经验模态分解EMD和边际谱计算这条技术路线算是把这副眼镜换成了自适应变焦的。这篇博文就围绕“从EMD到边际谱”这条主轴把原理、算法细节、实际跑通的完整流程、还有我踩过的坑一次讲清楚。适合刚入门时频分析的研究生、做故障诊断的工程师以及任何被非平稳信号折腾过的朋友。先说结论EMD能把任意复杂信号自适应地拆成若干本征模态函数IMF加上一个残差每个IMF再做Hilbert变换得到瞬时频率和瞬时幅值汇总成Hilbert谱把Hilbert谱对时间求积分得到的就是边际谱。边际谱的横轴是频率纵轴是该频率在所有时刻上累积的能量或幅值贡献它比傅里叶谱更能揭示非平稳信号里“真实存在过”的频率成分。下面我会从“为什么非它不可”讲到“代码怎么写、参数怎么调、坑怎么避”按照我自己的理解顺序来。1. 为什么不用傅里叶传统时频分析的两道坎1.1 傅里叶谱是“全局平均”掩盖了非平稳的本质多数人学完信号处理第一反应就是把信号扔给FFT看一眼频谱图。这没错但傅里叶变换有一个隐含假设信号是由若干固定频率、固定幅值的正弦波叠加而成的而且这些正弦波从负无穷到正无穷一直存在。对一段10秒的振动信号FFT算出来的谱峰只能说“这10秒里大概有这个频率成分”它完全丢掉了“这个频率是什么时候出现的、幅值怎么变化”的信息。举个实际例子某旋转机械在启动过程中转速从500转/分升到3000转/分振动信号里的转频分量从8Hz一路爬到50Hz。你用FFT去看这整段启动信号只能看到一段“糊在一起”的宽谱峰值落在某个折中频率上。你要是用这段FFT结果去做故障判断可能连“机器是匀速运行还是变速运行”都看不出来。问题的根源就是傅里叶谱做了时间平均非平稳信号的特征恰恰在“随时间改变”这一点上平均一搞特征就没了。1.2 短时傅里叶与小波不是不行是“窗口”太难伺候为了补救工程上常用短时傅里叶变换STFT和小波变换。STFT的思路是把信号切成一小段一小段对每段做FFT然后拼成一张时间-频率图。它最大的问题是窗口长度固定窗口短时间分辨率高但频率分辨率差窗口长频率分辨率好了但时间精度又丢了。“海森堡测不准原理”摆在那STFT就是在墙上的两边来回撞。小波变换稍微聪明一点它用可伸缩的母小波去匹配信号低频用长窗口、高频用短窗口实现了多分辨率分析。但小波也有个绕不开的问题需要预先选母小波还要定分解层数。选不同母小波结果差异之大足以让人怀疑自己是不是在做同一个实验。更麻烦的是小波基一旦选定分解方式对整段数据是统一的它仍然是“用固定基函数去拟合任意信号”的思路缺少自适应性。你说信号里有个突变成分它不属于任何光滑的小波基小波分解就会把它拆得稀碎散落到多个尺度上物理意义跟着变模糊。所以当信号是非平稳、非线性、且成分复杂到没法预先假设基函数时我宁愿选一条“让数据自己说话”的路子。这就是我转向经验模态分解的原因。2. EMD分解原理自适应的“拆信号”艺术2.1 IMF的数学定义与物理直觉经验模态分解的核心产物叫本征模态函数IMF。它的定义有两个条件看着简单但理解透了整个算法就好懂了整个数据序列里极值点数目与过零点数目相等或至多相差1个。在任意时间点上由局部极大值定义的上包络与局部极小值定义的下包络的均值为0。换句话说IMF要求信号局部对称上下包络关于零线基本镜像而且振荡形态“像”一个调幅调频的正弦波但不能出现“骑在趋势线上”的情况。这跟傅里叶里“固有模态”概念不太一样EMD不是用固定正弦基去分解而是通过信号自身的极值分布来定义“模态”。所以它分解出来的IMF通常具有物理意义比如第一个IMF往往是高频细节后面IMF频率逐步降低最后一个残差是趋势项。我用一个生活化的类比来理解假设你手里一堆混在一起的布料不同颜色的线粗细、长短都不一样。FFT的做法是你准备不同粗细的固定筛子强行按口径分EMD的做法是先用手指感知每根线的走向逐根抽出来抽出来之后还要确认它是不是“纯色纯粗细”的一根线如果不是就继续抽。抽到最后剩下的那团底线就是趋势项。2.2 筛分算法包络、均值、迭代EMD的筛分过程sifting process是核心循环我来拆解它每一步在干什么找到原始信号 ( x(t) ) 的所有局部极大值点和局部极小值点。用三次样条插值分别对极大值点和极小值点拟合得到上包络 ( e_{up}(t) ) 和下包络 ( e_{low}(t) )。计算上下包络的均值 ( m(t) (e_{up}(t) e_{low}(t)) / 2 )。从原信号里减去均值( h(t) x(t) - m(t) )。检查 ( h(t) ) 是否满足IMF的两个条件。不满足就把 ( h(t) ) 当作新的信号回到第1步继续筛满足就记为一个IMF ( c_1(t) )。用残差 ( r_1(t) x(t) - c_1(t) ) 重复以上过程得到 ( c_2(t) )、( c_3(t) )……直到残差是单调函数或幅度小于预设阈值为止。这里最容易让人疑惑的是“为什么要反复筛”而不是一次减去包络均值就结束。因为第一次做出来的 ( h(t) ) 往往还包含非对称的波动包络均值不等于0只有多次迭代才能让IMF逐步“纯化”局部均值趋于零。我试过只筛一两次就跑出来的IMF波形明显不对称做Hilbert变换后瞬时频率跳得厉害边际谱上一堆虚假峰值。所以筛分迭代别偷懒它是EMD质量的生命线。不过三次样条插值也有一个副作用——过冲。就是样条在极值点之间穿行时可能甩出比原始数据更高或更低的弧线引入虚假极值。这个在端点附近尤其明显后面我专门讲端点效应时再说。2.3 停止准则筛到什么时候算完筛分不可能无限迭代下去否则信号会被磨成纯正弦波失去物理意义。常用的停止准则有两类。一类是仿照Cauchy收敛定义标准差 ( SD \sum_{t} |h_{k-1}(t) - h_k(t)|^2 / \sum_{t} h_{k-1}^2(t) )当它小于某个阈值比如0.2到0.3就停止。这个准则关注的是相邻两次迭代结果变化是否足够小。另一类更常用是“S数”准则连续 S 次迭代中过零点数和极值点数保持相同就停止。我偏好S数准则因为IMF的物理形态极值-过零关系稳定了说明筛分已经收敛到稳定的振荡结构。通常S取3到10取3时分解结果更“激进”IMF数量多取10时更保守但可能过度光滑。实际做的时候我会先把S设成3跑一遍看IMF形态再调大跑第二遍对比。停止准则没有绝对“正确”值它取决于你之后要干什么。如果目标是做边际谱看频率成分稍微激进一点即S偏小通常问题不大如果目标是提取IMF做进一步时频分析平稳性更重要S取大一些。这里必须提醒一点EMD分解是逐级剥离高频分量但并非“先剥离的信号一定频率最高”。它剥离的是“包络均值最大的振荡结构”这个结构一般对应高频遇到模态混叠时情况就不纯粹了。这一点到第3节讲边际谱时还得往回看。3. 从IMF到边际谱Hilbert变换这一跳3.1 瞬时频率为什么非用Hilbert变换不可拿到IMF之后下一步要算“瞬时频率”。瞬时频率的定义有争议但从工程角度最实用的解析途径是Hilbert变换。对实信号 ( c(t) )它的Hilbert变换 ( \hat{c}(t) ) 是[ \hat{c}(t) \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{c(\tau)}{t - \tau} d\tau ]这个积分在 ( t \tau ) 处是奇异积分意思是取柯西主值。把实信号和它的Hilbert变换组合成解析信号[ z(t) c(t) j \hat{c}(t) a(t) e^{j \theta(t)} ]其中瞬时幅值 ( a(t) \sqrt{c^2(t) \hat{c}^2(t)} )瞬时相位 ( \theta(t) \arctan (\hat{c}(t)/c(t)) )瞬时频率 ( f(t) \frac{1}{2\pi} \frac{d\theta(t)}{dt} )。为什么要对IMF做Hilbert变换而不是直接对原信号做因为瞬时频率的定义在任意实信号上是“不唯一”且缺少物理意义的一个多分量信号算出来的瞬时频率是各个分量频率的加权乱七八糟混在一起没法解释。IMF是单分量信号近似瞬时频率才稳定、才有意义。所以EMD是Hilbert谱分析的前置条件得先分解再变换顺序反了全是噪声。我在最初学习时犯过一个错误直接对原始振动信号做Hilbert变换求瞬时频率结果频率曲线锯齿状乱跳完全看不出转频爬升的轨迹。后来想明白了Hilbert变换只对“局部窄带”信号友好宽带多分量信号一进去瞬时相位剧烈扭转求导后频率自然乱飞。这就是为什么Hilbert-Huang变换的核心是“先分解、再变换”。3.2 从Hilbert谱到边际谱时间维度上的累积对每个IMF都求出 ( a_i(t) ) 和 ( f_i(t) )然后把它们画在同一张时间-频率平面上就得到Hilbert谱 ( H(f, t) )它反映的是幅值随时间、频率的分布是一张三维图像。Hilbert谱可以看作是信号能量在时频平面上的“热力图”每一个时间点上的瞬时频率和瞬时幅值都被保留。边际谱是Hilbert谱对时间的积分[ h(f) \int_{0}^{T} H(f, t) , dt ]一句话理解把每个频率值在所有时刻的“存在感”累计起来。如果某个频率在信号生命周期里出现过较长时间或者幅值比较大它在边际谱上的值就高。反过来如果一个频率只是瞬间闪现在边际谱上只留下一个小峰。这正是边际谱跟傅里叶谱最大的区别。对比起来特别直观对同一段变速升频信号FFT的谱峰是一个宽缓的包把一个扫频过程“平均”成一个不存在的峰值边际谱则会在扫频范围内铺开一条接近矩形的谱带精确告诉你“这个信号从低到高扫过了哪些频率”能量分布范围一目了然。我在给某设备做启停机分析时用边际谱能清楚看到转频从8Hz爬到50Hz的过程FFT只能给出一坨模糊的宽峰谁更有利于诊断不用多说。3.3 边际谱的处理细节去端点、去筛尾边际谱虽然好用但直接从程序里“跑”出来的结果不能直接用有两个后处理要点第一瞬时频率在信号两端会因为Hilbert变换的边界效应产生异常值频率可能冲到负频或远高于奈奎斯特频率。在积分到边际谱之前我通常会把每个IMF在时间首尾砍掉约2到3个振荡周期的片段或者对这些边缘区域做加权平滑再计算瞬时频率。第二IMF末尾几个分量往往是低频残差趋势它们对应的瞬时频率极低在边际谱上会凝聚成一个覆盖超低频的巨大隆起淹没中高频信息。我一般做边际谱时只取前几个能量占比高的IMF残差直接忽略或者对低频段单独展示。做故障诊断时我经常只看前3到5个IMF的边际谱后面那些低频分量更像是数据趋势跟故障特征关联不强。还有一个细节边际谱的纵轴单位取决于 ( a_i(t) ) 的定义。如果你用的是瞬时幅值直接平方那就是能量谱用幅值本身则更接近“贡献幅值”。我习惯用幅值平方再归一化跟FFT功率谱对比时量纲更接近。4. 实操过程用Python跑通EMD边际谱全流程4.1 环境准备与数据构造我日常用的是Python核心库是信号处理库中的EMD模块另外还要用到傅里叶变换、数值积分和绘图的几个常规库。安装一句话的事不用我啰嗦。如果你是MATLAB阵营的也有现成的EMD和HHT函数流程基本一致只是接口名不同。这里构造一段模拟信号来演示它来自我做诊断实验时常用来验证算法的合成信号一个频率随时间线性升高的调频分量叠加上50Hz的固定高频分量再加一点随机噪声。为什么用这个信号因为它包含了非平稳特性变频分量和平稳特性定频分量非常适合看出EMD和边际谱的区分能力。import numpy as np from scipy.signal import hilbert from emd import sift fs 1000 t np.linspace(0, 10, 10*fs) f_scan 5 30*t/10 x 0.8*np.sin(2*np.pi*f_scan*t) x 0.4*np.sin(2*np.pi*50*t) x 0.05*np.random.randn(len(t))采样率选1000Hz是因为最高分析频率50Hz按照奈奎斯特条件1000Hz留了十倍余量瞬时频率的跟踪会更平滑。如果你实际信号的最高频率接近采样率一半建议把采样率提到至少5倍以上目标频率不然瞬时频率的求导噪声会很大。4.2 分解操作与参数选择调用EMD接口做分解imf sift.sift(x)sift函数默认会返回一个二维数组每一行是一个IMF最后一行是残差。默认参数对多数测试信号已经够用不过有几个参数我会根据信号特点调整max_imfs限制IMF的最大数量防止分解出太多低频小分量。做诊断时我一般限制到6到8个避免分析重点被琐碎分量带偏。sd_thresh停止准则阈值默认值偏严格可以适当放一点减少迭代次数。envelope_style包络样条类型默认是经典三次样条。如果信号含较多噪声我会换成带平滑的方式但代价是可能损失高频细节。做完之后先别急着算边际谱一定要看一遍分解结果。我通常把IMF波形和原始信号叠着画出来检查第一个IMF是不是捕获了最高频的振荡结构后几个IMF频率是否逐步走低有没有某个IMF出现明显的“双峰包络”或不对称形态如果有说明该分量可能混叠了两个源成分。4.3 Hilbert谱与边际谱的绘制计算每个IMF的瞬时幅值和瞬时频率imf imf[:-1] # 去掉残差或保留但单独处理 n_imf imf.shape[0] inst_amp np.zeros_like(imf) inst_freq np.zeros_like(imf) for i in range(n_imf): analytic hilbert(imf[i]) inst_amp[i] np.abs(analytic) inst_freq[i] np.gradient(np.unwrap(np.angle(analytic)), t) / (2*np.pi)这里用np.unwrap是必须的相位一旦包裹回落到[-π, π]区间直接求导会出现剧烈的负频率尖峰。unwrap把相位展开成连续曲线求导才稳定。拿到瞬时频率和幅值后绘制Hilbert谱把时频平面划分成网格把每个时间点的幅值“散落”到对应网格上得到二维分布图。边际谱则是把瞬时幅值平方在时间轴上按频率桶累加freq_bins np.linspace(0, 80, 400) margin_spectrum np.zeros(len(freq_bins)-1) for i in range(n_imf): bin_idx np.clip(np.searchsorted(freq_bins, inst_freq[i])-1, 0, len(freq_bins)-2) margin_spectrum np.bincount(bin_idx, weightsinst_amp[i]**2, minlengthlen(freq_bins)-1)这个累加过程就是数值化的积分。理论上积分是对连续时间做积分用等时间采样后变成对幅值平方求时间上的和。采样率均匀时这种累加等价于乘以采样间隔再求和。把边际谱画出来你会看到一个宽宽的峰位于5到35Hz对应扫频分量一个尖峰在50Hz对应固定频率分量。跟FFT谱对比着看FFT的扫频分量呈现为15到25Hz之间的一个鼓包边际谱则清楚标出5Hz到35Hz的完整扫描范围。4.4 结果解读的三个观察视角解读边际谱时我习惯按三个视角来看第一是频率定位。峰值数量和中心频率对应哪些可能存在振源。固定频率尖峰往往对应工频或结构固有频率扫频宽带往往是变速过程的特征比如转子转速爬升。第二是能量分布。边际谱里每个峰的“面积”反映了该频带在观测时段内累积的能量。两个峰值面积比可以作为不同振动贡献的定量参考。比如上面合成信号里0.8^2和0.4^2的比值大致是4比1你从边际谱峰值面积比里能看回来。第三是频谱形态。正常轴承或齿轮故障边际谱上高频范围会出现多个有规律的峰群早期磨损则往往表现为某个峰值旁的边带频率逐渐隆起。形态比数值更敏感所以我建议每次分析都先看谱形再看峰值。5. 常见问题与排查心得5.1 模态混叠近频成分“拧”在一起模态混叠是我在EMD实操里遇到最多的故障现象两个频率相近或者特征时间尺度接近的分量被分解到同一个IMF里表现为IMF包络剧烈震荡波形像调幅波。比如采样信号里同时有49Hz和51Hz两个分量理论上原始信号就是拍频形态EMD有可能把其中一个分量筛掉或者一个IMF里同时含两成分。模态混叠的根源在于EMD按“局部极值时间尺度”来切分信号当两个分量的极值点交替出现时包络均值就无法把它们分开。应对方法有几个先做带通滤波把待分析频带限制住再分解。这个听起来“不够自适应”但非常实用。我处理轴承外圈故障信号时会先划分低、中、高频三个子带分别EMD比全局分解干净得多。用集合经验模态分解EEMD或完整集合经验模态分解CEEMDAN处理。它们通过在信号里反复叠加有限幅值白噪声并多次分解取平均打散尺度接近的成分让每个IMF更纯净。代价是计算量大增取平均次数通常上百次。提高采样率和数据长度。数据够长极值点间隔的区分度更高混叠概率下降。模态混叠最直观的判定方法是做Hilbert谱观察如果某个IMF的瞬时频率在某一时间段内出现两个明显台阶八成是混叠了。这时别勉强用这个IMF换EEMD或滤波再跑一次。5.2 端点效应边界会“生长”出假振荡EMD筛分和Hilbert变换都依赖端点附近的极值插值信号开头和结尾的极值信息不完整三次样条外推往往大幅过冲产生“多余的”包络起伏。这种端点效应会沿筛分过程向内传播污染前几个IMF的波形和瞬时频率。我在工程应用里的处理习惯是数据采集时尽可能多采一段分析时把前后各10%到20%的长度切除只取中段稳定区。如果信号固定长度没法切使用特征波延拓或镜像延拓把端点附近的波形镜像翻转拼接让样条外推有据可依。多数EMD库提供了边界延拓参数比如镜像延拓是我的首选。对瞬时频率做同样的端点裁剪算边际谱时直接从有效时间段累加。从结果上看端点修正前后的边际谱在两端频率处差异很大。如果发现边际谱在0Hz附近有个意料之外的锐峰我会第一时间怀疑端点效应而不是物理信号真有一个极低频成分。5.3 停止准则与过筛筛过了会让“真信号”失真前面提到停止准则有SD准则和S数准则我实际用得最多的是S数准则。但S取太小会欠筛IMF里还残留其他成分S取太大容易过筛把一个真实信号中的突变细节磨平。我有个不太好听但很实际的体会EMD的结果对停止准则并不是特别鲁棒尤其在数据短、噪声大的情况下。不同停止准则跑出的边际谱峰值位置可能偏移几个赫兹。所以我的建议是不要只跑一版就下结论。把S从3到8各跑一遍看哪些谱峰在参数变化下仍然稳定存在那些“稳如泰山”的峰才是可信的物理特征对只出现一次的边缘峰宁可忽略。这个习惯帮我在一次某设备检测中免于误判当时用默认参数跑出一个24Hz的谱峰听起来像外圈故障特征。我换参数再跑24Hz消失同时出现了稳定的47Hz峰。后来拆机确认问题出在与47Hz对应的部件上。参数敏感性分析是真的能救命的。5.4 负频率、瞬时频率突变与不满足单分量用Hilbert变换求瞬时频率时最常见的异常就是负频率和频率突变。负频率的直接原因是相位出现非单调变化比如一个IMF波形有小的扰动导致相位曲线局部回绕unwrap之后出现局部下降段。瞬时频率突变一般发生在信号过零点密度异常、或者幅值接近零处此时相位对时间求导是“0/0型”数值噪声极大。排查思路先用瞬时幅值做掩膜瞬时幅值小于某阈值比如全局最大幅值的5%的时间点其瞬时频率置为无效不参与边际谱累积。这一招能压掉大量虚假谱线。对瞬时频率做中值滤波或Savitzky-Golay平滑后再进入边际谱计算但注意别过度平滑否则扫频分量的带宽被压扁。检查IMF是否真的满足单分量性质如果瞬时频率波动范围超过平均频率的30%建议回到分解阶段重新处理而不是在频率域硬扛。我也遇到过一种情况瞬时频率总在目标频率附近周期性往返像正弦波一样来回扫。后来发现那是把一个分量分解成了相邻两个IMF它们的频率交汇处产生了模态混叠式的边缘。把这两个IMF的瞬时幅值加起来作为合成信号重新算瞬时频率反而得到平滑跟踪。这个“先合再算”的小技巧在调制信号分析里偶尔会用到。为了让你排查起来顺手我把常见问题、症状和处置方式整理成一个速查表常见问题典型症状首要排查方向常用对策模态混叠IMF包络剧烈起伏、Hilbert谱中出现频率台阶原始信号频带是否过宽预滤波分频带处理、改用EEMD/CEEMDAN端点效应边际谱0Hz附近异常尖峰、IMF两端波形大幅弯曲数据长度与边界延拓方式镜像延拓、切除首尾段过筛/欠筛边际谱峰值随停止准则漂移停止准则参数敏感性多参数对比只取稳定峰负频率瞬时频率出现低于0的值相位包裹与幅值接近零点相位unwrap、幅值掩膜虚假谱峰边际谱出现无物理意义的孤立峰噪声水平与筛分迭代次数幅值加权、平滑瞬时频率结尾我从第一次跑通EMD到今天最大的感受是这条路线的核心价值不在“算法多炫”而在于它让非平稳、非线性信号终于有了一个不靠预设基函数也能拆解的路径。但它的代价也很实在——参数敏感、计算慢、需要人反复看结果。所以我给后来者的建议很朴素不要一上来追求完美参数先构造一段已知成分的合成信号把流程跑通用对照实验找手感再碰真实数据。我至今记得第一次用边际谱在扫频信号里看到那条干净频率带时的心情——原来那个一直躲在FFT模糊阴影里的成分是可以这样被看见的。后续如果还有精力我会接着写EEMD和CEEMDAN的实际对比以及边际谱在轴承故障诊断里怎么跟特征频率对上号。这条路还有不少坑但值得一个一个踩过来。