简介面向非线性、非平稳信号分析的学习者与工程技术人员这份希尔伯特黄变换HHT配套资料以经验模态分解与希尔伯特变换为主线覆盖从算法理解到MATLAB实现的完整链路适合生物医学信号、地震数据、机械故障诊断等场景的入门与进阶练习。压缩包共72个文件以40个m函数、17个c源文件、11个头文件为主辅以安装卸载脚本、示例目录与数据文件整体约94KB便于在MATLAB环境中直接调用和二次修改。已有5077人学习下载说明其在HHT实践群体中具有一定参考价值。资料围绕EMD迭代筛选、IMF与残差、包络与边界条件、希尔伯特谱及瞬时频率与振幅计算展开并配有可视化与示例模块可帮助读者理解信号多尺度分解过程、复现时频分析结果并积累函数组织与调试思路。1. 从一段升速振动数据说起HHT 能给你 FFT 给不了的那张图手上拿到过一段变转速设备的升速振动数据做 FFT 得到的是糊成一坨的谱线几个阶次成分的能量全挤在相邻频率栅格里你没法判断某个频率到底出现在第几秒。希尔伯特黄变换HHT就是冲这个场景来的先用经验模态分解EMD把一维信号拆成一族本征模态函数IMF再对每个 IMF 做希尔伯特变换拿到随时间变化的瞬时频率和瞬时幅值拼成时频图与边际谱。这份资源围绕的就是 EMD 加希尔伯特谱这条完整链路从分解、筛分停止准则、端点处理一直做到瞬时频率修正、边际谱和时频图渲染适合手上有一维非平稳序列、又不想在基函数和窗口宽度上反复试错的人做旋转机械诊断、地震波形、生理信号分析都能直接改参数上手。2. EMD 分解怎么跑通从原始序列到 IMF 族的可复现脚本2.1 为什么先分解再做希尔伯特而不是直接 FFTFFT 的基是铺满整个时间轴的正弦一段两秒的数据里频率从 20 Hz 扫到 120 Hz能量就会被摊到一大片频点上谱线看起来又宽又虚。短时傅里叶加窗能补一部分但窗口一宽时间分辨率就没了窗口一窄频率分辨率又崩窗口长度本质上是在两种模糊之间做取舍。小波变换往前走了一步可基函数选型和分解层数仍然是人为设定的选错了照样拆不出物理上有意义的分量。EMD 换了个思路不预设基函数直接从信号自身的极大值点和极小值点出发用三次样条拟合出上包络和下包络取均值当成低频趋势剥掉剩下的高频部分继续做同样的操作如此反复筛分一层层剥出 IMF。每个 IMF 需要满足两条约定极值点个数与过零点个数相差不超过 1上下包络的均值为零。这两条约定看着简单实际就是保证分解出来的分量在局部近似单分量后面用希尔伯特变换求瞬时频率才不会出现频率乱跳。代价也很明确。EMD 不是唯一分解同一段数据换一组停止准则、换一个端点处理策略IMF 的个数和形状都会变。所以参数必须写进报告不能只留一张图。这一点和 FFT 那种输入定了输出就定的变换完全不是一回事也是后面所有坑的根源。2.2 用 Python 跑通一次 EMD参数、停止准则和残差下面这段脚本可以直接跑用扫频信号加一个间歇分量做输入目的是让分解结果里同时出现连续分量和间歇分量方便观察模态混叠。import numpy as np from PyEMD import EMD fs 1000.0 # 采样率换成实际数据的值 t np.arange(0, 2, 1/fs) # 构造一段 20 Hz 扫到 120 Hz 的调频信号 f0, f1 20.0, 120.0 phase 2*np.pi*(f0*t 0.5*(f1-f0)/t[-1]*t**2) s1 np.sin(phase) # 间歇出现的 180 Hz 分量只在包络为正的时间段存在 s2 0.6*np.sin(2*np.pi*180*t) * (np.sin(2*np.pi*1.5*t) 0) sig s1 s2 0.05*np.random.randn(t.size) emd EMD() emd.FIXE_H 5 # 每个 IMF 的筛分次数上限控制分解深度与耗时 emd.nbsym 4 # 端点镜像延拓的极值点个数越大端点越稳但越容易吃掉真实极值 imfs emd.emd(sig, max_imf6) res sig - imfs.sum(axis0) print(IMF 矩阵形状:, imfs.shape) for i, imf in enumerate(imfs): # 极值点近似计数用来判断该分量是否还满足 IMF 条件 n_ext np.sum(np.diff(np.sign(np.diff(imf))) ! 0) print(fIMF{i1} 标准差{imf.std():.4f} 极值点数{n_ext}) print(残差能量占比: %.4f % (np.sum(res**2)/np.sum(sig**2)))这段代码里有三个参数值得单独说。FIXE_H决定筛分循环最多执行几次设得太小包络均值没消干净一个 IMF 里会混进相邻尺度设得太大不但耗时翻倍还会把同一物理分量越筛越碎最后拆出一堆没有物理意义的碎片。工程上一般取 3 到 10 之间先用 5 试看各 IMF 的极值点分布是否合理再微调。nbsym控制端点镜像延拓的范围单位是极值点个数。三次样条在两端外推本来就是发散的重灾区镜像点数少端点包络会剧烈上翘分解出的 IMF 首尾出现明显飞边镜像点数多了真实极值被镜像覆盖低频分量会被削掉一截。默认值 2 在平稳段足够用数据两端本身就有大跳变时再往上加到 4 或 6。max_imf是硬上限防止某些数据在噪声上无限剥下去。设成 6 意味着最多返回 6 个 IMF剩下的全部推给残差。调试阶段可以先设大一点看自然分解出多少层正式跑的时候再定死保证同一批数据的结果可比。残差能量占比是一个很直接的验收指标。理想的分解里残差应该只剩单调趋势能量占比通常在百分之几量级如果残差能量占到两成以上说明低频分量被提前截断了要么max_imf太小要么筛分次数不够。2.3 IMF 判读正交性、能量占比与伪分量剔除分解完拿到一堆 IMF第一件事不是画图是判断哪个该留、哪个该扔。常用的判据有三个。正交性指标衡量各 IMF 之间是否真的互不重叠按IO Σ_i Σ_{j≠i} Σ_n imf_i(n)·imf_j(n) / Σ_n sig²(n)计算实际数据里这个值落在 0.05 以内算正常超过 0.2 往往意味着两个 IMF 在时间上有大量重叠也就是发生了模态混叠。能量占比看每个 IMF 的方差占原始信号方差的比例排名前两三个分量通常承载主要物理信息排在末尾、能量不到百分之一的 IMF 基本都是噪声筛出来的伪分量。最实用的判据其实是最朴素的一个 IMF 如果极值点少于 3 个它已经不满足 IMF 的定义只是被算法硬凑出来的一段趋势一个 IMF 如果和原始信号的相关系数低于 0.05基本可以判定为伪分量。这两条组合起来用比任何复杂指标都快。def imf_report(sig, imfs): var_sig np.var(sig) rows [] for i, imf in enumerate(imfs): corr np.corrcoef(sig, imf)[0, 1] n_ext np.sum(np.diff(np.sign(np.diff(imf))) ! 0) rows.append((i1, imf.var()/var_sig, corr, n_ext)) return rows for r in imf_report(sig, imfs): print(IMF%d 能量占比%.4f 相关系数%.4f 极值点数%d % r)输出里如果某个 IMF 能量占比低于 0.005 且相关系数接近零它就是噪声。到这里可以做一个裁剪把伪分量剔掉再进入希尔伯特变换能省掉后面一大半的异常排查。注意裁剪是单向操作。剔掉的成分不影响已保留 IMF 的形状但会影响残差和边际谱的幅值归一化裁剪记录要跟着结果一起存档。3. 希尔伯特变换与瞬时频率把 IMF 拼成时频谱和边际谱3.1 解析信号的构造与瞬时频率的合法性边界对单个 IMF 做希尔伯特变换构造解析信号z(t) imf(t) j·H[imf(t)]包络是模相位是辐角相位对时间求导再除以 2π 就是瞬时频率。这一步数学上没难度难的是它成立的前提。经典的 Bedrosian 定理要求信号是窄带单分量否则解析信号的实部和虚部会出现频谱重叠解出来的瞬时频率可能出现负值或者在两个真实分量之间来回跳。EMD 拆出来的 IMF 只是近似满足单分量条件所以负频率在中后段 IMF 上出现得很常见尤其是当分量幅值很小、信噪比很低的时候。工程上的处理策略有两条。一是直接剔除负频样本把瞬时频率低于零的点丢弃代价是时频图上会留空但不会污染分布二是做归一化希尔伯特变换先把 IMF 除以自身的包络让每个点变成单位幅值的调频信号再做微分求瞬时频率负频率会大幅减少。两条都行做故障诊断类分析时我一般走第二条因为后段 IMF 的负频比例经常能到 5% 以上直接剔除会在关键时刻丢数据。from scipy.signal import hilbert def instantaneous(imf, fs, normalizeFalse): if normalize: env np.abs(hilbert(imf)) env[env 1e-10] 1e-10 # 防止除零包络过零点很常见 imf imf / env z hilbert(imf) amp np.abs(z) ph np.unwrap(np.angle(z)) freq np.diff(ph, prependph[0]) * fs / (2*np.pi) return amp, freq for i, imf in enumerate(imfs[:3]): amp, f instantaneous(imf, fs, normalizeTrue) print(IMF%d 频率中位数%.2f Hz 负频比例%.4f % (i1, np.median(f), np.mean(f 0)))np.unwrap不能省。相位在 ±π 处会发生跳变不做展开差分出来的瞬时频率会出现一两个采样点的尖峰幅度大到能把整个图的色标撑爆。prependph[0]是为了让差分后的长度和原序列对齐方便后面按时间索引累加。归一化分支里的1e-10是保护值包络过零时直接除会得到 inf污染整条频率曲线。负频比例可以作为分解质量的旁证指标。前两个 IMF 的负频比例通常在 0.01 以下到了第四第五个往上跑是正常的说明这些分量已经接近噪声基底做谱分析时权重该往下压。3.2 边际谱与瞬时能量谱的计算脚本希尔伯特谱是一个二维分布H(f, t)把每个时刻每个 IMF 的瞬时幅值按瞬时频率投到频率轴上再沿时间轴累加得到的就是边际谱。它和 FFT 幅值谱形式相似但含义不同边际谱反映的是每个频率上的累计幅值贡献没有傅里叶变换里能量乘以 2 再开根号的那套约定所以两者对比前必须做归一化处理。def hilbert_spectrum(imfs, fs, n_bins512): f_max fs / 2 edges np.linspace(0, f_max, n_bins 1) centers 0.5 * (edges[:-1] edges[1:]) H np.zeros((n_bins, imfs.shape[1])) for imf in imfs: amp, f instantaneous(imf, fs, normalizeTrue) idx np.digitize(f, edges) - 1 ok (idx 0) (idx n_bins) np.add.at(H, (idx[ok], np.arange(f.size)[ok]), amp[ok]) marginal H.sum(axis1) / imfs.shape[1] # 按样本点归一化便于和幅值谱对照 return centers, H, marginalnp.digitize把连续的瞬时频率落进离散栅格np.add.at是因为同一频率格在相邻时刻会被反复命中普通索引赋值会把前面累加的值覆盖掉这是最容易被忽略的一行。如果图上的谱线出现莫名其妙的空洞八成就是这里用了H[idx, t] amp的写法。频率栅格数n_bins和 FFT 的点数没有关系它只是一个显示和统计粒度。想要和 FFT 幅值谱并排比较就把n_bins设为 FFT 点数的一半让两者的频率分辨率对齐否则一边分辨率高一边低形状永远对不上白折腾。瞬时能量谱是另一种看法把每个 IMF 的包络平方沿频率轴累加得到能量随时间的分布。它对冲击类信号特别敏感设备出现周期性冲击时能量谱上会出现等间隔的峰值间隔对应冲击重复频率这个信息在边际谱里是看不到的。3.3 时频图的频率栅格、色标与端点截断时频图最容易出问题的地方不是算法是显示。三个参数决定图能不能看。频率轴上限别偷懒设成 fs/2。高阶 IMF 的瞬时频率经常贴着奈奎斯特频率乱跳直接把主轴撑到 500 Hz真正的低频成分全挤在最下面几个像素里。按信号主频的 3 到 5 倍设上限或者干脆用分位数截断把 99 分位以上的频率值剔掉。色标必须用对数或者分位数截断。幅值动态范围动辄三四十年线性色标下只有最强的几个点可见其余全是背景色。用np.percentile(H[H0], 99)作为上限再np.log10(H eps)低幅值区域的结构才会浮出来。时间轴两端各截掉 5% 到 10%。端点效应是 EMD 改不掉的问题即便做了镜像延拓首尾各几十个采样点的瞬时频率仍然会异常偏高这段数据在图上表现为两条竖直的高频亮带非常干扰判读。截断是最经济的做法比各种端点延拓算法都稳。import matplotlib.pyplot as plt centers, H, marginal hilbert_spectrum(imfs[:4], fs, n_bins256) vmax np.percentile(H[H 0], 99) cut int(0.05 * H.shape[1]) plt.pcolormesh(t[cut:-cut], centers, np.log10(H[:, cut:-cut] 1e-6), vminnp.log10(vmax) - 3, vmaxnp.log10(vmax), shadingauto) plt.ylim(0, 300) # 频率上限按主频的 3~5 倍设置 plt.xlabel(时间 / s); plt.ylabel(频率 / Hz) plt.colorbar(labellog10 幅值)vmin设成vmax减 3 个数量级是为了让 0.1% 以下的弱成分沉底避免噪声把整张图染成雪花。shadingauto在大数据量下能避免渲染卡顿这个细节在处理几万点数据时差别很明显。4. EEMD 与 CEEMDAN模态混叠、端点效应与参数怎么定4.1 模态混叠是怎么被白噪声打散的模态混叠的典型表现是同一个物理分量被拆到两个相邻 IMF 里或者两个时间尺度差别很大的分量挤在同一个 IMF 里。根子在 EMD 的极值点提取机制上间歇分量的存在会让包络在间歇段突然失去极值点筛分过程被迫中断剩下的信息只能塞给下一个 IMF。EEMD 的做法是往原始信号里叠加一组白噪声再分解重复若干次取平均。白噪声在频域上是均匀分布的它把原本断开的极值点区域填满让包络在间歇段也能连续拟合各次分解的随机性在平均之后相互抵消。CEEMDAN 更进一步在每一阶筛分时都加入自适应的噪声并在分解出该阶 IMF 后把噪声残差扣掉重构误差比 EEMD 小一个量级代价是计算量翻倍。两者都不是免费的午餐。EEMD 的平均不能完全消除噪声残留集成次数不够时结果里会带一层薄薄的噪声底CEEMDAN 的参数更敏感epsilon设大了会引入新的伪分量设小了又打不散混叠。from PyEMD import EEMD, CEEMDAN np.random.seed(42) # 固定种子保证结果可复现 eemd EEMD(trials100, noise_width0.2) imfs_eemd eemd.eemd(sig, max_imf6) ceemdan CEEMDAN(trials100, epsilon0.005) imfs_ceemdan ceemdan.ceemdan(sig, max_imf6) print(EEMD 层数:, imfs_eemd.shape[0], CEEMDAN 层数:, imfs_ceemdan.shape[0])noise_width是白噪声标准差相对于信号标准差的倍数取 0.2 意味着加入的噪声比信号本身弱一些。这个值的经验区间是 0.1 到 0.3信号里高频成分占主导时取小值低频趋势明显时取大值。取到 0.5 以上噪声会自己变成新的 IMF分解层数凭空多出一两层。trials是集成次数100 次是工程上的起点噪声残留大约降到 1/√100 即 10% 左右。要做到 1% 需要一万次计算时间不可接受所以实际项目里很少追求完全消除残留而是接受一个小噪声底靠后面的相关性筛选把伪分量剔掉。4.2 集成次数、噪声幅值、种子这三个参数怎么定这三个参数不是独立的调整时要一起看。集成次数主要由噪声残留决定和信号长度基本无关。数据长、单次分解慢的时候宁可减到 50 次也不要缩短数据因为截断数据会引入新的端点效应得不偿失。噪声幅值要跟模态混叠的严重程度对齐。判断方法是先做一次普通 EMD看相邻两个 IMF 的中心频率比值如果接近 1 说明它们本来就该是一个分量这时候噪声幅值往 0.2 以上加如果各阶中心频率呈明显的二倍关系说明分解已经比较干净噪声幅值降到 0.1 以下就够加多了反而把原本清楚的分量搅碎。随机种子必须固定。同一个种子在不同机器上能复现出完全一致的结果换种子之后 IMF 的形状会有细微差别做批量对比实验时如果不固定种子你会以为是数据本身有差异白白排查半天。提示噪声种子、集成次数、噪声幅值、最大 IMF 数这四个值建议做成一个配置字典跟着每个数据文件一起存事后回溯时能省掉大量重复劳动。4.3 分解结果的稳定性验收三次重跑与一致性比对EEMD 和 CEEMDAN 带随机性所以验收的标准不是「结果是否唯一」而是「结果在多次重跑之间是否稳定」。做法是固定种子跑三次每次都换一个种子比较三件事。第一件是 IMF 层数。三次结果层数相差超过 1说明参数设置对数据太敏感需要降低集成次数或者调整噪声幅值再试。层数稳定的前提下再比较每层的中心频率同一序号 IMF 的中心频率漂移超过 5%就意味着该层的物理意义不可靠不能拿来当特征。第二件是能量分布的稳定性。把每层 IMF 的方差占比算出来三组结果的占比曲线应该基本重合如果某个分量的占比在三次之间从 30% 跳到 10%这条分量就不能进后续模型。第三件是重构误差。把全部 IMF 加残差和原始信号相减相对误差应该在 1e-10 量级这是数值精度层面的事一旦超出就说明分解过程本身出了问题比如数据里有 NaN 或者输入维度不对。def stability_check(sig, seeds(1, 2, 3), trials100, width0.2): out [] for sd in seeds: np.random.seed(sd) imfs EEMD(trialstrials, noise_widthwidth).eemd(sig, max_imf6) recon imfs.sum(axis0) out.append((imfs.shape[0], np.var(imfs, axis1) / np.var(sig), np.linalg.norm(recon - sig) / np.linalg.norm(sig))) return out for row in stability_check(sig): print(层数%d 重构相对误差%.2e 能量占比%s % (row[0], row[2], np.round(row[1], 3)))这一步跑完参数基本就定死了。后面所有分析都用这套参数不再逐条数据调整否则批量处理时会陷入「每条数据都要单独调参」的泥潭。5. 避坑与排查HHT 落地最常见的五个翻车现场5.1 分解阶段的三个典型问题现象一不做 EMD直接对原始信号做希尔伯特变换瞬时频率在 0 和奈奎斯特之间乱跳。原因很直接原始信号是多分量叠加解析信号的实部和虚部频谱重叠瞬时频率的物理意义根本不成立。这不是代码写错了是用法越界。解决方式就是老老实实先分解对每个 IMF 单独做变换且只对前几阶高能量 IMF 求瞬时频率后段噪声主导的分量直接跳过。判断该跳过哪一阶用 2.3 里的能量占比和相关系数双阈值比凭感觉切靠谱。现象二同一批数据跑两次IMF 层数和形状都不一样报告没法写。EEMD 和 CEEMDAN 内部有随机噪声不固定种子就必然不可复现普通 EMD 虽然确定但如果FIXE_H没设筛分按标准差阈值动态停止浮点误差累积也会让结果漂移。解决方法是在脚本入口处统一np.random.seed()同时把FIXE_H、nbsym、max_imf全部显式赋值不要依赖库的默认值。血泪经验是默认值换一个库版本就变只有显式写死的参数才是你自己的参数。现象三分解出的 IMF 首尾各有一段幅值异常放大的飞边时频图两端竖直亮带。三次样条在数据端点没有约束包络会向外发散镜像延拓能缓解但不能根治。nbsym从默认的 2 加到 4 或 6 会有改善但更彻底的做法是分析时截掉首尾各 5% 到 10% 的数据只在中间段做统计。如果截断会丢掉关键事件就先把数据向前后各延拓一段比如用 AR 模型预测或者简单的镜像拼接分解完再裁回原长度这个后悔药比事后修图管用。5.2 谱分析阶段的三个典型问题现象四边际谱和 FFT 幅值谱放在一起对比形状相似但幅值差了好几倍怎么调都对不上。边际谱的定义是幅值沿时间的累加再按样本数归一化FFT 幅值谱是各频率分量的幅值两者量纲本来就不同。要对齐得把边际谱除以样本点数同时把 FFT 结果乘以 2 再除以点数并且两边用同一套频率栅格。更常见的错误是频率分辨率没对齐n_bins设了 256 而 FFT 用了 1024 点形状再怎么像也没法逐点比较。做对比实验前先把两边的频率轴打印出来看一眼。现象五时频图上出现大量负频率或者超高频的散点图面全是噪点。负频来自单分量条件不满足超高频来自端点效应和噪声分量的瞬时频率发散。处理顺序是先剔除能量占比低于 0.005 的 IMF再对保留分量做归一化希尔伯特变换最后把频率超过主频 5 倍的样本点置零。三步做完还压不住说明分解本身有问题回到第 4 章重新调噪声幅值。这里没有什么玄学就是逐层过滤。还有一类不容易察觉的问题数据里混进了直流偏置或者工频干扰EMD 会把它们当成第一阶 IMF 剥出来后面的分量全部错位。预处理阶段先做去均值和高通滤波能省掉后面一大半的解释成本。6. 进阶用法把 HHT 接进实际分析流程的验收习惯单次跑通和批量落地之间隔着一段距离。真正把 HHT 用进日常分析核心不是算法本身而是把「参数固化 自动验收」做成固定动作。我的习惯是先给整批数据定一套统一参数然后写一个批量检查函数把每份数据的 IMF 层数、残差能量占比、前两阶负频比例、重构误差四项指标打出来任何一项超阈值就标红单独看。这样一批几十份数据跑下来异常的那几份会自动浮出来不用一条条画图。指标合理区间超阈值时的第一反应IMF 层数与同批数据偏差不超过 1检查max_imf与噪声幅值残差能量占比小于 0.05增大max_imf或筛分次数前两阶负频比例小于 0.02改用归一化希尔伯特变换重构相对误差小于 1e-9查 NaN、查输入维度EEMD 和 CEEMDAN 的选型也有个大致边界数据量小、对重构精度要求高、能接受十几倍计算时间选 CEEMDAN数据量大、只要时频分布的形状、噪声残留可以后面筛掉选 EEMD。两者的分解层数在多数工程数据上差别不大别在这上面反复纠结。最后一个技巧把时频图和边际谱的生成封装成一个函数输入只有信号、采样率和参数字典输出固定为三张图加一份指标 JSON。从那以后我每次拿到新数据都强制走一遍这个封装不再手写绘图代码。手写看着灵活但改一个色标参数忘了同步另一张图报告里两张图的量纲就对不上这种错误在交付前才被发现代价太高。希望帮到你。本文还有配套的精品资源点击获取