企业官网建设流程全解析

简介：希尔伯特-黄变换（HHT）时频图是分析非线性、非平稳信号的重要工具。这份压缩包提供一段MATLAB实现代码，面向信号处理方向的学生、科研人员与工程开发者，可用于快速理解HHT时频图的完整生成流程，并迁移至生物医学信号、机械故障诊断、地震数据分析等应用场景。包体十分精简，仅含1个m文件，大小约1KB，代码短小但覆盖核心步骤：经验模态分解（EMD）得到IMF分量，希尔伯特变换计算瞬时频率与幅度，再组合绘制成时频分布图，便于对照理论逐一学习。已有441人学习下载。通过运行并修改这段代码，可直观观察不同参数和IMF分量对时频图的影响，掌握非线性非平稳信号分析的实现要点，避开HHT应用中的常见误区，为后续深入研究和工程实践打下基础。

1. HHT时频图的坑：从“毛刺满屏”说起

从网上下载一个写着“HHT时频图”的压缩包，解压跑完，得到的往往不是论文里那种干净利落的频率脊线，而是一张布满毛刺的彩色噪点图。这不是你下载错了代码，而是HHT时频图本身对分解参数和绘制方式极度敏感。HHT（Hilbert-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, endpoint=False) # 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], weights=aa[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, shading='auto', cmap='jet') ax.set_ylim(0, 150) ax.set_xlabel('Time (s)') ax.set_ylabel('Frequency (Hz)') fig.colorbar(mesh, ax=ax, label='Energy') plt.tight_layout() plt.savefig('hht_spectrum.png', dpi=300)

用shading='auto'可以消除pcolormesh对网格边界的抱怨。频率上限先放到150Hz，120Hz干扰那一条可以看见；下限不硬截，留给后续端点处理。到这里，一张最原始的HHT时频图已经能跑出来了，但大概率你会看到三个问题：左右两端出现几条斜冲上去的竖线，120Hz那条横线抖成波浪，低频区域有一片连续色块。这三个问题对应第3章的三类修正。

3. 画好HHT时频图的参数修正表：端点效应、模态混叠与负频率

3.1 端点发散：时频图两端的“飞翼”与镜像延拓

时频图最显眼的毛病是两端出现向上或向下的频率飞翼，频率值瞬间冲到几百赫兹，颜色异常亮。原因是样条包络在信号首尾处没有足够极值点约束，包络线会过冲，导致IMF在端点处产生形变。这类形变在做Hilbert变换后直接反映为瞬时频率的剧烈波动。

常见做法是给信号做镜像延拓：以首尾极值点所在时刻为镜面，把信号往外翻一段，让样条包络在端点处有数据可依。PyEMD的EMD类已经内置了镜像延拓选项，通过extrema_detection和NEW_MIRROR参数控制。实测中，只延拓两端各一个周期的数据，就能消除大部分飞翼。

emd_ext = EMD(extrema_detection='parabol') 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, trials=50, noise_width=0.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)))算瞬时频率，有一个数学上的副作用：相位差分会放大高频噪声，噪声点处的相位抖动会把瞬时频率压成负值或推向极高值。表现在时频图上，就是低频区域出现整片色块，或者频率轴上出现垂直的亮条纹。

处理这个问题的优先顺序是：

1. 先对Hilbert变换前的IMF做轻平滑，比如3点中值滤波，去除孤立极值点。
2. 再对瞬时频率做移动平均，窗口长度取信号周期的1/10左右。
3. 最后做频率掩膜，只保留全时间轴能量占比前95%的频率范围。

from scipy.ndimage import median_filter freq_smooth = median_filter(freq, size=3)

中值滤波对尖峰脉冲的抑制能力比均值滤波强，且不会模糊真实频率跳变点。如果你的信号本身存在频率突变（比如转速阶跃），不要在瞬时频率上做时间窗过长的平滑，否则阶跃会被拉成斜坡。

3.4 已踩过的参数整理成速查表

表格里每一项改动都不应该一次全做。我一般一次只动一个参数，画两次图对比，否则改了六个参数之后检出问题，根本不知道是哪一步修复了它。如果你时间紧，优先改飞翼和混叠两项，其余用绘图层面的掩盖手段处理更省事。

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, aspect='auto', origin='lower', extent=[t[0], t[-1], freq_show[0], freq_show[-1]], norm=LogNorm(vmin=1e-3, vmax=energy.max()), interpolation='bilinear', cmap='jet') fig.colorbar(img, ax=ax, label='Log Energy')

extent参数定义了图像的四边分别在数据坐标里的位置，aspect='auto'让x和y方向独立缩放，不会因为频率轴单位不同而压扁波形。LogNorm是给能量图配的对数色标，适合动态范围跨了几个数量级的场景。vmin不要设为0，对数坐标对0无定义，这里用1e-3做下限。

提示：imshow默认把数组当图像逐像素显示，不做任何聚合。如果你不希望看到计算噪声，打开interpolation='bilinear'做平滑，别用nearest。

4.3 综合显示函数：归一化、透明背景与矢量导出

把上面几个动作收拢成一个函数，方便直接换数据调用：

def plot_hht(t, imfs, f_max=1000, vmin=1e-3, cmap='jet'): 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, aspect='auto', origin='lower', extent=[t[0], t[-1], freq_bins[:-1][mask][0], f_max], norm=LogNorm(vmin=vmin, vmax=energy.max()), cmap=cmap) ax.set_xlim(t[0], t[-1]) ax.set_ylabel('Frequency (Hz)') fig.colorbar(img, ax=ax, label='Energy') return fig, ax fig, ax = plot_hht(t, imf_list, f_max=1000, cmap='turbo') fig.savefig('hht_optimized.pdf', bbox_inches='tight')

矢量图导出用PDF格式，正文插图用PNG。PDF保留全部细节，适合写报告和论文；PNG按bbox_inches='tight'裁剪白边，适合直接嵌入网页。此时你会注意到，数据里的主要频率脊线已经清晰可见，但颜色最亮的位置到底是什么频率？这需要第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, bins=200, range=(0, 150), weights=amp_all**2)

对第2.2节的合成信号，边际谱的峰值应该出现在20Hz到50Hz范围的两端附近，120Hz处有一个小尖峰。因为扫频信号在每个频率上停留时间相同，理想情况是20Hz和50Hz两端能量密度最高。如果边际谱峰值完全偏离真实频率，说明EMD把信号拆成了没有物理意义的纯数学分量。

5.2 用一组阈值判断分解质量

这个表格不是拍脑袋定的，它对应的是“分解结果偏离物理事实”的最小可接受范围。真实数据里没有理论频率可对比，此时检查第三项：瞬时频率的中位数应该在信号的已知频带内；如果连中位数都跑到不可信区间，最稳妥的做法是回到第2步，降低SD阈值重新分解，并检查是否该换EEMD。

测完边际谱后，把marginal_spectrum的结果存成CSV，和时频图放在一起。下次再改任何参数，先对比边际谱峰值偏移量再决定是否保留改动。这套自检流程跑通后，HHT时频图的问题就只剩下审美了。

本文还有配套的精品资源，点击获取