简介希尔伯特-黄变换HHT时频图是分析非线性、非平稳信号的重要工具。这份压缩包提供一段MATLAB实现代码面向信号处理方向的学生、科研人员与工程开发者可用于快速理解HHT时频图的完整生成流程并迁移至生物医学信号、机械故障诊断、地震数据分析等应用场景。包体十分精简仅含1个m文件大小约1KB代码短小但覆盖核心步骤经验模态分解EMD得到IMF分量希尔伯特变换计算瞬时频率与幅度再组合绘制成时频分布图便于对照理论逐一学习。已有441人学习下载。通过运行并修改这段代码可直观观察不同参数和IMF分量对时频图的影响掌握非线性非平稳信号分析的实现要点避开HHT应用中的常见误区为后续深入研究和工程实践打下基础。1. HHT时频图的坑从“毛刺满屏”说起从网上下载一个写着“HHT时频图”的压缩包解压跑完得到的往往不是论文里那种干净利落的频率脊线而是一张布满毛刺的彩色噪点图。这不是你下载错了代码而是HHT时频图本身对分解参数和绘制方式极度敏感。HHTHilbert-Huang Transform由EMD分解和Hilbert变换两步组成前一步决定频率成分分得干不干净后一步决定画出来的瞬时频率有没有物理意义。这篇文章不堆公式讲怎么把一张能用的HHT时频图画出来以及当它画坏时从哪个参数开始调。适合做机械故障诊断、地震信号分析、脑电或振动信号处理的工程师。2. 手写HHT时频图EMD分解与Hilbert变换的最小实现2.1 为什么HHT时频图比STFT更“锐”短时傅里叶变换STFT受窗函数限制时间分辨率和频率分辨率不能同时提高。HHT的思路完全不同先用EMD把信号分解成若干个本征模态函数IMF每个IMF是单分量信号然后对每个IMF做Hilbert变换求瞬时频率。瞬时频率是逐点定义的所以时频图理论上可以达到任意时间分辨率频率也随信号自适应变化。代价是EMD的自适应性带来了不确定性分解份数、端点效应、模态混叠全部会反映在最终的时频图上。“HT”部分的常见做法是对IMF做Hilbert变换得到解析信号再对相位做差分得到瞬时频率。这里有两个坑一是相位差分不是唯一的瞬时频率定义二是差分运算会放大噪声。后面第3章会专门处理。先跑通一幅不加花哨修饰的时频图再说。2.2 最小实现从合成信号到一张可用的HHT时频图2.2.1 生成一个带时变频率的测试信号先用一段调频加调幅的合成信号做演示频率从20Hz线性扫到50Hz叠加一个120Hz的低幅正弦作为干扰。这样可以清晰看出时频图上哪些成分是信号本身的哪些是分解产生的假象。import numpy as np from scipy.signal import hilbert from PyEMD import EMD import matplotlib.pyplot as plt fs 1000 t np.linspace(0, 1, fs, endpointFalse) # 20Hz - 50Hz 线性扫频幅值同时做慢调制 f_inst 20 30 * t phase 2 * np.pi * np.cumsum(f_inst) / fs x np.sin(phase) * (1 0.3 * np.sin(2 * np.pi * 5 * t)) x x 0.1 * np.sin(2 * np.pi * 120 * t) # 高频干扰EMD对象需要时间轴参数t建议用等间隔浮点数组不要用整数索引。PyEMD内部会依据t的间隔做样条插值时间轴不均匀会直接影响包络拟合质量。2.2.2 EMD分解与IMF提取emd EMD() imfs emd.emd(x, t) # 形状 (n, len(t))最后一行是残差 n_imfs imfs.shape[0] imf_list imfs[:-1] # 去掉残差趋势项 residue imfs[-1]emd.emd()是标准EMD流程内部默认用三次样条包络sifting次数由SD准则控制。PyEMD默认的SD_THRESHOLD是0.1容忍度越低分解出的IMF数量越多耗时也越长。对长度为1000点的信号n_imfs一般在5到8之间如果超过10说明信号噪声偏大或SD阈值设得太严。分解后立即检查两点每个IMF的过零点数和极值点数是否相等或差1残差是否单调。这两条是IMF定义的基本判据任何一条不满足后面画出的时频图都会出现不连续跳变。2.2.3 用Hilbert变换计算瞬时频率再统计成时频图def inst_freq(imf, t): h hilbert(imf) phase np.unwrap(np.angle(h)) freq np.diff(phase) / (2 * np.pi * np.diff(t)) amp np.abs(h)[:-1] return t[:-1], freq, amp freq_bins np.linspace(0, 200, 400) # 频率轴覆盖到200Hz time_bins t energy np.zeros((len(freq_bins) - 1, len(time_bins) - 1)) for imf in imf_list: tt, ff, aa inst_freq(imf, t) keep (ff 0) (ff 200) energy np.histogram2d(tt[keep], ff[keep], bins[time_bins, freq_bins], weightsaa[keep] ** 2)[0].T这段代码做的事情可以拆开看hilbert(imf)返回解析信号实部是原IMF虚部是它的Hilbert变换np.angle取相位角np.unwrap把相位从(-π, π]展开成连续曲线避免在±π处产生2π跳变。瞬时频率是相位的一阶导数。histogram2d按时间、频率二维分箱用瞬时幅值的平方作为权重和“能量”对应。最后的.T转置是把形状从(时间, 频率)转成(频率, 时间)方便pcolormesh直接画。提示瞬时频率差分后长度比原信号少1所以t[:-1]和amp[:-1]必须对齐。习惯性把三者压缩到同一长度能避免后面画图时坐标错位。2.2.4 绘制时频图fig, ax plt.subplots(figsize(10, 5)) mesh ax.pcolormesh(time_bins, freq_bins, energy, shadingauto, cmapjet) ax.set_ylim(0, 150) ax.set_xlabel(Time (s)) ax.set_ylabel(Frequency (Hz)) fig.colorbar(mesh, axax, labelEnergy) plt.tight_layout() plt.savefig(hht_spectrum.png, dpi300)用shadingauto可以消除pcolormesh对网格边界的抱怨。频率上限先放到150Hz120Hz干扰那一条可以看见下限不硬截留给后续端点处理。到这里一张最原始的HHT时频图已经能跑出来了但大概率你会看到三个问题左右两端出现几条斜冲上去的竖线120Hz那条横线抖成波浪低频区域有一片连续色块。这三个问题对应第3章的三类修正。3. 画好HHT时频图的参数修正表端点效应、模态混叠与负频率3.1 端点发散时频图两端的“飞翼”与镜像延拓时频图最显眼的毛病是两端出现向上或向下的频率飞翼频率值瞬间冲到几百赫兹颜色异常亮。原因是样条包络在信号首尾处没有足够极值点约束包络线会过冲导致IMF在端点处产生形变。这类形变在做Hilbert变换后直接反映为瞬时频率的剧烈波动。常见做法是给信号做镜像延拓以首尾极值点所在时刻为镜面把信号往外翻一段让样条包络在端点处有数据可依。PyEMD的EMD类已经内置了镜像延拓选项通过extrema_detection和NEW_MIRROR参数控制。实测中只延拓两端各一个周期的数据就能消除大部分飞翼。emd_ext EMD(extrema_detectionparabol) emd_ext.NEW_MIRROR True imfs_ext emd_ext.emd(x, t)extrema_detection可选项有simple和parabol后者用抛物线拟合极值点的位置和取值对噪声更鲁棒但也更慢。如果在你的数据上两种方式结果差别不大保留simple即可。判断标准飞翼如果只出现在全长时间轴的前后2%区域内可以接受如果侵入到中间区域就需要改延拓或增加端点处的数据长度。3.2 模态混叠频率脊线十字交叉与EEMD替代模态混叠的典型特征是某一条IMF的瞬时频率跳到另一条IMF的频率范围内时频图上出现“十字交叉”或两条脊线互相穿插。成因是信号中包含间歇性高频分量比如第2.2节里120Hz那个分量在EMD筛分过程中会干扰低频IMF的包络极值点分布。处理模态混叠最直接的手段是换用EEMD集合经验模态分解或CEEMDAN。EEMD给信号加多次白噪声利用噪声的统计特性把不同尺度的分量分离到不同IMF中然后对多次结果取平均。from PyEMD import EEMD eemd EEMD() imfs_ee eemd.eemd(x, t, trials50, noise_width0.05)trials是集成次数noise_width是噪声幅值相对于信号标准差的比值。noise_width取0.05到0.2之间比较安全太小起不到抑制混叠的作用太大会把噪声带进IMF。集成次数50次以上时多次平均的残差影响可以忽略。代价是计算时间线性增长10000点以上的信号建议先把trials调到30做试验确认分解质量后再拉满。注意问题出在“间歇性高频”时EEMD效果显著但如果你只有纯线性扫频信号普通EMD就够了。不要为了“听起来更高级”强行加噪声那会让本来干净的瞬时频率曲线出现随机抖动。3.3 负频率与低频鼓包相位解缠的副作用用np.diff(np.unwrap(np.angle(h)))算瞬时频率有一个数学上的副作用相位差分会放大高频噪声噪声点处的相位抖动会把瞬时频率压成负值或推向极高值。表现在时频图上就是低频区域出现整片色块或者频率轴上出现垂直的亮条纹。处理这个问题的优先顺序是先对Hilbert变换前的IMF做轻平滑比如3点中值滤波去除孤立极值点。再对瞬时频率做移动平均窗口长度取信号周期的1/10左右。最后做频率掩膜只保留全时间轴能量占比前95%的频率范围。from scipy.ndimage import median_filter freq_smooth median_filter(freq, size3)中值滤波对尖峰脉冲的抑制能力比均值滤波强且不会模糊真实频率跳变点。如果你的信号本身存在频率突变比如转速阶跃不要在瞬时频率上做时间窗过长的平滑否则阶跃会被拉成斜坡。3.4 已踩过的参数整理成速查表症状优先调整参数位置失效后的备选方案两端飞翼镜像延拓 截断边界NEW_MIRRORTrue绘图set_xlim信号两端补一段衰减窗脊线交叉提升sifting容忍度SD_THRESHOLD0.2换EEMDtrials50低频糊成一片中值滤波瞬时频率median_filter(size3)对能量图做高斯模糊高频亮线降低色标上限绘图的vmax频率轴截断到关注区间整图布满颗粒检查噪声幅值数据预处理带通滤波EEMD的noise_width调低表格里每一项改动都不应该一次全做。我一般一次只动一个参数画两次图对比否则改了六个参数之后检出问题根本不知道是哪一步修复了它。如果你时间紧优先改飞翼和混叠两项其余用绘图层面的掩盖手段处理更省事。4. 大样本下HHT时频图的显示优化频率截断、归一化与下采样4.1 频率截断不要让能量全堆在低频色块里实测数据里能量往往集中在低频段直接画出图像时关注频段的颜色会被整体“压暗”。先把频率轴截到目标区间比如只看0到1kHz再画图视觉对比度会立刻改善。不要依赖set_ylim去做这件事它会保留全部计算量只改变显示范围。f_max 1000 mask freq_bins[:-1] f_max energy_show energy[mask, :] freq_show freq_bins[:-1][mask]这段代码计算对象不变只是取出需要的频率带。如果你的数据本身是高频振动信号目标频段可能只有总量的20%截断后还能顺便把分布稀疏的高频噪声从色标范围中剔除。4.2 imshow重采样把时频图从“能看”变成“能发”pcolormesh在点数少于10万时表现没问题但当信号长度达到几十万点、频率分箱超过1000个时绘图开销会明显增大生成的PDF也会卡顿。这时可以换成imshow它按像素网格重采样内存占用低一个量级。from matplotlib.colors import LogNorm fig, ax plt.subplots(figsize(10, 5)) img ax.imshow(energy_show, aspectauto, originlower, extent[t[0], t[-1], freq_show[0], freq_show[-1]], normLogNorm(vmin1e-3, vmaxenergy.max()), interpolationbilinear, cmapjet) fig.colorbar(img, axax, labelLog Energy)extent参数定义了图像的四边分别在数据坐标里的位置aspectauto让x和y方向独立缩放不会因为频率轴单位不同而压扁波形。LogNorm是给能量图配的对数色标适合动态范围跨了几个数量级的场景。vmin不要设为0对数坐标对0无定义这里用1e-3做下限。提示imshow默认把数组当图像逐像素显示不做任何聚合。如果你不希望看到计算噪声打开interpolationbilinear做平滑别用nearest。4.3 综合显示函数归一化、透明背景与矢量导出把上面几个动作收拢成一个函数方便直接换数据调用def plot_hht(t, imfs, f_max1000, vmin1e-3, cmapjet): energy build_spectrum(t, imfs) # 复用 2.2.3 的统计逻辑 mask freq_bins[:-1] f_max energy_show energy[mask, :] fig, ax plt.subplots(figsize(12, 5)) img ax.imshow(energy_show, aspectauto, originlower, extent[t[0], t[-1], freq_bins[:-1][mask][0], f_max], normLogNorm(vminvmin, vmaxenergy.max()), cmapcmap) ax.set_xlim(t[0], t[-1]) ax.set_ylabel(Frequency (Hz)) fig.colorbar(img, axax, labelEnergy) return fig, ax fig, ax plot_hht(t, imf_list, f_max1000, cmapturbo) fig.savefig(hht_optimized.pdf, bbox_inchestight)矢量图导出用PDF格式正文插图用PNG。PDF保留全部细节适合写报告和论文PNG按bbox_inchestight裁剪白边适合直接嵌入网页。此时你会注意到数据里的主要频率脊线已经清晰可见但颜色最亮的位置到底是什么频率这需要第5章的边际谱来验证而不是拿鼠标在图上猜。5. 用边际谱校验HHT时频图分解质量的自检方法5.1 归一化边际谱与瞬时频率自检时频图好看不等于分解正确。要验证HHT时频图反映的是真实物理频率标准做法是算边际谱把所有时刻的瞬时频率做直方图统计并按能量加权。边际谱的物理含义是“每个频率上信号累计消耗的能量”和傅里叶幅度谱类似但不要求线性平稳。def marginal_spectrum(imfs, t): freq_cum [] amp_cum [] for imf in imfs: h hilbert(imf) phase np.unwrap(np.angle(h)) f np.diff(phase) / (2 * np.pi * np.diff(t)) a np.abs(h)[:-1] freq_cum.append(f) amp_cum.append(a) freq_all np.concatenate(freq_cum) amp_all np.concatenate(amp_cum) return freq_all, amp_all freq_all, amp_all marginal_spectrum(imf_list, t) hist, bins np.histogram(freq_all, bins200, range(0, 150), weightsamp_all**2)对第2.2节的合成信号边际谱的峰值应该出现在20Hz到50Hz范围的两端附近120Hz处有一个小尖峰。因为扫频信号在每个频率上停留时间相同理想情况是20Hz和50Hz两端能量密度最高。如果边际谱峰值完全偏离真实频率说明EMD把信号拆成了没有物理意义的纯数学分量。5.2 用一组阈值判断分解质量检查项通过阈值处理动作边际谱主峰频率与理论中心频率误差5%通过保留当前参数120Hz分量峰与基波峰能量比干扰峰能量低于主峰不满足则调noise_width或trials瞬时频率中位数与理论扫频范围重叠中位数在25-45Hz区间不满足则回查IMF判据负频率占比1%大于1%则加大中值滤波窗口这个表格不是拍脑袋定的它对应的是“分解结果偏离物理事实”的最小可接受范围。真实数据里没有理论频率可对比此时检查第三项瞬时频率的中位数应该在信号的已知频带内如果连中位数都跑到不可信区间最稳妥的做法是回到第2步降低SD阈值重新分解并检查是否该换EEMD。测完边际谱后把marginal_spectrum的结果存成CSV和时频图放在一起。下次再改任何参数先对比边际谱峰值偏移量再决定是否保留改动。这套自检流程跑通后HHT时频图的问题就只剩下审美了。本文还有配套的精品资源点击获取 HHT时频图绘制与优化：从EMD到边际谱的完整实践指南
简介希尔伯特-黄变换HHT时频图是分析非线性、非平稳信号的重要工具。这份压缩包提供一段MATLAB实现代码面向信号处理方向的学生、科研人员与工程开发者可用于快速理解HHT时频图的完整生成流程并迁移至生物医学信号、机械故障诊断、地震数据分析等应用场景。包体十分精简仅含1个m文件大小约1KB代码短小但覆盖核心步骤经验模态分解EMD得到IMF分量希尔伯特变换计算瞬时频率与幅度再组合绘制成时频分布图便于对照理论逐一学习。已有441人学习下载。通过运行并修改这段代码可直观观察不同参数和IMF分量对时频图的影响掌握非线性非平稳信号分析的实现要点避开HHT应用中的常见误区为后续深入研究和工程实践打下基础。1. HHT时频图的坑从“毛刺满屏”说起从网上下载一个写着“HHT时频图”的压缩包解压跑完得到的往往不是论文里那种干净利落的频率脊线而是一张布满毛刺的彩色噪点图。这不是你下载错了代码而是HHT时频图本身对分解参数和绘制方式极度敏感。HHTHilbert-Huang Transform由EMD分解和Hilbert变换两步组成前一步决定频率成分分得干不干净后一步决定画出来的瞬时频率有没有物理意义。这篇文章不堆公式讲怎么把一张能用的HHT时频图画出来以及当它画坏时从哪个参数开始调。适合做机械故障诊断、地震信号分析、脑电或振动信号处理的工程师。2. 手写HHT时频图EMD分解与Hilbert变换的最小实现2.1 为什么HHT时频图比STFT更“锐”短时傅里叶变换STFT受窗函数限制时间分辨率和频率分辨率不能同时提高。HHT的思路完全不同先用EMD把信号分解成若干个本征模态函数IMF每个IMF是单分量信号然后对每个IMF做Hilbert变换求瞬时频率。瞬时频率是逐点定义的所以时频图理论上可以达到任意时间分辨率频率也随信号自适应变化。代价是EMD的自适应性带来了不确定性分解份数、端点效应、模态混叠全部会反映在最终的时频图上。“HT”部分的常见做法是对IMF做Hilbert变换得到解析信号再对相位做差分得到瞬时频率。这里有两个坑一是相位差分不是唯一的瞬时频率定义二是差分运算会放大噪声。后面第3章会专门处理。先跑通一幅不加花哨修饰的时频图再说。2.2 最小实现从合成信号到一张可用的HHT时频图2.2.1 生成一个带时变频率的测试信号先用一段调频加调幅的合成信号做演示频率从20Hz线性扫到50Hz叠加一个120Hz的低幅正弦作为干扰。这样可以清晰看出时频图上哪些成分是信号本身的哪些是分解产生的假象。import numpy as np from scipy.signal import hilbert from PyEMD import EMD import matplotlib.pyplot as plt fs 1000 t np.linspace(0, 1, fs, endpointFalse) # 20Hz - 50Hz 线性扫频幅值同时做慢调制 f_inst 20 30 * t phase 2 * np.pi * np.cumsum(f_inst) / fs x np.sin(phase) * (1 0.3 * np.sin(2 * np.pi * 5 * t)) x x 0.1 * np.sin(2 * np.pi * 120 * t) # 高频干扰EMD对象需要时间轴参数t建议用等间隔浮点数组不要用整数索引。PyEMD内部会依据t的间隔做样条插值时间轴不均匀会直接影响包络拟合质量。2.2.2 EMD分解与IMF提取emd EMD() imfs emd.emd(x, t) # 形状 (n, len(t))最后一行是残差 n_imfs imfs.shape[0] imf_list imfs[:-1] # 去掉残差趋势项 residue imfs[-1]emd.emd()是标准EMD流程内部默认用三次样条包络sifting次数由SD准则控制。PyEMD默认的SD_THRESHOLD是0.1容忍度越低分解出的IMF数量越多耗时也越长。对长度为1000点的信号n_imfs一般在5到8之间如果超过10说明信号噪声偏大或SD阈值设得太严。分解后立即检查两点每个IMF的过零点数和极值点数是否相等或差1残差是否单调。这两条是IMF定义的基本判据任何一条不满足后面画出的时频图都会出现不连续跳变。2.2.3 用Hilbert变换计算瞬时频率再统计成时频图def inst_freq(imf, t): h hilbert(imf) phase np.unwrap(np.angle(h)) freq np.diff(phase) / (2 * np.pi * np.diff(t)) amp np.abs(h)[:-1] return t[:-1], freq, amp freq_bins np.linspace(0, 200, 400) # 频率轴覆盖到200Hz time_bins t energy np.zeros((len(freq_bins) - 1, len(time_bins) - 1)) for imf in imf_list: tt, ff, aa inst_freq(imf, t) keep (ff 0) (ff 200) energy np.histogram2d(tt[keep], ff[keep], bins[time_bins, freq_bins], weightsaa[keep] ** 2)[0].T这段代码做的事情可以拆开看hilbert(imf)返回解析信号实部是原IMF虚部是它的Hilbert变换np.angle取相位角np.unwrap把相位从(-π, π]展开成连续曲线避免在±π处产生2π跳变。瞬时频率是相位的一阶导数。histogram2d按时间、频率二维分箱用瞬时幅值的平方作为权重和“能量”对应。最后的.T转置是把形状从(时间, 频率)转成(频率, 时间)方便pcolormesh直接画。提示瞬时频率差分后长度比原信号少1所以t[:-1]和amp[:-1]必须对齐。习惯性把三者压缩到同一长度能避免后面画图时坐标错位。2.2.4 绘制时频图fig, ax plt.subplots(figsize(10, 5)) mesh ax.pcolormesh(time_bins, freq_bins, energy, shadingauto, cmapjet) ax.set_ylim(0, 150) ax.set_xlabel(Time (s)) ax.set_ylabel(Frequency (Hz)) fig.colorbar(mesh, axax, labelEnergy) plt.tight_layout() plt.savefig(hht_spectrum.png, dpi300)用shadingauto可以消除pcolormesh对网格边界的抱怨。频率上限先放到150Hz120Hz干扰那一条可以看见下限不硬截留给后续端点处理。到这里一张最原始的HHT时频图已经能跑出来了但大概率你会看到三个问题左右两端出现几条斜冲上去的竖线120Hz那条横线抖成波浪低频区域有一片连续色块。这三个问题对应第3章的三类修正。3. 画好HHT时频图的参数修正表端点效应、模态混叠与负频率3.1 端点发散时频图两端的“飞翼”与镜像延拓时频图最显眼的毛病是两端出现向上或向下的频率飞翼频率值瞬间冲到几百赫兹颜色异常亮。原因是样条包络在信号首尾处没有足够极值点约束包络线会过冲导致IMF在端点处产生形变。这类形变在做Hilbert变换后直接反映为瞬时频率的剧烈波动。常见做法是给信号做镜像延拓以首尾极值点所在时刻为镜面把信号往外翻一段让样条包络在端点处有数据可依。PyEMD的EMD类已经内置了镜像延拓选项通过extrema_detection和NEW_MIRROR参数控制。实测中只延拓两端各一个周期的数据就能消除大部分飞翼。emd_ext EMD(extrema_detectionparabol) emd_ext.NEW_MIRROR True imfs_ext emd_ext.emd(x, t)extrema_detection可选项有simple和parabol后者用抛物线拟合极值点的位置和取值对噪声更鲁棒但也更慢。如果在你的数据上两种方式结果差别不大保留simple即可。判断标准飞翼如果只出现在全长时间轴的前后2%区域内可以接受如果侵入到中间区域就需要改延拓或增加端点处的数据长度。3.2 模态混叠频率脊线十字交叉与EEMD替代模态混叠的典型特征是某一条IMF的瞬时频率跳到另一条IMF的频率范围内时频图上出现“十字交叉”或两条脊线互相穿插。成因是信号中包含间歇性高频分量比如第2.2节里120Hz那个分量在EMD筛分过程中会干扰低频IMF的包络极值点分布。处理模态混叠最直接的手段是换用EEMD集合经验模态分解或CEEMDAN。EEMD给信号加多次白噪声利用噪声的统计特性把不同尺度的分量分离到不同IMF中然后对多次结果取平均。from PyEMD import EEMD eemd EEMD() imfs_ee eemd.eemd(x, t, trials50, noise_width0.05)trials是集成次数noise_width是噪声幅值相对于信号标准差的比值。noise_width取0.05到0.2之间比较安全太小起不到抑制混叠的作用太大会把噪声带进IMF。集成次数50次以上时多次平均的残差影响可以忽略。代价是计算时间线性增长10000点以上的信号建议先把trials调到30做试验确认分解质量后再拉满。注意问题出在“间歇性高频”时EEMD效果显著但如果你只有纯线性扫频信号普通EMD就够了。不要为了“听起来更高级”强行加噪声那会让本来干净的瞬时频率曲线出现随机抖动。3.3 负频率与低频鼓包相位解缠的副作用用np.diff(np.unwrap(np.angle(h)))算瞬时频率有一个数学上的副作用相位差分会放大高频噪声噪声点处的相位抖动会把瞬时频率压成负值或推向极高值。表现在时频图上就是低频区域出现整片色块或者频率轴上出现垂直的亮条纹。处理这个问题的优先顺序是先对Hilbert变换前的IMF做轻平滑比如3点中值滤波去除孤立极值点。再对瞬时频率做移动平均窗口长度取信号周期的1/10左右。最后做频率掩膜只保留全时间轴能量占比前95%的频率范围。from scipy.ndimage import median_filter freq_smooth median_filter(freq, size3)中值滤波对尖峰脉冲的抑制能力比均值滤波强且不会模糊真实频率跳变点。如果你的信号本身存在频率突变比如转速阶跃不要在瞬时频率上做时间窗过长的平滑否则阶跃会被拉成斜坡。3.4 已踩过的参数整理成速查表症状优先调整参数位置失效后的备选方案两端飞翼镜像延拓 截断边界NEW_MIRRORTrue绘图set_xlim信号两端补一段衰减窗脊线交叉提升sifting容忍度SD_THRESHOLD0.2换EEMDtrials50低频糊成一片中值滤波瞬时频率median_filter(size3)对能量图做高斯模糊高频亮线降低色标上限绘图的vmax频率轴截断到关注区间整图布满颗粒检查噪声幅值数据预处理带通滤波EEMD的noise_width调低表格里每一项改动都不应该一次全做。我一般一次只动一个参数画两次图对比否则改了六个参数之后检出问题根本不知道是哪一步修复了它。如果你时间紧优先改飞翼和混叠两项其余用绘图层面的掩盖手段处理更省事。4. 大样本下HHT时频图的显示优化频率截断、归一化与下采样4.1 频率截断不要让能量全堆在低频色块里实测数据里能量往往集中在低频段直接画出图像时关注频段的颜色会被整体“压暗”。先把频率轴截到目标区间比如只看0到1kHz再画图视觉对比度会立刻改善。不要依赖set_ylim去做这件事它会保留全部计算量只改变显示范围。f_max 1000 mask freq_bins[:-1] f_max energy_show energy[mask, :] freq_show freq_bins[:-1][mask]这段代码计算对象不变只是取出需要的频率带。如果你的数据本身是高频振动信号目标频段可能只有总量的20%截断后还能顺便把分布稀疏的高频噪声从色标范围中剔除。4.2 imshow重采样把时频图从“能看”变成“能发”pcolormesh在点数少于10万时表现没问题但当信号长度达到几十万点、频率分箱超过1000个时绘图开销会明显增大生成的PDF也会卡顿。这时可以换成imshow它按像素网格重采样内存占用低一个量级。from matplotlib.colors import LogNorm fig, ax plt.subplots(figsize(10, 5)) img ax.imshow(energy_show, aspectauto, originlower, extent[t[0], t[-1], freq_show[0], freq_show[-1]], normLogNorm(vmin1e-3, vmaxenergy.max()), interpolationbilinear, cmapjet) fig.colorbar(img, axax, labelLog Energy)extent参数定义了图像的四边分别在数据坐标里的位置aspectauto让x和y方向独立缩放不会因为频率轴单位不同而压扁波形。LogNorm是给能量图配的对数色标适合动态范围跨了几个数量级的场景。vmin不要设为0对数坐标对0无定义这里用1e-3做下限。提示imshow默认把数组当图像逐像素显示不做任何聚合。如果你不希望看到计算噪声打开interpolationbilinear做平滑别用nearest。4.3 综合显示函数归一化、透明背景与矢量导出把上面几个动作收拢成一个函数方便直接换数据调用def plot_hht(t, imfs, f_max1000, vmin1e-3, cmapjet): energy build_spectrum(t, imfs) # 复用 2.2.3 的统计逻辑 mask freq_bins[:-1] f_max energy_show energy[mask, :] fig, ax plt.subplots(figsize(12, 5)) img ax.imshow(energy_show, aspectauto, originlower, extent[t[0], t[-1], freq_bins[:-1][mask][0], f_max], normLogNorm(vminvmin, vmaxenergy.max()), cmapcmap) ax.set_xlim(t[0], t[-1]) ax.set_ylabel(Frequency (Hz)) fig.colorbar(img, axax, labelEnergy) return fig, ax fig, ax plot_hht(t, imf_list, f_max1000, cmapturbo) fig.savefig(hht_optimized.pdf, bbox_inchestight)矢量图导出用PDF格式正文插图用PNG。PDF保留全部细节适合写报告和论文PNG按bbox_inchestight裁剪白边适合直接嵌入网页。此时你会注意到数据里的主要频率脊线已经清晰可见但颜色最亮的位置到底是什么频率这需要第5章的边际谱来验证而不是拿鼠标在图上猜。5. 用边际谱校验HHT时频图分解质量的自检方法5.1 归一化边际谱与瞬时频率自检时频图好看不等于分解正确。要验证HHT时频图反映的是真实物理频率标准做法是算边际谱把所有时刻的瞬时频率做直方图统计并按能量加权。边际谱的物理含义是“每个频率上信号累计消耗的能量”和傅里叶幅度谱类似但不要求线性平稳。def marginal_spectrum(imfs, t): freq_cum [] amp_cum [] for imf in imfs: h hilbert(imf) phase np.unwrap(np.angle(h)) f np.diff(phase) / (2 * np.pi * np.diff(t)) a np.abs(h)[:-1] freq_cum.append(f) amp_cum.append(a) freq_all np.concatenate(freq_cum) amp_all np.concatenate(amp_cum) return freq_all, amp_all freq_all, amp_all marginal_spectrum(imf_list, t) hist, bins np.histogram(freq_all, bins200, range(0, 150), weightsamp_all**2)对第2.2节的合成信号边际谱的峰值应该出现在20Hz到50Hz范围的两端附近120Hz处有一个小尖峰。因为扫频信号在每个频率上停留时间相同理想情况是20Hz和50Hz两端能量密度最高。如果边际谱峰值完全偏离真实频率说明EMD把信号拆成了没有物理意义的纯数学分量。5.2 用一组阈值判断分解质量检查项通过阈值处理动作边际谱主峰频率与理论中心频率误差5%通过保留当前参数120Hz分量峰与基波峰能量比干扰峰能量低于主峰不满足则调noise_width或trials瞬时频率中位数与理论扫频范围重叠中位数在25-45Hz区间不满足则回查IMF判据负频率占比1%大于1%则加大中值滤波窗口这个表格不是拍脑袋定的它对应的是“分解结果偏离物理事实”的最小可接受范围。真实数据里没有理论频率可对比此时检查第三项瞬时频率的中位数应该在信号的已知频带内如果连中位数都跑到不可信区间最稳妥的做法是回到第2步降低SD阈值重新分解并检查是否该换EEMD。测完边际谱后把marginal_spectrum的结果存成CSV和时频图放在一起。下次再改任何参数先对比边际谱峰值偏移量再决定是否保留改动。这套自检流程跑通后HHT时频图的问题就只剩下审美了。本文还有配套的精品资源点击获取