简介基于经验模态分解法EMD的Hilbert-Huang变换HHTMATLAB程序包面向信号处理、故障诊断及振动分析领域的研究人员和工程师用于将非平稳信号转化为平稳信号进而完成时频分析。程序通过EMD将信号分解为若干本征模态函数IMF分量再对IMF分量进行Hilbert变换以获取瞬时频率与幅值最终累加重构得平稳化结果核心逻辑清晰层次直观便于理解HHT算法链路并做二次开发。压缩包内共3个文件均为.m源码包含主程序HHT.m以及频谱计算、瞬时频率估计所需的hhspectrum.m与instfreq.m函数整体仅2KB体量极小运行轻便需配合已安装的EMD工具箱中的emd函数使用。目前已有1277人学习/浏览适合正在学习HHT原理、搭建MATLAB分析流程或复现场景的读者。借助这份代码可以快速掌握IMF分量提取、Hilbert谱构建、瞬时频率计算等关键步骤并直接嵌入自己的实验或项目中减少从零编写的时间成本。1. EMD与HHT从非平稳信号里拆出可解释的时频结构轴承振动、脑电波、风速序列、结构应变这些信号有个共同点频率成分随时间漂移。拿一段转速爬升过程中的振动数据做FFT频谱上能量铺成一片峰值对应的频率说不清是哪个时刻的加窗做短时傅里叶变换又要在时间分辨率和频率分辨率之间做取舍窗宽一固定整个时频面只能用同一把尺子量。经验模态分解Empirical Mode Decomposition, EMD搭配Hilbert-Huang变换HHT走的是另一条路先按信号自身的极值分布把数据分解成若干本征模态函数Intrinsic Mode Function, IMF再对每个IMF做Hilbert变换求出有物理意义的瞬时频率从而得到一条沿时间轴变化的频率曲线。整个过程不预设基函数不要求信号平稳也不需要选择窗函数这两点让它和FFT、小波变换在思路上彻底分开。这套方法适合做振动故障诊断、结构健康监测、生物医学信号分析的技术人员也适合所有对「频率随时间怎么变」有量化需求的场景。下面从EMD的筛分逻辑讲起一路落到能跑的代码和参数调优。2. EMD核心流程从原始信号到IMF的筛分步骤与停止准则2.1 IMF的两个定义性条件过零点约束与包络对称把一个信号拆成多个IMF前提是得先定义清楚「拆出来的每一层到底长什么样」。IMF必须同时满足两个条件一是整个数据段内过零点的数量和极值点的数量相等或至多相差一个二是在任意时间点上由局部极大值拟合出的上包络和由局部极小值拟合出的下包络的均值必须为零也就是上下包络关于时间轴对称。第一条保证了IMF是窄带振荡第二条保证瞬时频率不会因为包络不对称而出现无意义的跳动。直观理解就是IMF不能像原始信号那样既包含慢变的大幅趋势又叠着快速的小幅抖动它必须是一段「围绕零点、上下包络均衡」的振荡波形。这里可以对比一下傅里叶分解FFT把信号投影到无限长的正弦波基上任何一个局部扰动都会摊到整个频谱上而EMD的IMF是从数据本身长出来的极值点密集的地方就分解出高频IMF极值稀疏的地方就落到低频IMF基函数跟着数据走这是「经验」二字的由来。2.2 筛分迭代极值定位、三次样条包络与均值剥离2.2.1 单轮筛分的三个动作EMD对信号的处理不是一次性分解而是反复迭代。每一轮筛分做三件事找出当前信号的所有局部极大值和局部极小值点用三次样条曲线分别拟合上包络和下包络计算包络均值用原信号减去包络均值得到筛分后的候选分量。把这个候选分量作为新信号重复上述过程直到满足IMF的两个条件就取走这一层IMF。用公式表达单轮操作就是 h(t) x(t) - m(t)其中m(t)是上下包络的均值h(t)是本次筛分结果。实际分解时第一次筛分得到的h(t)通常已经很像一个IMF但往往还残留一些局部不对称需要反复剥离包络均值。每筛一轮极值点和过零点的数量越来越接近包络均值越来越小波形越来越对称。整个过程可以用几句伪代码概括def sift_once(x, spline_kindakima): # 找到局部极值点需要保证首尾各保留一个端点 maxima_pos, maxima_val find_local_maxima(x) minima_pos, minima_val find_local_minima(x) # 用样条分别拟合上下包络 upper interp1d(maxima_pos, maxima_val, kindspline_kind) lower interp1d(minima_pos, minima_val, kindspline_kind) # 包络均值 mean_env (upper lower) / 2 # 剥离均值返回候选IMF和包络均值供停止准则判据用 return x - mean_env, mean_env这段逻辑对应的是每一次「剥离」的动作。实际操作中要注意极值点首尾两端没有包络值必须做端点延拓处理否则包络会在端点处迅速发散污染整层IMF。PyEMD默认用镜像延拓nbsym2也就是把端点处的极值按对称方式补齐一组虚拟极值这比直接截断要稳得多。三次样条插值的阶数也不是越高越好高阶样条会引入额外的振荡常见做法是三次或Akima样条——Akima样条对过冲更敏感的场景表现更稳。2.2.2 残余分量与分解终止条件提取出一层IMF之后把原始信号减去IMF得到残余信号 r(t) x(t) - IMF1(t)然后对 r(t) 重复整个筛分过程提取IMF2、IMF3……每提走一层残余信号里的极值点数量就减少一批频率成分也越来越低。这和FFT的「一次把所有频率都算出来」的风格完全不同EMD更像是逐层剥离每次只拿走最不平滑的那个成分。分解进行到什么时候停两个条件满足其一即可残余信号变成单调函数再也找不出新的极值对残余信号足够小小到低于预先设定的幅值阈值。第一个条件更常用因为单调信号已经没有振荡可拆再继续筛分只能得到无意义的数值噪声。最终剩下的残余信号表示信号的趋势项它通常不是零均值的振荡而是一个缓慢变化的基线。要注意EMD分解不保证正交性各层IMF之间可能存在能量泄漏这也是为什么后面要引入EEMD、CEEMDAN这类改进算法来对冲模态混叠和端点效应的影响。2.3 停止准则与典型参数SD阈值、最大筛分次数与残差终止筛分迭代不是无限进行的。Huang在提出EMD时给出的停止准则是连续两次筛分结果的差值标准差Sd按公式算的话就是相邻两次候选IMF的差的平方和除以本次候选IMF的平方和再开方。经典建议值是 Sd 在 0.2 到 0.3 之间。但这里有个工程中常踩的坑如果Sd设得偏小也就是要求包络均值归零归得很苛刻筛分次数会急剧增加每一轮都在磨平幅值包络最终得到的IMF会退化成一个频率调制分量瞬时频率倒是光滑了但包络幅值被削成常数幅值调制信息彻底丢失。所以很多实现里除了Sd还会同时限制最大筛分次数比如每层IMF最多迭代50次谁先触发谁停。再看几个关键参数的典型范围参数含义常见取值范围设置不当的表现Sd相邻筛分结果的归一化标准差收敛判据0.1 ~ 0.3偏小则过筛IMF幅值包络被削平偏大则IMF不满足窄带条件max_iter单层IMF的最大筛分次数50 ~ 100限制过小则IMF对称性差瞬时频率出现负值nbsym端点延拓方式0/1/2对应不延拓/镜像/周期复制2不延拓时端点效应污染严重且误差逐层扩散spline_kind包络插值方式linear/cubic/akimacubic与akimalinear包络粗糙高阶样条容易过冲残差终止的判断也要关注当残余信号的极值点数量小于2时直接结束分解。不要等到残余信号的幅值趋近于零才停因为噪声和数值误差会在最后几层制造虚假极值拖慢计算速度且不增加任何信息量。3. Hilbert变换与瞬时频率从IMF到可解释的时频结构3.1 Hilbert变换的90度相移与解析信号Hilbert变换从数学上看是一个线性算子它对信号做的是±90度的相移滤波对正频率成分滞后90度对负频率成分超前90度幅值保持。因为它的频率响应是恒定的相移所以不会改变信号的幅值谱。对一个实信号x(t)做完Hilbert变换后得到x̂(t)可以构造解析信号 z(t) x(t) j·x̂(t)这个复数信号把实信号和它的相移版本组合在一起相当于把原来的单边频谱搬到了正频率轴上。解析信号的物理意义在于它可以唯一地定义每个时刻的瞬时幅值和瞬时相位。瞬时幅值就是解析信号的模瞬时相位就是解析信号的辐角对瞬时相位求时间导数再除以2π就得到瞬时频率。这一步是HHT真正区别于FFT的地方FFT给出的是「整个时间段里存在哪些频率」而Hilbert变换给出的是「每个时刻的那个频率值」。对一段振动信号做故障诊断时瞬时频率曲线能直接反映转频随时间的爬升过程而FFT频谱只能告诉你这段时间内转频大致在哪个范围波动。在Python里用现成库就能完成这一步。scipy.signal.hilbert返回的正是解析信号之后对相位做unwrap防止相位折叠再求梯度就能得到瞬时频率序列。代码如下from scipy.signal import hilbert import numpy as np # imf是EMD分解出来的某个本征模态函数dt为采样间隔 analytic hilbert(imf) inst_amp np.abs(analytic) # 瞬时幅值包络 inst_phase np.unwrap(np.angle(analytic)) # 展开后的瞬时相位 inst_freq np.diff(inst_phase) / (2.0 * np.pi * dt) # 瞬时频率单位Hz瞬时频率序列里偶尔会出现负值这是判断分解质量最直接的信号。负频率没有任何物理意义出现它的常见原因就是输入信号不是窄带的IMF或者是端点效应把相位扰动了。做故障诊断时我会先用这段代码把每层IMF的瞬时频率画出来扫一眼比直接看谱图更快地暴露分解失败的问题。3.2 直接对原始信号求瞬时频率为什么不行理论上Hilbert变换对任何信号都能算出一个瞬时频率但算出来的结果不一定有物理意义。问题出在相位求导这个过程上如果信号不是窄带的相位轨迹会出现局部滞留甚至回折瞬时频率就会大幅震荡甚至出现负数。举个例子一段同时包含工频50Hz和2kHz振动冲击的信号直接做Hilbert变换瞬时频率会在两个频段之间来回跳频率曲线看起来像毛刺序列毫无分析价值。更根本的约束来自Bedrosian定理和Nuttall定理。Bedrosian定理说只有当幅值包络的频谱和相位项的频谱不重叠时解析信号的幅值才能还原真实包络Nuttall定理进一步给出了瞬时频率无偏的条件。实际信号很难天然满足这些约束所以必须先用EMD把信号拆成窄带分量让每个IMF满足过零点数与极值点数接近这一近似窄带条件再做Hilbert变换才可靠。这也是「Hilbert-Huang变换」这个名字的由来Huang的贡献不是Hilbert变换本身而是用EMD给Hilbert变换铺好了可行的输入条件。3.3 Hilbert谱与边际谱的表征逻辑把每一层IMF的瞬时幅值按时间-频率平面铺开就得到Hilbert谱它是二维分布横轴是时间纵轴是瞬时频率颜色深浅表示瞬时幅值能量。这个谱图能清楚看到频率成分随时间的变化轨迹比如变频电机的启动过程中频率曲线从低到高连续爬升谱图上就是一条上升的亮线。作为对比短时傅里叶变换的时频谱是块状的受窗函数限制在低频段频率分辨率差小波变换的时频谱虽然多尺度但母小波的选择会影响结果。Hilbert谱没有固定窗宽频率值直接由相位求导得到因此在频率变化平滑时具有极高的精度。对Hilbert谱沿时间方向积分得到边际谱# 假设inst_freq和inst_amp均为数组freq_bins为预先指定的频率分箱 marginal_spectrum, _ np.histogram(inst_freq, binsfreq_bins, weightsinst_amp)边际谱的物理含义是在整个观测时间范围内每个频率点上累积的能量密度。它和FFT幅值谱最大的区别在于FFT谱的纵轴代表该频率分量在时域上持续振荡的幅值边际谱的纵轴代表该瞬时频率的出现总时长乘以相应幅值所以边际谱能分辨出FFT谱上被平均掉的瞬态事件。对冲击类故障来说边际谱上会在特征频率处出现明显峰值这是FFT频谱很难做到的。4. 用PyEMD在本地跑通EMD-HHT的最小流程4.1 运行环境与测试信号手边有Python环境就能直接上手。需要的核心库是PyEMD、numpy、scipy和matplotlib安装命令一条带过pip install PyEMD scipy numpy matplotlib即可。PyEMD提供的是纯Python的EMD实现没有重型运行时依赖在Windows、macOS、Linux下行为一致这也是我推荐拿它做方案验证的原因。测试信号不适合用纯正弦叠加那体现不出EMD的优势。我会构造一个频率随时间变化的调频分量叠上一个固定频率的低幅分量import numpy as np from PyEMD import EMD from scipy.signal import hilbert import matplotlib.pyplot as plt # 采样参数 fs 500 t np.linspace(0, 2, fs * 2, endpointFalse) # 分量1频率从20Hz扫到60Hz的调频信号 f_inst 20 20 * t x1 1.0 * np.cos(2 * np.pi * (20 * t 10 * t**2)) # 分量2固定80Hz的低幅信号 x2 0.3 * np.sin(2 * np.pi * 80 * t) x x1 x2这段信号的时频结构中有一个明显的扫频轨迹还有一个恒定高频弱分量做FFT只能看到20到60Hz之间的能量涂抹成一片而EMD-HHT应能清楚分出两条频率轨迹。4.2 完整代码EMD分解、瞬时频率计算与谱图绘制EMD分解只用一行但参数需要按需设置。Sd、nbsym、spline_kind三个参数是控制分解质量的关键。我这里用Akima样条配镜像延拓能够明显减少包络过冲# 初始化EMDSd0.2镜像延拓Akima样条 emd EMD(Sd0.2, nbsym2, spline_kindakima) imfs emd(x) # 每一行是一个IMF最后一行为残余 # 计算各IMF的Hilbert谱相关量 dt 1 / fs t_grid t[: -1] # 瞬时频率比原始序列少一个点 for i, imf in enumerate(imfs): analytic hilbert(imf) inst_amp np.abs(analytic) inst_phase np.unwrap(np.angle(analytic)) inst_freq np.diff(inst_phase) / (2 * np.pi * dt) # 这里可以保存每一层的瞬时频率和幅值用于后续绘图说明一下参数行为Sd0.2是Huang论文里的经典收敛阈值工程上很少需要调到0.1以下低于0.1会显著增加筛分次数并削平幅值包络nbsym2指定镜像延拓把端点极值对称地复制到序列外部防止包络在端点发散这是抑制端点效应最基础的手段spline_kindakima相比cubic样条对局部过冲更克制尤其适合数据含噪或幅值变化剧烈的场景。如果IMFs数量明显多于预期优先怀疑是端点延拓没有生效或采样率设置不合理。把瞬时频率按时间画成散点图颜色映射瞬时幅值就是一个基础版的Hilbert谱fig, ax plt.subplots(figsize(12, 4)) for i in range(len(imfs) - 1): # 残余项不做Hilbert处理 analytic hilbert(imfs[i]) inst_amp np.abs(analytic) inst_phase np.unwrap(np.angle(analytic)) inst_freq np.diff(inst_phase) / (2 * np.pi * dt) ax.scatter(t_grid, inst_freq, cinst_amp, s1, cmaphot, alpha0.6) ax.set_xlabel(Time (s)) ax.set_ylabel(Instantaneous Frequency (Hz)) ax.set_ylim(0, 150) plt.tight_layout() plt.show()散点图的读法要说明一下横轴时间纵轴频率颜色代表瞬时能量。如果分解正确图上应该看到一条从20Hz平滑升到60Hz的亮色轨迹和一条稳定在80Hz左右、颜色较暗的水平轨迹。如果出现大量散乱点或负频率说明该层IMF不纯需要回到EMD参数上做调整。边际谱则可以直接对每层瞬时幅值按频率分箱累加得到不需要额外库。4.3 参数设置与结果判读跑通一次分解之后重要的不是「看见了曲线」而是能判断曲线是否可信。我的做法是先分解再分层核验重点看三件事。第一件检查IMF层数和各层之间的频率关系。正常分解下IMF1频率最高后续各层逐步变低残余是单调趋势。如果IMF1里混着明显的高低频成分说明Sd偏大或分解提前终止了如果某一层IMF的频率轨迹出现大段抖动多半是端点效应传染此时把nbsym参数从2换成周期延拓做对照实验。第二件看瞬时频率是否越界。前文代码中瞬时频率的下限和上限应当限制在0到采样率的一半之间。一旦出现负频率几乎可以断定这一层没有完全满足窄带条件需要提高筛分收敛精度也就是调小Sd或调大最大迭代次数。第三件把重构信号和原始信号做差检查残差。把全部IMF连同残余相加应当能完全重建原始信号reconstructed np.sum(imfs, axis0) recon_error np.max(np.abs(reconstructed - x)) print(fMaximum reconstruction error: {recon_error:.3e})如果重构误差大于1e-10量级说明PyEMD内部做了降采样或端点截断这时要核对数据长度是否为2的幂、首尾是否有剧变这两点是重构误差最常见的原因。实际的EMD-HHT参数调优经验可以汇总成一张速查表症状优先调整的参数调整方向备选方案IMF层数偏少分解粗糙Sd / max_iter调小Sd到0.15调大max_iter到200改用spline_kindcubic瞬时频率出现负值Sd / nbsym调小Sdnbsym设为2对IMF做归一化Hilbert变换端点处频率大幅摆动nbsym从2改为周期延拓先用信号延拓算法扩展两端模态混叠明显改用EEMD添加白噪声用CEEMDAN替代EMD包络过冲导致幅值失真spline_kind改为akima提高采样率或做低通预滤波5. 端点效应与模态混叠EMD-HHT的扩散抑制与验证5.1 端点误差如何向IMF内部蔓延EMD的包络拟合依赖极值点而信号两端往往不是极值位置三次样条在外部没有约束点包络在端点附近会翘起或下垂这是端点效应的根源。更麻烦的是EMD是逐层剥离的第一层IMF端点处的误差会原样留在残余信号里下一层IMF继续在同样的位置产生新的误差误差一层层向下传导。所以端点效应不是只在序列两头的小段范围出现而是在每一层IMF里都留下痕迹越到后面的IMF污染区域反而越宽因为高频IMF先被剥离剩余信号振荡周期变大端点误差影响的时间尺度也随之加大。端点效应的抑制有几种常见做法。最省事的是镜像延拓也就是前面提到的nbsym2它把端点最近的一个极值镜像到序列外部相当于人为构造了一段虚拟信号供样条拟合约束。另一种是多项式拟合延拓用信号前几拍的走势外推出端点外的极值位置适合趋势明显的信号。还有AR模型预测延拓在端点处建立自回归模型向前预测几十个点这种方法对平稳段有效对突变段容易放大误差。实际项目里极少有人会为了端点效应单独写一套延拓算法PyEMD的镜像延拓配合序列两端各剔除5到10个点再分析是性价比最高的组合。5.2 模态混叠与EEMD/CEEMDAN的降噪经验模态混叠指的是本应落在不同频段的成分在某一层IMF里混在一起典型例子是低频慢变信号上叠加了一串间歇高频冲击。EMD分解高频冲击时极值点分布不均匀样条包络在冲击段和空闲段之间大幅波动结果在这段区间内的IMF时频轨迹来回跳。这就是Huang后来提出EEMD的直接动机给原始信号反复叠加有限幅值的白噪声利用噪声的均匀极值分布把不同尺度的信号强制分配到不同IMF中最后对所有分解结果做集合平均噪声因为是随机的会相互抵消。EEMD有两个参数与降噪效果直接相关叠加白噪声的幅值倍数Nstd和集合平均次数NE。经验规则是Nstd取原始信号标准差的0.2倍左右NE取200次以上。噪声太小起不到均匀化极值的作用太大则会把真实信号成分打散导致IMF失真。CEEMDAN是EEMD的改进版本它不是每次都在完整信号上加独立噪声而是在每一层分解残差上添加自适应噪声集合平均后残余噪声大幅下降IMF的正交性也更好。在Python里用PyEMD切换到EEMD只需修改三行代码from PyEMD import EEMD eemd EEMD(trials200, noise_width0.2) # trials对应NEnoise_width对应Nstd emds eemd(x) # 返回三维数组不同trials的分解结果按需取均值这里noise_width是噪声标准差与信号标准差的比值而不是噪声的绝对值。对同一份数据做EEMD时的经验是把这个值设为0.1到0.3之间如果信号本身信噪比低就取上限0.3。trial次数的影响比较单调越多越稳但计算时间线性增长。一般200次够用追求更高稳定性可以调到500但很少有必要超过1000。5.3 用归一化的Hilbert变换做可量化的自检最后落一个可以立刻动手的验证技巧。即使参数调得再仔细EMD分解出的IMF也很难是理想窄带信号直接对IMF做Hilbert变换求瞬时频率在幅值包络变化剧烈的区段仍会出现负频点。规范的工程做法是做一步归一化处理先从IMF提取幅值包络E(t)然后把IMF除以E(t)得到纯调频信号F(t) IMF / E(t)对F(t)做Hilbert变换求瞬时频率。因为F(t)的幅值恒为1Hilbert变换受幅值调制的影响被彻底消除得到的瞬时频率曲线更加平滑负值频点数量可以降到接近零。def normalized_hilbert_freq(imf, dt): analytic_env np.abs(hilbert(imf)) # 先估计幅值包络 fm_part imf / (analytic_env 1e-12) # 归一化避免除零 analytic_fm hilbert(fm_part) # 对纯调频分量做Hilbert inst_phase np.unwrap(np.angle(analytic_fm)) return np.diff(inst_phase) / (2 * np.pi * dt)判断分解质量时统计一下这段代码算出的负频率点占比如果在1%以下说明该层IMF质量过关如果超过5%先把这层IMF单独画出来和上一节的参数速查表对照多半是出现了模态混叠或过筛。这个归一化Hilbert变换与直接Hilbert变换之间的差异本身就是对IMF纯度的量化估计两者频率曲线差异越小说明EMD的分解越接近理论理想也说明你手里这套参数对这组数据靠谱。对一批传感器数据批量做HHT分析时我会先抽一条数据做这个自检再决定后续是直接用EMD还是升级到EEMD/CEEMDAN而不是一上来就开高参数跑全量数据。本文还有配套的精品资源点击获取
EMD与HHT详解：经验模态分解、IMF提取与瞬时频率分析
简介基于经验模态分解法EMD的Hilbert-Huang变换HHTMATLAB程序包面向信号处理、故障诊断及振动分析领域的研究人员和工程师用于将非平稳信号转化为平稳信号进而完成时频分析。程序通过EMD将信号分解为若干本征模态函数IMF分量再对IMF分量进行Hilbert变换以获取瞬时频率与幅值最终累加重构得平稳化结果核心逻辑清晰层次直观便于理解HHT算法链路并做二次开发。压缩包内共3个文件均为.m源码包含主程序HHT.m以及频谱计算、瞬时频率估计所需的hhspectrum.m与instfreq.m函数整体仅2KB体量极小运行轻便需配合已安装的EMD工具箱中的emd函数使用。目前已有1277人学习/浏览适合正在学习HHT原理、搭建MATLAB分析流程或复现场景的读者。借助这份代码可以快速掌握IMF分量提取、Hilbert谱构建、瞬时频率计算等关键步骤并直接嵌入自己的实验或项目中减少从零编写的时间成本。1. EMD与HHT从非平稳信号里拆出可解释的时频结构轴承振动、脑电波、风速序列、结构应变这些信号有个共同点频率成分随时间漂移。拿一段转速爬升过程中的振动数据做FFT频谱上能量铺成一片峰值对应的频率说不清是哪个时刻的加窗做短时傅里叶变换又要在时间分辨率和频率分辨率之间做取舍窗宽一固定整个时频面只能用同一把尺子量。经验模态分解Empirical Mode Decomposition, EMD搭配Hilbert-Huang变换HHT走的是另一条路先按信号自身的极值分布把数据分解成若干本征模态函数Intrinsic Mode Function, IMF再对每个IMF做Hilbert变换求出有物理意义的瞬时频率从而得到一条沿时间轴变化的频率曲线。整个过程不预设基函数不要求信号平稳也不需要选择窗函数这两点让它和FFT、小波变换在思路上彻底分开。这套方法适合做振动故障诊断、结构健康监测、生物医学信号分析的技术人员也适合所有对「频率随时间怎么变」有量化需求的场景。下面从EMD的筛分逻辑讲起一路落到能跑的代码和参数调优。2. EMD核心流程从原始信号到IMF的筛分步骤与停止准则2.1 IMF的两个定义性条件过零点约束与包络对称把一个信号拆成多个IMF前提是得先定义清楚「拆出来的每一层到底长什么样」。IMF必须同时满足两个条件一是整个数据段内过零点的数量和极值点的数量相等或至多相差一个二是在任意时间点上由局部极大值拟合出的上包络和由局部极小值拟合出的下包络的均值必须为零也就是上下包络关于时间轴对称。第一条保证了IMF是窄带振荡第二条保证瞬时频率不会因为包络不对称而出现无意义的跳动。直观理解就是IMF不能像原始信号那样既包含慢变的大幅趋势又叠着快速的小幅抖动它必须是一段「围绕零点、上下包络均衡」的振荡波形。这里可以对比一下傅里叶分解FFT把信号投影到无限长的正弦波基上任何一个局部扰动都会摊到整个频谱上而EMD的IMF是从数据本身长出来的极值点密集的地方就分解出高频IMF极值稀疏的地方就落到低频IMF基函数跟着数据走这是「经验」二字的由来。2.2 筛分迭代极值定位、三次样条包络与均值剥离2.2.1 单轮筛分的三个动作EMD对信号的处理不是一次性分解而是反复迭代。每一轮筛分做三件事找出当前信号的所有局部极大值和局部极小值点用三次样条曲线分别拟合上包络和下包络计算包络均值用原信号减去包络均值得到筛分后的候选分量。把这个候选分量作为新信号重复上述过程直到满足IMF的两个条件就取走这一层IMF。用公式表达单轮操作就是 h(t) x(t) - m(t)其中m(t)是上下包络的均值h(t)是本次筛分结果。实际分解时第一次筛分得到的h(t)通常已经很像一个IMF但往往还残留一些局部不对称需要反复剥离包络均值。每筛一轮极值点和过零点的数量越来越接近包络均值越来越小波形越来越对称。整个过程可以用几句伪代码概括def sift_once(x, spline_kindakima): # 找到局部极值点需要保证首尾各保留一个端点 maxima_pos, maxima_val find_local_maxima(x) minima_pos, minima_val find_local_minima(x) # 用样条分别拟合上下包络 upper interp1d(maxima_pos, maxima_val, kindspline_kind) lower interp1d(minima_pos, minima_val, kindspline_kind) # 包络均值 mean_env (upper lower) / 2 # 剥离均值返回候选IMF和包络均值供停止准则判据用 return x - mean_env, mean_env这段逻辑对应的是每一次「剥离」的动作。实际操作中要注意极值点首尾两端没有包络值必须做端点延拓处理否则包络会在端点处迅速发散污染整层IMF。PyEMD默认用镜像延拓nbsym2也就是把端点处的极值按对称方式补齐一组虚拟极值这比直接截断要稳得多。三次样条插值的阶数也不是越高越好高阶样条会引入额外的振荡常见做法是三次或Akima样条——Akima样条对过冲更敏感的场景表现更稳。2.2.2 残余分量与分解终止条件提取出一层IMF之后把原始信号减去IMF得到残余信号 r(t) x(t) - IMF1(t)然后对 r(t) 重复整个筛分过程提取IMF2、IMF3……每提走一层残余信号里的极值点数量就减少一批频率成分也越来越低。这和FFT的「一次把所有频率都算出来」的风格完全不同EMD更像是逐层剥离每次只拿走最不平滑的那个成分。分解进行到什么时候停两个条件满足其一即可残余信号变成单调函数再也找不出新的极值对残余信号足够小小到低于预先设定的幅值阈值。第一个条件更常用因为单调信号已经没有振荡可拆再继续筛分只能得到无意义的数值噪声。最终剩下的残余信号表示信号的趋势项它通常不是零均值的振荡而是一个缓慢变化的基线。要注意EMD分解不保证正交性各层IMF之间可能存在能量泄漏这也是为什么后面要引入EEMD、CEEMDAN这类改进算法来对冲模态混叠和端点效应的影响。2.3 停止准则与典型参数SD阈值、最大筛分次数与残差终止筛分迭代不是无限进行的。Huang在提出EMD时给出的停止准则是连续两次筛分结果的差值标准差Sd按公式算的话就是相邻两次候选IMF的差的平方和除以本次候选IMF的平方和再开方。经典建议值是 Sd 在 0.2 到 0.3 之间。但这里有个工程中常踩的坑如果Sd设得偏小也就是要求包络均值归零归得很苛刻筛分次数会急剧增加每一轮都在磨平幅值包络最终得到的IMF会退化成一个频率调制分量瞬时频率倒是光滑了但包络幅值被削成常数幅值调制信息彻底丢失。所以很多实现里除了Sd还会同时限制最大筛分次数比如每层IMF最多迭代50次谁先触发谁停。再看几个关键参数的典型范围参数含义常见取值范围设置不当的表现Sd相邻筛分结果的归一化标准差收敛判据0.1 ~ 0.3偏小则过筛IMF幅值包络被削平偏大则IMF不满足窄带条件max_iter单层IMF的最大筛分次数50 ~ 100限制过小则IMF对称性差瞬时频率出现负值nbsym端点延拓方式0/1/2对应不延拓/镜像/周期复制2不延拓时端点效应污染严重且误差逐层扩散spline_kind包络插值方式linear/cubic/akimacubic与akimalinear包络粗糙高阶样条容易过冲残差终止的判断也要关注当残余信号的极值点数量小于2时直接结束分解。不要等到残余信号的幅值趋近于零才停因为噪声和数值误差会在最后几层制造虚假极值拖慢计算速度且不增加任何信息量。3. Hilbert变换与瞬时频率从IMF到可解释的时频结构3.1 Hilbert变换的90度相移与解析信号Hilbert变换从数学上看是一个线性算子它对信号做的是±90度的相移滤波对正频率成分滞后90度对负频率成分超前90度幅值保持。因为它的频率响应是恒定的相移所以不会改变信号的幅值谱。对一个实信号x(t)做完Hilbert变换后得到x̂(t)可以构造解析信号 z(t) x(t) j·x̂(t)这个复数信号把实信号和它的相移版本组合在一起相当于把原来的单边频谱搬到了正频率轴上。解析信号的物理意义在于它可以唯一地定义每个时刻的瞬时幅值和瞬时相位。瞬时幅值就是解析信号的模瞬时相位就是解析信号的辐角对瞬时相位求时间导数再除以2π就得到瞬时频率。这一步是HHT真正区别于FFT的地方FFT给出的是「整个时间段里存在哪些频率」而Hilbert变换给出的是「每个时刻的那个频率值」。对一段振动信号做故障诊断时瞬时频率曲线能直接反映转频随时间的爬升过程而FFT频谱只能告诉你这段时间内转频大致在哪个范围波动。在Python里用现成库就能完成这一步。scipy.signal.hilbert返回的正是解析信号之后对相位做unwrap防止相位折叠再求梯度就能得到瞬时频率序列。代码如下from scipy.signal import hilbert import numpy as np # imf是EMD分解出来的某个本征模态函数dt为采样间隔 analytic hilbert(imf) inst_amp np.abs(analytic) # 瞬时幅值包络 inst_phase np.unwrap(np.angle(analytic)) # 展开后的瞬时相位 inst_freq np.diff(inst_phase) / (2.0 * np.pi * dt) # 瞬时频率单位Hz瞬时频率序列里偶尔会出现负值这是判断分解质量最直接的信号。负频率没有任何物理意义出现它的常见原因就是输入信号不是窄带的IMF或者是端点效应把相位扰动了。做故障诊断时我会先用这段代码把每层IMF的瞬时频率画出来扫一眼比直接看谱图更快地暴露分解失败的问题。3.2 直接对原始信号求瞬时频率为什么不行理论上Hilbert变换对任何信号都能算出一个瞬时频率但算出来的结果不一定有物理意义。问题出在相位求导这个过程上如果信号不是窄带的相位轨迹会出现局部滞留甚至回折瞬时频率就会大幅震荡甚至出现负数。举个例子一段同时包含工频50Hz和2kHz振动冲击的信号直接做Hilbert变换瞬时频率会在两个频段之间来回跳频率曲线看起来像毛刺序列毫无分析价值。更根本的约束来自Bedrosian定理和Nuttall定理。Bedrosian定理说只有当幅值包络的频谱和相位项的频谱不重叠时解析信号的幅值才能还原真实包络Nuttall定理进一步给出了瞬时频率无偏的条件。实际信号很难天然满足这些约束所以必须先用EMD把信号拆成窄带分量让每个IMF满足过零点数与极值点数接近这一近似窄带条件再做Hilbert变换才可靠。这也是「Hilbert-Huang变换」这个名字的由来Huang的贡献不是Hilbert变换本身而是用EMD给Hilbert变换铺好了可行的输入条件。3.3 Hilbert谱与边际谱的表征逻辑把每一层IMF的瞬时幅值按时间-频率平面铺开就得到Hilbert谱它是二维分布横轴是时间纵轴是瞬时频率颜色深浅表示瞬时幅值能量。这个谱图能清楚看到频率成分随时间的变化轨迹比如变频电机的启动过程中频率曲线从低到高连续爬升谱图上就是一条上升的亮线。作为对比短时傅里叶变换的时频谱是块状的受窗函数限制在低频段频率分辨率差小波变换的时频谱虽然多尺度但母小波的选择会影响结果。Hilbert谱没有固定窗宽频率值直接由相位求导得到因此在频率变化平滑时具有极高的精度。对Hilbert谱沿时间方向积分得到边际谱# 假设inst_freq和inst_amp均为数组freq_bins为预先指定的频率分箱 marginal_spectrum, _ np.histogram(inst_freq, binsfreq_bins, weightsinst_amp)边际谱的物理含义是在整个观测时间范围内每个频率点上累积的能量密度。它和FFT幅值谱最大的区别在于FFT谱的纵轴代表该频率分量在时域上持续振荡的幅值边际谱的纵轴代表该瞬时频率的出现总时长乘以相应幅值所以边际谱能分辨出FFT谱上被平均掉的瞬态事件。对冲击类故障来说边际谱上会在特征频率处出现明显峰值这是FFT频谱很难做到的。4. 用PyEMD在本地跑通EMD-HHT的最小流程4.1 运行环境与测试信号手边有Python环境就能直接上手。需要的核心库是PyEMD、numpy、scipy和matplotlib安装命令一条带过pip install PyEMD scipy numpy matplotlib即可。PyEMD提供的是纯Python的EMD实现没有重型运行时依赖在Windows、macOS、Linux下行为一致这也是我推荐拿它做方案验证的原因。测试信号不适合用纯正弦叠加那体现不出EMD的优势。我会构造一个频率随时间变化的调频分量叠上一个固定频率的低幅分量import numpy as np from PyEMD import EMD from scipy.signal import hilbert import matplotlib.pyplot as plt # 采样参数 fs 500 t np.linspace(0, 2, fs * 2, endpointFalse) # 分量1频率从20Hz扫到60Hz的调频信号 f_inst 20 20 * t x1 1.0 * np.cos(2 * np.pi * (20 * t 10 * t**2)) # 分量2固定80Hz的低幅信号 x2 0.3 * np.sin(2 * np.pi * 80 * t) x x1 x2这段信号的时频结构中有一个明显的扫频轨迹还有一个恒定高频弱分量做FFT只能看到20到60Hz之间的能量涂抹成一片而EMD-HHT应能清楚分出两条频率轨迹。4.2 完整代码EMD分解、瞬时频率计算与谱图绘制EMD分解只用一行但参数需要按需设置。Sd、nbsym、spline_kind三个参数是控制分解质量的关键。我这里用Akima样条配镜像延拓能够明显减少包络过冲# 初始化EMDSd0.2镜像延拓Akima样条 emd EMD(Sd0.2, nbsym2, spline_kindakima) imfs emd(x) # 每一行是一个IMF最后一行为残余 # 计算各IMF的Hilbert谱相关量 dt 1 / fs t_grid t[: -1] # 瞬时频率比原始序列少一个点 for i, imf in enumerate(imfs): analytic hilbert(imf) inst_amp np.abs(analytic) inst_phase np.unwrap(np.angle(analytic)) inst_freq np.diff(inst_phase) / (2 * np.pi * dt) # 这里可以保存每一层的瞬时频率和幅值用于后续绘图说明一下参数行为Sd0.2是Huang论文里的经典收敛阈值工程上很少需要调到0.1以下低于0.1会显著增加筛分次数并削平幅值包络nbsym2指定镜像延拓把端点极值对称地复制到序列外部防止包络在端点发散这是抑制端点效应最基础的手段spline_kindakima相比cubic样条对局部过冲更克制尤其适合数据含噪或幅值变化剧烈的场景。如果IMFs数量明显多于预期优先怀疑是端点延拓没有生效或采样率设置不合理。把瞬时频率按时间画成散点图颜色映射瞬时幅值就是一个基础版的Hilbert谱fig, ax plt.subplots(figsize(12, 4)) for i in range(len(imfs) - 1): # 残余项不做Hilbert处理 analytic hilbert(imfs[i]) inst_amp np.abs(analytic) inst_phase np.unwrap(np.angle(analytic)) inst_freq np.diff(inst_phase) / (2 * np.pi * dt) ax.scatter(t_grid, inst_freq, cinst_amp, s1, cmaphot, alpha0.6) ax.set_xlabel(Time (s)) ax.set_ylabel(Instantaneous Frequency (Hz)) ax.set_ylim(0, 150) plt.tight_layout() plt.show()散点图的读法要说明一下横轴时间纵轴频率颜色代表瞬时能量。如果分解正确图上应该看到一条从20Hz平滑升到60Hz的亮色轨迹和一条稳定在80Hz左右、颜色较暗的水平轨迹。如果出现大量散乱点或负频率说明该层IMF不纯需要回到EMD参数上做调整。边际谱则可以直接对每层瞬时幅值按频率分箱累加得到不需要额外库。4.3 参数设置与结果判读跑通一次分解之后重要的不是「看见了曲线」而是能判断曲线是否可信。我的做法是先分解再分层核验重点看三件事。第一件检查IMF层数和各层之间的频率关系。正常分解下IMF1频率最高后续各层逐步变低残余是单调趋势。如果IMF1里混着明显的高低频成分说明Sd偏大或分解提前终止了如果某一层IMF的频率轨迹出现大段抖动多半是端点效应传染此时把nbsym参数从2换成周期延拓做对照实验。第二件看瞬时频率是否越界。前文代码中瞬时频率的下限和上限应当限制在0到采样率的一半之间。一旦出现负频率几乎可以断定这一层没有完全满足窄带条件需要提高筛分收敛精度也就是调小Sd或调大最大迭代次数。第三件把重构信号和原始信号做差检查残差。把全部IMF连同残余相加应当能完全重建原始信号reconstructed np.sum(imfs, axis0) recon_error np.max(np.abs(reconstructed - x)) print(fMaximum reconstruction error: {recon_error:.3e})如果重构误差大于1e-10量级说明PyEMD内部做了降采样或端点截断这时要核对数据长度是否为2的幂、首尾是否有剧变这两点是重构误差最常见的原因。实际的EMD-HHT参数调优经验可以汇总成一张速查表症状优先调整的参数调整方向备选方案IMF层数偏少分解粗糙Sd / max_iter调小Sd到0.15调大max_iter到200改用spline_kindcubic瞬时频率出现负值Sd / nbsym调小Sdnbsym设为2对IMF做归一化Hilbert变换端点处频率大幅摆动nbsym从2改为周期延拓先用信号延拓算法扩展两端模态混叠明显改用EEMD添加白噪声用CEEMDAN替代EMD包络过冲导致幅值失真spline_kind改为akima提高采样率或做低通预滤波5. 端点效应与模态混叠EMD-HHT的扩散抑制与验证5.1 端点误差如何向IMF内部蔓延EMD的包络拟合依赖极值点而信号两端往往不是极值位置三次样条在外部没有约束点包络在端点附近会翘起或下垂这是端点效应的根源。更麻烦的是EMD是逐层剥离的第一层IMF端点处的误差会原样留在残余信号里下一层IMF继续在同样的位置产生新的误差误差一层层向下传导。所以端点效应不是只在序列两头的小段范围出现而是在每一层IMF里都留下痕迹越到后面的IMF污染区域反而越宽因为高频IMF先被剥离剩余信号振荡周期变大端点误差影响的时间尺度也随之加大。端点效应的抑制有几种常见做法。最省事的是镜像延拓也就是前面提到的nbsym2它把端点最近的一个极值镜像到序列外部相当于人为构造了一段虚拟信号供样条拟合约束。另一种是多项式拟合延拓用信号前几拍的走势外推出端点外的极值位置适合趋势明显的信号。还有AR模型预测延拓在端点处建立自回归模型向前预测几十个点这种方法对平稳段有效对突变段容易放大误差。实际项目里极少有人会为了端点效应单独写一套延拓算法PyEMD的镜像延拓配合序列两端各剔除5到10个点再分析是性价比最高的组合。5.2 模态混叠与EEMD/CEEMDAN的降噪经验模态混叠指的是本应落在不同频段的成分在某一层IMF里混在一起典型例子是低频慢变信号上叠加了一串间歇高频冲击。EMD分解高频冲击时极值点分布不均匀样条包络在冲击段和空闲段之间大幅波动结果在这段区间内的IMF时频轨迹来回跳。这就是Huang后来提出EEMD的直接动机给原始信号反复叠加有限幅值的白噪声利用噪声的均匀极值分布把不同尺度的信号强制分配到不同IMF中最后对所有分解结果做集合平均噪声因为是随机的会相互抵消。EEMD有两个参数与降噪效果直接相关叠加白噪声的幅值倍数Nstd和集合平均次数NE。经验规则是Nstd取原始信号标准差的0.2倍左右NE取200次以上。噪声太小起不到均匀化极值的作用太大则会把真实信号成分打散导致IMF失真。CEEMDAN是EEMD的改进版本它不是每次都在完整信号上加独立噪声而是在每一层分解残差上添加自适应噪声集合平均后残余噪声大幅下降IMF的正交性也更好。在Python里用PyEMD切换到EEMD只需修改三行代码from PyEMD import EEMD eemd EEMD(trials200, noise_width0.2) # trials对应NEnoise_width对应Nstd emds eemd(x) # 返回三维数组不同trials的分解结果按需取均值这里noise_width是噪声标准差与信号标准差的比值而不是噪声的绝对值。对同一份数据做EEMD时的经验是把这个值设为0.1到0.3之间如果信号本身信噪比低就取上限0.3。trial次数的影响比较单调越多越稳但计算时间线性增长。一般200次够用追求更高稳定性可以调到500但很少有必要超过1000。5.3 用归一化的Hilbert变换做可量化的自检最后落一个可以立刻动手的验证技巧。即使参数调得再仔细EMD分解出的IMF也很难是理想窄带信号直接对IMF做Hilbert变换求瞬时频率在幅值包络变化剧烈的区段仍会出现负频点。规范的工程做法是做一步归一化处理先从IMF提取幅值包络E(t)然后把IMF除以E(t)得到纯调频信号F(t) IMF / E(t)对F(t)做Hilbert变换求瞬时频率。因为F(t)的幅值恒为1Hilbert变换受幅值调制的影响被彻底消除得到的瞬时频率曲线更加平滑负值频点数量可以降到接近零。def normalized_hilbert_freq(imf, dt): analytic_env np.abs(hilbert(imf)) # 先估计幅值包络 fm_part imf / (analytic_env 1e-12) # 归一化避免除零 analytic_fm hilbert(fm_part) # 对纯调频分量做Hilbert inst_phase np.unwrap(np.angle(analytic_fm)) return np.diff(inst_phase) / (2 * np.pi * dt)判断分解质量时统计一下这段代码算出的负频率点占比如果在1%以下说明该层IMF质量过关如果超过5%先把这层IMF单独画出来和上一节的参数速查表对照多半是出现了模态混叠或过筛。这个归一化Hilbert变换与直接Hilbert变换之间的差异本身就是对IMF纯度的量化估计两者频率曲线差异越小说明EMD的分解越接近理论理想也说明你手里这套参数对这组数据靠谱。对一批传感器数据批量做HHT分析时我会先抽一条数据做这个自检再决定后续是直接用EMD还是升级到EEMD/CEEMDAN而不是一上来就开高参数跑全量数据。本文还有配套的精品资源点击获取