简介这是一套面向 HHT希尔伯特-黄变换的 MATLAB 实现与可视化代码聚焦 EMD 经验模态分解和希尔伯特谱分析适合需要处理非线性、非平稳信号的研究人员、工程师及相关专业学生使用。压缩包采用 rar 格式仅含 1 个 m 文件整体大小 1KB核心脚本将数据预处理、EMD 迭代分解、IMF 提取、残余项计算与希尔伯特谱绘制整合在一起并提供原始信号、各 IMF 分量及谱图的直观展示。已有 643 人学习下载。运行该脚本可以清晰看到 EMD 如何通过局部极值与 spline 插值逐层剥离本征模态函数理解 HSA 如何获得瞬时频率与幅度从而为地震波、生物医学信号、经济数据等复杂信号的时频分析搭建可直接复用的实验框架对正在学习 HHT 原理或希望快速搭建分析工具的人来说是一份紧凑、能上手的参考代码。1. 把 EMD、HHT、时频图焊成一条线的可视化思路emd_visu_hht_EMD_emd_visu_这个像函数名一样的长标题其实对应着一条稳定的信号分析链先用经验模态分解EMD把非平稳信号拆成一组频率由高到低的 IMF再对每个 IMF 做 Hilbert 变换并叠出 HHT 时频谱最后把时间域、IMF 列表和频谱结果放到同一幅图上互相印证。我在处理轴承振动、风速突变这类信号时最常回答的问题是“FFT 已经够用了为什么还要 HHT”。答案是瞬时频率比固定窗频谱更贴近物理过程而emd_visu类可视化能让分解质量和频率变化一眼看出问题。接下来这套方案不依赖某份源码是对这条分析链最短、最可复现的落地写法适合做故障诊断、气象数据和生物电信号的工程技术人员。2. EMD 分解原理与数据准备先让 sifting 过程能复现2.1 IMF 的两个硬条件与包络均值EMD 的假设很朴素任意复杂信号都可以看成若干本征模态函数IMFIntrinsic Mode Function叠加。一个序列要被称为 IMF必须同时满足两个条件。第一在整个数据段里极值点个数与过零点个数相等或至多相差 1第二在任意时刻由局部极大值拟合的上包络与局部极小值拟合的下包络的均值为 0。这两个条件比“看起来像正弦波”严格得多它保证了每个 IMF 都能做有意义的 Hilbert 变换。找出 IMF 的过程叫筛分sifting每轮做三件事用三次样条分别连接所有局部极大值和局部极小值形成上下包络计算上下包络的平均值让原始序列减去这个平均包络得到新序列然后重复。这个减法在去除低频骑行分量的同时也在不断修正波形直到上下包络对称到肉眼几乎看不出偏差。工程上常用停止条件 SD 来卡住迭代SD 等于相邻两次筛分结果的能量差归一化值取值在 0.1 到 0.3 之间。SD 设得越小IMF 越光滑、瞬时频率越连续但计算量和对端点振荡的敏感度都会上升。2.2 一个能用的 sifting 最小实现很多封装库把筛分过程封装成了黑盒排错时反而难下手。我习惯先在本地用 NumPy 加 SciPy 把 sifting 循环写出来再替换成 PyEMD这样每一步都能打印包络和均值。下面这段不是完整生产代码但可以直接跑通并验证“每条 IMF 都满足包络均值接近 0”的结论。import numpy as np from scipy.signal import argrelmax, argrelmin from scipy.interpolate import CubicSpline def envelope_mean(x): n len(x) max_idx argrelmax(x, order2)[0] min_idx argrelmin(x, order2)[0] # 极值点不足两个时直接返回原序列避免端点处样条外插失控 if len(max_idx) 2 or len(min_idx) 2: return x.copy(), np.zeros(n) up CubicSpline(max_idx, x[max_idx])(np.arange(n)) low CubicSpline(min_idx, x[min_idx])(np.arange(n)) return x - 0.5 * (up low), up - low def sift_emd(x, max_imfs6, sd0.2): residue x.copy() imfs [] for _ in range(max_imfs): proto residue.copy() for _ in range(200): prev proto.copy() proto, _ envelope_mean(proto) # 用相邻两次能量差判断是否收敛而不是只看包络均值 energy np.sum((prev - proto) ** 2) / (np.sum(prev ** 2) 1e-12) if energy sd: break imfs.append(proto) residue residue - proto if np.std(residue) 1e-8: break return np.array(imfs), residue if __name__ __main__: fs 1000 t np.linspace(0, 1, fs, endpointFalse) # 模拟 3Hz 调幅的 50Hz 载波叠加上 7Hz 正弦 x (1 0.3 * np.sin(2 * np.pi * 3 * t)) * np.sin(2 * np.pi * 50 * t) \ 0.8 * np.sin(2 * np.pi * 7 * t) imfs, r sift_emd(x, max_imfs6, sd0.2) print(IMF 数量:, len(imfs), 残差能量占比:, f{np.mean(r**2) / np.mean(x**2):.2e})代码里的order2用于限制极值点间距避免噪声制造出密集伪极值。envelope_mean返回均值和包络差方便调试时画出上下包络曲线。停止条件取相邻能量差小于sd这比固定迭代次数更合理因为对采样率高的信号可能 30 次就收敛对低频主导的信号却要跑 100 次以上。max_imfs只是硬性上限实际分解会在残差能量趋近零时提前结束。2.3 数据准备单位、长度、采样率如何影响 HHT 后续结果EMD 对数据格式的要求比 FFT 更苛刻。采样率必须写进后续的 Hilbert 变换和频谱绘制里否则瞬时频率算出来是无量纲的弧度步长。更隐蔽的问题是数据长度少于 200 个点时包络样条会非常不稳定我一般要求至少 1000 点也就是 10 个以上完整周期数据过长也会让 EMD 计算负担变大超过 10 万点时应先分段或降采样。单位也值得统一。加速度计输出的是重力加速度 g气象站输出的是米每秒电压传感器输出的是伏特。EMD 对幅值单位不敏感但后续计算边际谱能量时幅值平方意味着单位会变成 g²/Hz 或 (m/s)²/Hz这个单位要写进图标注里否则报告上的 Y 轴只是一个没有物理含义的数字。常见做法是先把信号减去均值再除以标准差做 Z-Score 标准化。标准化会让瞬时频率保持不变但 HHT 谱的幅值会变成相对强度好处是不同测点之间可以直接比较。注意不要对信号做带通滤波后再交给 EMD。滤波会改变极值分布导致分解出的 IMF 与原始物理模态脱节。如果必须去噪只做 50Hz 工频陷波或去趋势且把处理过程记录清楚。3. 用 emd_visu_hht_EMD 绘制 HHT 时频谱3.1 Hilbert 变换与瞬时频率为什么不能直接用 FFT对第 2 章分解出的每条 IMF 做 Hilbert 变换得到一个解析信号就能在任意时刻算出幅值和相位。瞬时频率定义为相位的导数它反映的是“这一刻振荡多快”而不是一段窗口里的平均频率。FFT 把信号当成无限周期信号的正弦叠加对突变型信号会摊平能量短时傅里叶变换STFT用固定窗长窗口短则频率分辨率差窗口长则时间定位差。HHT 的路径不同先筛掉非平稳成分再做相位求导频率和时间都没有被窗函数约束。emd_visu_hht_EMD的核心不是某个特殊算法而是把“瞬时频率会不会出现负值”“边缘有没有频率跳变”这类结果可视化出来。瞬时频率为负值是常见问题原因通常不是 Hilbert 变换出错而是该 IMF 没有真正满足包络均值条件。这时要回头调大 sifting 迭代次数或减小 SD。可视化是把诊断和结果展示合二为一的关键步骤。3.2 最小可运行绘图代码IMF 波形与 HHT 谱同画下面代码读取上一步分解出的imfs把每条 IMF 的时间波形、瞬时幅值包络和时频谱放进同一张图。import numpy as np import matplotlib.pyplot as plt from scipy.signal import hilbert def plot_emd_hht(imfs, fs1000, tNone, freq_bins128): n_imfs len(imfs) n imfs.shape[1] if t is None: t np.arange(n) / fs fig plt.figure(figsize(12, 2 * n_imfs 2)) spec_ax None for i, imf in enumerate(imfs): analytic hilbert(imf) amp np.abs(analytic) phase np.unwrap(np.angle(analytic)) # 用中央差分计算瞬时频率首尾各补一个点 inst_freq np.zeros_like(phase) inst_freq[1:-1] np.diff(phase, 2) / (2 * np.pi / fs) inst_freq[0] inst_freq[1] inst_freq[-1] inst_freq[-2] ax_wave fig.add_subplot(n_imfs 1, 3, i * 3 1) ax_wave.plot(t, imf, lw0.8) ax_wave.plot(t, amp, lw0.6, linestyle--, colorred) ax_wave.set_ylabel(fIMF{i 1}) # 将 (时间, 瞬时频率, 幅值) 统计到二维直方图即 HHT 谱的雏形 freq_ok np.isfinite(inst_freq) (inst_freq 0) H, xedges, yedges np.histogram2d( t[freq_ok], inst_freq[freq_ok], bins[n // 100, freq_bins], weightsamp[freq_ok] ** 2 ) if spec_ax is None: ax_spec fig.add_subplot(n_imfs 1, 3, (i 1) * 3 - 1) spec_ax ax_spec else: ax_spec fig.add_subplot(n_imfs 1, 3, (i 1) * 3 - 1) pcm ax_spec.pcolormesh( xedges, yedges * fs / (2 * np.pi), H.T, shadingauto, cmapmagma ) ax_spec.set_xlim(t[0], t[-1]) ax_spec.set_ylim(0, fs / 2) ax_spec.set_yticks(np.linspace(0, fs / 2, 5).astype(int)) fig.suptitle(emd_visu_hht_EMD: IMF waveforms and HHT spectrum) fig.colorbar(pcm, axspec_ax, labelEnergy (a.u.)) return fig if __name__ __main__: fs 1000 t np.linspace(0, 1, fs, endpointFalse) x (1 0.3 * np.sin(2 * np.pi * 3 * t)) * np.sin(2 * np.pi * 50 * t) \ 0.8 * np.sin(2 * np.pi * 7 * t) imfs, _ sift_emd(x, max_imfs6, sd0.2) fig plot_emd_hht(imfs, fsfs, tt) fig.savefig(emd_hht_spectrum.png, dpi150)这段代码把每条 IMF 的瞬时幅值包络画成红虚线便于观察包络是否平稳。np.histogram2d的作用是把不规则的瞬时频率点聚合成规则网格时间方向以每 10 个样本为一个网格频率方向由freq_bins指定。权重取幅值平方对应能量密度所以 HHT 谱的颜色深浅代表瞬时能量的强弱。yedges * fs / (2 * np.pi)是把“每样本弧度”换算成 Hz 的关键一步漏掉这个换算会让频率轴全部偏小。3.3 频率分辨率与网格参数的取舍HHT 谱的频率分辨率不是固定值它由瞬时频率的离散化步长决定。瞬时频率本身是连续的但直方图会按freq_bins量化。freq_bins偏大时谱图出现大量空洞偏小时会把 49Hz 和 51Hz 混在一起。经验上把freq_bins设为采样率除以 2 再除以最小可分辨带宽。比如采样率 1000Hz目标分辨 2Hz就取 250 个频点。时间方向的网格也有讲究。n // 100表示把 1000 个样本分成 10 段每段有 100 个样本参与统计。对冲击性故障信号时间网格应该更细比如n // 200对平稳旋转机械可以粗到n // 50。网格过细会让同一个瞬时频率点被打散谱图看起来像麻点。另一个常见问题是瞬时频率阵首尾各补了一个值这两个点如果落在 0 频率附近会把谱图底部拉出一条横线。处理办法是绘制时舍去前 5% 和后 5% 的时间区间代价是图谱两侧各少一小段。4. EMD 可视化中的边界条件、模式混叠与排错4.1 端点效应的三种处理方式EMD 的样条包络在数据首尾没有极值约束拟合曲线会大幅外插导致首尾的 IMF 明显摆动HHT 谱图两端出现非物理的高频或负频率。处理端点效应有三种常见做法。第一种是镜像延拓把数据左端向右翻转复制、右端向左翻转复制让极值点变成周期性的延拓长度通常取两端各一个显著极值点的距离。第二种是添加特征波形在首尾各接上一小段与数据自身周期相近的正弦波接缝处用余弦窗加权平滑。第三种最简单直接在 sifting 循环里把首尾各 5% 的样本排除在包络拟合之外只做校验不参与差值计算。我一般优先用镜像延拓因为它的参数最少不会像正弦延拓那样需要估计频率。但镜像延拓要求数据两端斜率不要太大否则镜像出来的极值位置仍不准确。如果 HHT 谱图两端还是上下乱跳就在绘制时主动裁掉前 5% 和后 5%工程报告里注明“边缘截断窗口”即可。没有延拓能彻底消除端点效应这个认知比任何参数都重要。4.2 模式混叠mask 与集合 EMD 的取舍模式混叠是 EMD 最典型的失败模式一段频率接近的分量被拆进两条 IMF或者一条 IMF 里同时装着相差较大的两个频率。直观表现是 IMF 波形间歇性消失又出现HHT 谱上原本该是连续的一条频率带被打断成碎片。混叠的根源是极值点在时间轴上被强信号“淹没”样条包络无法分辨出弱信号的极值。解决模式混叠有两条路线。mask 方法加入一个高频正弦试探信号让强信号的极值密度被重新分割分解完再把试探信号减去集合经验模态分解EEMD 或 CEEMDAN则反复注入白噪声并取平均用它来填满频率缝隙。mask 方法速度快、可解释性强缺点是掩膜幅值和频率要靠试一般取目标高频率分量的 2 到 5 倍EEMD 更自动化但会引入噪声残留且计算量是原始 EMD 的几十倍。在emd_visu_hht_EMD工作流里先用 HHT 谱看碎片位置再决定上哪种方案比直接套用 EEMD 更省算力。4.3 emd_visu 参数对照表与常见误用写代码很容易调参才是 HHT 可视化里真正的成本。下面这张表来自我多次试错后的经验值适用于采样率在 256Hz 到 10kHz 的振动与生物电信号。症状根因调整方向首尾 IMF 幅值突然放大端点效应增加镜像延拓长度或绘图时裁掉首尾 5% 区间瞬时频率出现负值SD 设置过大导致包络不均将 SD 从 0.3 降到 0.1提高最大迭代次数HHT 谱出现横向条纹freq_bins设置过大减小freq_bins使每个频格覆盖更宽带宽高频分量和低频分量挤在同一条 IMF模式混叠对目标频带做 mask EMD或改用 CEEMDAN谱图上能量判断不出来幅值未归一化对输入信号做 Z-Score 标准化后再分解分解出的 IMF 数量超过 8 条将噪声当信号分解设定max_imfsmin(8, log2(n))并检查残差占比另一个常见误用是把 EMD 当滤波器希望靠丢弃前面几条 IMF 来去除高频噪声。实际上前几条 IMF 很可能包含真实物理特征直接丢弃会把短时冲击的尖峰削平。正确做法是保留所有 IMF把 HHT 谱的边带噪声单独做阈值处理例如只显示能量大于最大能量 2% 的区域。注意如果两条 IMF 的瞬时频率曲线几乎重合说明分解过分解或停止条件太松。此时不应人工合并 IMF而是重新设置更小的 SD 值增加 sifting 迭代次数。5. 用正交性与边际谱验证再进入故障特征提取5.1 快速验证 IMF 是否“干净”分解完成后别急着看图说话先跑一段正交性检查。好的 EMD 分解要求任意两条 IMF 之间近似正交也就是内积量级远小于各自能量。同时所有 IMF 加残差应该能还原原始信号还原误差能量占比通常在 1e-6 量级。如果还原误差明显偏大说明 sifting 循环提前退出或某条 IMF 被错误迭代。def check_emd_quality(x, imfs, residue): recon np.sum(imfs, axis0) residue rel_err np.sqrt(np.mean((x - recon) ** 2)) / (np.sqrt(np.mean(x ** 2)) 1e-12) n len(imfs) orth 0.0 for i in range(n): for j in range(i 1, n): # 用 Pearson 相关系数代替内积单位无关且更直观 coef np.corrcoef(imfs[i], imfs[j])[0, 1] orth abs(coef) return rel_err, orth / (n * (n - 1) / 2) rel_err, avg_coef check_emd_quality(x, imfs, r) print(重构误差:, f{rel_err:.1e}, IMF平均相关系数:, f{avg_coef:.3f})这段代码返回两个指标重构误差和平均相关系数。重构误差大于 1e-4 时先检查残差标准差平均相关系数大于 0.1 时大概率存在模式混叠。注意相关系数不等于正交性但在同一数据尺度下它比内积更容易被非技术人员理解。报告里把这两个数字贴出来比贴十张 HHT 谱图更有说服力。5.2 边际谱与 FFT 幅值谱的差别HHT 谱是时间和频率的二维函数把时间方向积分得到的就是边际谱。它表示每个频率点在观测时间内积累的总能量。边际谱与 FFT 谱很像但含义不同FFT 幅值谱描述的是固定频率分量的幅值而边际谱描述的是瞬时频率落在某个区间的概率密度乘能量。对频率随时间变化的信号FFT 会在基频周围展开一段边带边际谱则把能量压缩在一个更窄的频率带内这个压缩特性正是故障诊断里区分松旷和磨损的关键指标。计算边际谱时不要直接对每条 IMF 的频率直方图求和而是对 HHT 谱矩阵做np.sum(H, axis0)然后按频率网格画阶梯图。这个计算思路简单却能暴露一个隐藏陷阱瞬时频率小于 0 的点会被histogram2d丢弃如果丢弃比例超过 5%说明 Hilbert 变换阶段就不干净边际谱会整体少一块能量。5.3 一个有用的双层拆解实验最后一招是把可视化工具变成验证工具。对同一段数据做两次分解第一次用常规sd0.2第二次用更严格的sd0.1并把最大迭代次数翻倍。然后把两次 HHT 谱相减得到差分谱。差分谱上能量集中的位置就是对停止条件最敏感的区域那里通常不够稳定也是现场分析真正需要关注的位置。比如轴承早期故障差分谱会在某个特征频率附近出现白色条带而其他位置接近零。用这个方法emd_visu_hht_EMD就不只是展示工具而是能反推算法参数合理性的诊断装置。本文还有配套的精品资源点击获取 EMD-HHT时频谱可视化实战：从Sifting到HHT谱图绘制
简介这是一套面向 HHT希尔伯特-黄变换的 MATLAB 实现与可视化代码聚焦 EMD 经验模态分解和希尔伯特谱分析适合需要处理非线性、非平稳信号的研究人员、工程师及相关专业学生使用。压缩包采用 rar 格式仅含 1 个 m 文件整体大小 1KB核心脚本将数据预处理、EMD 迭代分解、IMF 提取、残余项计算与希尔伯特谱绘制整合在一起并提供原始信号、各 IMF 分量及谱图的直观展示。已有 643 人学习下载。运行该脚本可以清晰看到 EMD 如何通过局部极值与 spline 插值逐层剥离本征模态函数理解 HSA 如何获得瞬时频率与幅度从而为地震波、生物医学信号、经济数据等复杂信号的时频分析搭建可直接复用的实验框架对正在学习 HHT 原理或希望快速搭建分析工具的人来说是一份紧凑、能上手的参考代码。1. 把 EMD、HHT、时频图焊成一条线的可视化思路emd_visu_hht_EMD_emd_visu_这个像函数名一样的长标题其实对应着一条稳定的信号分析链先用经验模态分解EMD把非平稳信号拆成一组频率由高到低的 IMF再对每个 IMF 做 Hilbert 变换并叠出 HHT 时频谱最后把时间域、IMF 列表和频谱结果放到同一幅图上互相印证。我在处理轴承振动、风速突变这类信号时最常回答的问题是“FFT 已经够用了为什么还要 HHT”。答案是瞬时频率比固定窗频谱更贴近物理过程而emd_visu类可视化能让分解质量和频率变化一眼看出问题。接下来这套方案不依赖某份源码是对这条分析链最短、最可复现的落地写法适合做故障诊断、气象数据和生物电信号的工程技术人员。2. EMD 分解原理与数据准备先让 sifting 过程能复现2.1 IMF 的两个硬条件与包络均值EMD 的假设很朴素任意复杂信号都可以看成若干本征模态函数IMFIntrinsic Mode Function叠加。一个序列要被称为 IMF必须同时满足两个条件。第一在整个数据段里极值点个数与过零点个数相等或至多相差 1第二在任意时刻由局部极大值拟合的上包络与局部极小值拟合的下包络的均值为 0。这两个条件比“看起来像正弦波”严格得多它保证了每个 IMF 都能做有意义的 Hilbert 变换。找出 IMF 的过程叫筛分sifting每轮做三件事用三次样条分别连接所有局部极大值和局部极小值形成上下包络计算上下包络的平均值让原始序列减去这个平均包络得到新序列然后重复。这个减法在去除低频骑行分量的同时也在不断修正波形直到上下包络对称到肉眼几乎看不出偏差。工程上常用停止条件 SD 来卡住迭代SD 等于相邻两次筛分结果的能量差归一化值取值在 0.1 到 0.3 之间。SD 设得越小IMF 越光滑、瞬时频率越连续但计算量和对端点振荡的敏感度都会上升。2.2 一个能用的 sifting 最小实现很多封装库把筛分过程封装成了黑盒排错时反而难下手。我习惯先在本地用 NumPy 加 SciPy 把 sifting 循环写出来再替换成 PyEMD这样每一步都能打印包络和均值。下面这段不是完整生产代码但可以直接跑通并验证“每条 IMF 都满足包络均值接近 0”的结论。import numpy as np from scipy.signal import argrelmax, argrelmin from scipy.interpolate import CubicSpline def envelope_mean(x): n len(x) max_idx argrelmax(x, order2)[0] min_idx argrelmin(x, order2)[0] # 极值点不足两个时直接返回原序列避免端点处样条外插失控 if len(max_idx) 2 or len(min_idx) 2: return x.copy(), np.zeros(n) up CubicSpline(max_idx, x[max_idx])(np.arange(n)) low CubicSpline(min_idx, x[min_idx])(np.arange(n)) return x - 0.5 * (up low), up - low def sift_emd(x, max_imfs6, sd0.2): residue x.copy() imfs [] for _ in range(max_imfs): proto residue.copy() for _ in range(200): prev proto.copy() proto, _ envelope_mean(proto) # 用相邻两次能量差判断是否收敛而不是只看包络均值 energy np.sum((prev - proto) ** 2) / (np.sum(prev ** 2) 1e-12) if energy sd: break imfs.append(proto) residue residue - proto if np.std(residue) 1e-8: break return np.array(imfs), residue if __name__ __main__: fs 1000 t np.linspace(0, 1, fs, endpointFalse) # 模拟 3Hz 调幅的 50Hz 载波叠加上 7Hz 正弦 x (1 0.3 * np.sin(2 * np.pi * 3 * t)) * np.sin(2 * np.pi * 50 * t) \ 0.8 * np.sin(2 * np.pi * 7 * t) imfs, r sift_emd(x, max_imfs6, sd0.2) print(IMF 数量:, len(imfs), 残差能量占比:, f{np.mean(r**2) / np.mean(x**2):.2e})代码里的order2用于限制极值点间距避免噪声制造出密集伪极值。envelope_mean返回均值和包络差方便调试时画出上下包络曲线。停止条件取相邻能量差小于sd这比固定迭代次数更合理因为对采样率高的信号可能 30 次就收敛对低频主导的信号却要跑 100 次以上。max_imfs只是硬性上限实际分解会在残差能量趋近零时提前结束。2.3 数据准备单位、长度、采样率如何影响 HHT 后续结果EMD 对数据格式的要求比 FFT 更苛刻。采样率必须写进后续的 Hilbert 变换和频谱绘制里否则瞬时频率算出来是无量纲的弧度步长。更隐蔽的问题是数据长度少于 200 个点时包络样条会非常不稳定我一般要求至少 1000 点也就是 10 个以上完整周期数据过长也会让 EMD 计算负担变大超过 10 万点时应先分段或降采样。单位也值得统一。加速度计输出的是重力加速度 g气象站输出的是米每秒电压传感器输出的是伏特。EMD 对幅值单位不敏感但后续计算边际谱能量时幅值平方意味着单位会变成 g²/Hz 或 (m/s)²/Hz这个单位要写进图标注里否则报告上的 Y 轴只是一个没有物理含义的数字。常见做法是先把信号减去均值再除以标准差做 Z-Score 标准化。标准化会让瞬时频率保持不变但 HHT 谱的幅值会变成相对强度好处是不同测点之间可以直接比较。注意不要对信号做带通滤波后再交给 EMD。滤波会改变极值分布导致分解出的 IMF 与原始物理模态脱节。如果必须去噪只做 50Hz 工频陷波或去趋势且把处理过程记录清楚。3. 用 emd_visu_hht_EMD 绘制 HHT 时频谱3.1 Hilbert 变换与瞬时频率为什么不能直接用 FFT对第 2 章分解出的每条 IMF 做 Hilbert 变换得到一个解析信号就能在任意时刻算出幅值和相位。瞬时频率定义为相位的导数它反映的是“这一刻振荡多快”而不是一段窗口里的平均频率。FFT 把信号当成无限周期信号的正弦叠加对突变型信号会摊平能量短时傅里叶变换STFT用固定窗长窗口短则频率分辨率差窗口长则时间定位差。HHT 的路径不同先筛掉非平稳成分再做相位求导频率和时间都没有被窗函数约束。emd_visu_hht_EMD的核心不是某个特殊算法而是把“瞬时频率会不会出现负值”“边缘有没有频率跳变”这类结果可视化出来。瞬时频率为负值是常见问题原因通常不是 Hilbert 变换出错而是该 IMF 没有真正满足包络均值条件。这时要回头调大 sifting 迭代次数或减小 SD。可视化是把诊断和结果展示合二为一的关键步骤。3.2 最小可运行绘图代码IMF 波形与 HHT 谱同画下面代码读取上一步分解出的imfs把每条 IMF 的时间波形、瞬时幅值包络和时频谱放进同一张图。import numpy as np import matplotlib.pyplot as plt from scipy.signal import hilbert def plot_emd_hht(imfs, fs1000, tNone, freq_bins128): n_imfs len(imfs) n imfs.shape[1] if t is None: t np.arange(n) / fs fig plt.figure(figsize(12, 2 * n_imfs 2)) spec_ax None for i, imf in enumerate(imfs): analytic hilbert(imf) amp np.abs(analytic) phase np.unwrap(np.angle(analytic)) # 用中央差分计算瞬时频率首尾各补一个点 inst_freq np.zeros_like(phase) inst_freq[1:-1] np.diff(phase, 2) / (2 * np.pi / fs) inst_freq[0] inst_freq[1] inst_freq[-1] inst_freq[-2] ax_wave fig.add_subplot(n_imfs 1, 3, i * 3 1) ax_wave.plot(t, imf, lw0.8) ax_wave.plot(t, amp, lw0.6, linestyle--, colorred) ax_wave.set_ylabel(fIMF{i 1}) # 将 (时间, 瞬时频率, 幅值) 统计到二维直方图即 HHT 谱的雏形 freq_ok np.isfinite(inst_freq) (inst_freq 0) H, xedges, yedges np.histogram2d( t[freq_ok], inst_freq[freq_ok], bins[n // 100, freq_bins], weightsamp[freq_ok] ** 2 ) if spec_ax is None: ax_spec fig.add_subplot(n_imfs 1, 3, (i 1) * 3 - 1) spec_ax ax_spec else: ax_spec fig.add_subplot(n_imfs 1, 3, (i 1) * 3 - 1) pcm ax_spec.pcolormesh( xedges, yedges * fs / (2 * np.pi), H.T, shadingauto, cmapmagma ) ax_spec.set_xlim(t[0], t[-1]) ax_spec.set_ylim(0, fs / 2) ax_spec.set_yticks(np.linspace(0, fs / 2, 5).astype(int)) fig.suptitle(emd_visu_hht_EMD: IMF waveforms and HHT spectrum) fig.colorbar(pcm, axspec_ax, labelEnergy (a.u.)) return fig if __name__ __main__: fs 1000 t np.linspace(0, 1, fs, endpointFalse) x (1 0.3 * np.sin(2 * np.pi * 3 * t)) * np.sin(2 * np.pi * 50 * t) \ 0.8 * np.sin(2 * np.pi * 7 * t) imfs, _ sift_emd(x, max_imfs6, sd0.2) fig plot_emd_hht(imfs, fsfs, tt) fig.savefig(emd_hht_spectrum.png, dpi150)这段代码把每条 IMF 的瞬时幅值包络画成红虚线便于观察包络是否平稳。np.histogram2d的作用是把不规则的瞬时频率点聚合成规则网格时间方向以每 10 个样本为一个网格频率方向由freq_bins指定。权重取幅值平方对应能量密度所以 HHT 谱的颜色深浅代表瞬时能量的强弱。yedges * fs / (2 * np.pi)是把“每样本弧度”换算成 Hz 的关键一步漏掉这个换算会让频率轴全部偏小。3.3 频率分辨率与网格参数的取舍HHT 谱的频率分辨率不是固定值它由瞬时频率的离散化步长决定。瞬时频率本身是连续的但直方图会按freq_bins量化。freq_bins偏大时谱图出现大量空洞偏小时会把 49Hz 和 51Hz 混在一起。经验上把freq_bins设为采样率除以 2 再除以最小可分辨带宽。比如采样率 1000Hz目标分辨 2Hz就取 250 个频点。时间方向的网格也有讲究。n // 100表示把 1000 个样本分成 10 段每段有 100 个样本参与统计。对冲击性故障信号时间网格应该更细比如n // 200对平稳旋转机械可以粗到n // 50。网格过细会让同一个瞬时频率点被打散谱图看起来像麻点。另一个常见问题是瞬时频率阵首尾各补了一个值这两个点如果落在 0 频率附近会把谱图底部拉出一条横线。处理办法是绘制时舍去前 5% 和后 5% 的时间区间代价是图谱两侧各少一小段。4. EMD 可视化中的边界条件、模式混叠与排错4.1 端点效应的三种处理方式EMD 的样条包络在数据首尾没有极值约束拟合曲线会大幅外插导致首尾的 IMF 明显摆动HHT 谱图两端出现非物理的高频或负频率。处理端点效应有三种常见做法。第一种是镜像延拓把数据左端向右翻转复制、右端向左翻转复制让极值点变成周期性的延拓长度通常取两端各一个显著极值点的距离。第二种是添加特征波形在首尾各接上一小段与数据自身周期相近的正弦波接缝处用余弦窗加权平滑。第三种最简单直接在 sifting 循环里把首尾各 5% 的样本排除在包络拟合之外只做校验不参与差值计算。我一般优先用镜像延拓因为它的参数最少不会像正弦延拓那样需要估计频率。但镜像延拓要求数据两端斜率不要太大否则镜像出来的极值位置仍不准确。如果 HHT 谱图两端还是上下乱跳就在绘制时主动裁掉前 5% 和后 5%工程报告里注明“边缘截断窗口”即可。没有延拓能彻底消除端点效应这个认知比任何参数都重要。4.2 模式混叠mask 与集合 EMD 的取舍模式混叠是 EMD 最典型的失败模式一段频率接近的分量被拆进两条 IMF或者一条 IMF 里同时装着相差较大的两个频率。直观表现是 IMF 波形间歇性消失又出现HHT 谱上原本该是连续的一条频率带被打断成碎片。混叠的根源是极值点在时间轴上被强信号“淹没”样条包络无法分辨出弱信号的极值。解决模式混叠有两条路线。mask 方法加入一个高频正弦试探信号让强信号的极值密度被重新分割分解完再把试探信号减去集合经验模态分解EEMD 或 CEEMDAN则反复注入白噪声并取平均用它来填满频率缝隙。mask 方法速度快、可解释性强缺点是掩膜幅值和频率要靠试一般取目标高频率分量的 2 到 5 倍EEMD 更自动化但会引入噪声残留且计算量是原始 EMD 的几十倍。在emd_visu_hht_EMD工作流里先用 HHT 谱看碎片位置再决定上哪种方案比直接套用 EEMD 更省算力。4.3 emd_visu 参数对照表与常见误用写代码很容易调参才是 HHT 可视化里真正的成本。下面这张表来自我多次试错后的经验值适用于采样率在 256Hz 到 10kHz 的振动与生物电信号。症状根因调整方向首尾 IMF 幅值突然放大端点效应增加镜像延拓长度或绘图时裁掉首尾 5% 区间瞬时频率出现负值SD 设置过大导致包络不均将 SD 从 0.3 降到 0.1提高最大迭代次数HHT 谱出现横向条纹freq_bins设置过大减小freq_bins使每个频格覆盖更宽带宽高频分量和低频分量挤在同一条 IMF模式混叠对目标频带做 mask EMD或改用 CEEMDAN谱图上能量判断不出来幅值未归一化对输入信号做 Z-Score 标准化后再分解分解出的 IMF 数量超过 8 条将噪声当信号分解设定max_imfsmin(8, log2(n))并检查残差占比另一个常见误用是把 EMD 当滤波器希望靠丢弃前面几条 IMF 来去除高频噪声。实际上前几条 IMF 很可能包含真实物理特征直接丢弃会把短时冲击的尖峰削平。正确做法是保留所有 IMF把 HHT 谱的边带噪声单独做阈值处理例如只显示能量大于最大能量 2% 的区域。注意如果两条 IMF 的瞬时频率曲线几乎重合说明分解过分解或停止条件太松。此时不应人工合并 IMF而是重新设置更小的 SD 值增加 sifting 迭代次数。5. 用正交性与边际谱验证再进入故障特征提取5.1 快速验证 IMF 是否“干净”分解完成后别急着看图说话先跑一段正交性检查。好的 EMD 分解要求任意两条 IMF 之间近似正交也就是内积量级远小于各自能量。同时所有 IMF 加残差应该能还原原始信号还原误差能量占比通常在 1e-6 量级。如果还原误差明显偏大说明 sifting 循环提前退出或某条 IMF 被错误迭代。def check_emd_quality(x, imfs, residue): recon np.sum(imfs, axis0) residue rel_err np.sqrt(np.mean((x - recon) ** 2)) / (np.sqrt(np.mean(x ** 2)) 1e-12) n len(imfs) orth 0.0 for i in range(n): for j in range(i 1, n): # 用 Pearson 相关系数代替内积单位无关且更直观 coef np.corrcoef(imfs[i], imfs[j])[0, 1] orth abs(coef) return rel_err, orth / (n * (n - 1) / 2) rel_err, avg_coef check_emd_quality(x, imfs, r) print(重构误差:, f{rel_err:.1e}, IMF平均相关系数:, f{avg_coef:.3f})这段代码返回两个指标重构误差和平均相关系数。重构误差大于 1e-4 时先检查残差标准差平均相关系数大于 0.1 时大概率存在模式混叠。注意相关系数不等于正交性但在同一数据尺度下它比内积更容易被非技术人员理解。报告里把这两个数字贴出来比贴十张 HHT 谱图更有说服力。5.2 边际谱与 FFT 幅值谱的差别HHT 谱是时间和频率的二维函数把时间方向积分得到的就是边际谱。它表示每个频率点在观测时间内积累的总能量。边际谱与 FFT 谱很像但含义不同FFT 幅值谱描述的是固定频率分量的幅值而边际谱描述的是瞬时频率落在某个区间的概率密度乘能量。对频率随时间变化的信号FFT 会在基频周围展开一段边带边际谱则把能量压缩在一个更窄的频率带内这个压缩特性正是故障诊断里区分松旷和磨损的关键指标。计算边际谱时不要直接对每条 IMF 的频率直方图求和而是对 HHT 谱矩阵做np.sum(H, axis0)然后按频率网格画阶梯图。这个计算思路简单却能暴露一个隐藏陷阱瞬时频率小于 0 的点会被histogram2d丢弃如果丢弃比例超过 5%说明 Hilbert 变换阶段就不干净边际谱会整体少一块能量。5.3 一个有用的双层拆解实验最后一招是把可视化工具变成验证工具。对同一段数据做两次分解第一次用常规sd0.2第二次用更严格的sd0.1并把最大迭代次数翻倍。然后把两次 HHT 谱相减得到差分谱。差分谱上能量集中的位置就是对停止条件最敏感的区域那里通常不够稳定也是现场分析真正需要关注的位置。比如轴承早期故障差分谱会在某个特征频率附近出现白色条带而其他位置接近零。用这个方法emd_visu_hht_EMD就不只是展示工具而是能反推算法参数合理性的诊断装置。本文还有配套的精品资源点击获取