HHT时频图实战：EMD分解与MATLAB实现及调参技巧

简介：HHT时频图（希尔伯特-黄变换）是一种针对非线性、非平稳信号的时频分析方法，常用于机械故障诊断、生物医学信号和地震数据分析。压缩包提供一个基于MATLAB的HHT时频图实现脚本，全包共1个m文件，容量仅1KB，适合希望快速上手EMD分解与瞬时频率绘制的学习者参考。代码通过经验模态分解（EMD）将信号拆分为本征模态函数（IMF），再经希尔伯特变换获得瞬时幅度与频率，最终组合为时间-频率-幅度的时频分布图；对于初次接触HHT或想排查时频图绘制问题的读者，运行该脚本能直观理解完整算法流程。资源目前已有441人浏览学习，作为轻量级示例，覆盖了从数据预处理、EMD分解到时频图后处理的完整链路，也便于在此基础上修改以适配自身的非平稳信号数据。

1. HHT时频图为什么值得自己调

接手过轴承振动数据分析的人大都遇到过这种尴尬：FFT只能告诉你哪些频率成分存在，却说不清这些频率是什么时候出现的。面对变转速、冲击、裂纹扩展这类非线性非平稳信号，频谱图上的峰值往往是一片模糊的平均结果。HHT时频图之所以在机械故障诊断、生物医学信号分析、地震波处理里被反复提起，是因为它把“频率随时间变化”这件事直接画成了一张二维图——横轴时间、纵轴频率、颜色深度代表瞬时能量。与短时傅里叶变换需要提前选窗函数不同，HHT先通过经验模态分解把信号拆成本征模态函数，再对每个IMF做希尔伯特变换得到瞬时频率。这样做的好处是基函数来自信号本身，不预设固定窗宽，对瞬态冲击和频率调制更敏感。不过HHT的实现细节比教科书上写的要敏感得多，端点效应、筛分停止条件、IMF判据都会直接影响时频图的形态。如果你手上正好有那份包含Untitled.m的HHT时频图示例，接下来这套从原理到排错的拆解，能帮你把这个工具真正落到自己的数据上。

2. 特征尺度分离：EMD分解的机理与IMF边界

2.1 从包络局部均值到筛分迭代

经验模态分解最核心的思路是“借由信号的局部极值特征分离不同尺度的振荡”。给定一个离散信号 (x(t))，EMD先找出所有局部极大值和极小值，用三次样条插值分别构造上包络 (e_{\max}(t)) 和下包络 (e_{\min}(t))，取二者的均值作为局部均值包络 (m_1(t)=(e_{\max}(t)+e_{\min}(t))/2)。从原信号中减去这个包络得到第一个候选分量 (h_1(t)=x(t)-m_1(t))。问题是，一次相减往往不够，因为样条插值会产生新的极值点，所以必须对 (h_1(t)) 重复上述过程，直到满足IMF条件。这个反复相减的操作叫筛分（sifting），它本质上是把信号中叠加的高频振荡逐层剥离出来。

在MATLAB中，如果你用的是R2018a以上版本，可以直接调用内置的emd函数。以一段仿真的调幅调频信号为例：

fs = 2000; % 采样率 2000 Hz t = (0:1/fs:1)'; % 1秒时长 f1 = 50; % 基频 50 Hz x = sin(2*pi*f1*t).*(1 + 0.5*cos(2*pi*2*t)) + 0.3*sin(2*pi*350*t + 0.1*t.^2); [imf, residual] = emd(x, 'Display', 1);

emd默认返回IMF矩阵和残余项，每一列是一个IMF。Display参数会打印筛分迭代次数，方便你观察分解过程。这里的第一个IMF对应350 Hz附近的调频分量，第二个IMF对应50 Hz的调幅分量，残余项则是直流或低频趋势。

2.2 希尔伯特变换为什么能给出瞬时频率

每个IMF都是窄带信号，才能用希尔伯特变换定义有物理意义的瞬时频率。对IMF (c_i(t)) 做希尔伯特变换得到 (c_i(t) + j\mathcal{H}[c_i(t)])，它的解析信号幅度是 (a_i(t)=\sqrt{c_i(t)^2+\mathcal{H}[c_i(t)]^2})，相位是 (\theta_i(t)=\arctan(\mathcal{H}[c_i(t)]/c_i(t)))。瞬时频率定义为相位对时间的导数除以 (2\pi)，也就是 (f_i(t)=\frac{1}{2\pi}\frac{d\theta_i(t)}{dt})。

MATLAB的hilbert函数直接返回解析信号，不需要手动构造。计算瞬时频率时要特别注意相位差分。直接对unwrap后的相位用diff会产生一个比原信号少一个点的频率序列，画图前需要对齐时间轴：

analytic = hilbert(imf(:,1)); % 取第一个IMF inst_amp = abs(analytic); inst_phase = unwrap(angle(analytic)); inst_freq = diff(inst_phase)/(2*pi) * fs; % 结果长度为 N-1 inst_freq = [inst_freq; inst_freq(end)]; % 末值填充

代码里先对相位做unwrap消除周期性跳变，再除以采样间隔得到瞬时频率。末尾用最后一个值填充，是为了让inst_freq和inst_amp长度一致，方便后续画图或矩阵拼接。如果跳过这一步，plot会因为长度不匹配报错或画出错位的曲线。

2.3 IMF判据的两种误区

第一个误区是把“极值点个数与过零点个数相等或最多差一个”当成唯一判定条件。这个条件只是必要条件，实际上还需要局部均值趋于零，也就是上下包络对称。很多自实现EMD代码只看极值点数，导致分解出来的“IMF”根本不满足窄带要求，希尔伯特谱上出现负频率或交叉频率。第二个误区是认为筛分次数越多越好。筛分次数过多会把振幅调制抹平，让瞬时幅度失去物理意义。MATLAB内置emd默认相对容差是0.2，最大筛分次数100，你可以在代码中通过SiftRelativeTolerance和MaxNumSifting调整。调参时观察每个IMF的包络振幅是否还有明显波动，如果包络变成一条平滑直线，就说明筛分过头了。

调整这些参数时，建议先用内置默认值分解一次，画出所有IMF，观察哪些分量在物理意义上对应目标特征，再针对性地收紧或放松容差。我曾经处理一组带有强低频趋势的振动信号时，把SiftRelativeTolerance从0.2改到0.05后，第二个IMF从趋势项中分离出了原本被吞掉的10 Hz转频边带。代价是分解时间从不到1秒增加到3秒，但对于离线数据的离线分析，这个成本完全可接受。

3. 从EMD到HHT时频图：MATLAB完整实现

3.1 数据预处理：去趋势与端点处理

直接对原始加速度计信号做EMD，经常会被直流分量和低频漂移干扰。预处理的第一步是去掉均值，必要时用高通滤波器滤除0.1 Hz以下的趋势项。注意不要用陷波器滤除工频，因为陷波滤波器在频域会产生群延迟，导致时频图中的瞬时频率在工频附近扭曲。更稳妥的做法是使用detrend函数，它默认去除线性趋势，对非线性漂移则需要先用EMD分解出残余项再扣除。我通常的做法：

x = x - mean(x); % 去均值 x = detrend(x, 'constant'); % 再次确认均值归零 % 若有缓慢趋势，先做一次EMD取最后一个IMF作为趋势 [~, res] = emd(x); x = x - res;

这段代码先去除均值和常数趋势，然后利用EMD自身把残余项作为慢变趋势从原信号里扣除。注意这里第二次EMD只取残余项，不关心中间IMF，所以不需要保存完整分解结果。如果数据采样率很高，比如5 kHz以上，建议先做低通抗混叠滤波到下采样，否则EMD会把噪声分解成大量低能量IMF。

3.2 对每个IMF计算瞬时频率与幅值

分解完成后要对每个IMF依次调用hilbert。实际项目中IMF数量通常为5到10个，但并不是所有IMF都值得纳入时频图。低频残余项和能量占比极小的IMF既消耗内存，又会把时频图的颜色动态范围拉低。所以我一般在计算完解析信号后，统计每个IMF的RMS能量，只保留能量超过总能量1%的分量。

[imf, ~] = emd(x); num_imf = size(imf, 2); t_freq = cell(1, num_imf); t_amp = cell(1, num_imf); for k = 1:num_imf analytic = hilbert(imf(:,k)); phase = unwrap(angle(analytic)); freq = diff(phase) * fs / (2*pi); freq = [freq; freq(end)]; t_freq{k} = freq; t_amp{k} = abs(analytic); end

这里用diff计算瞬时频率，会放大小相位扰动。噪声大的IMF会出现频率尖刺，后续可以用中值滤波器对频率曲线做平滑，但平滑窗口不要超过信号最小周期的十分之一。比如信号最高分析频率为500 Hz，对应周期2 ms，窗口取0.2 ms，在采样率2000 Hz下就是不到1个点，实际很少做平滑，而是靠IMF本身窄带特性保证频率曲线光滑。

3.3 绘制时频图的三条实现路径

绘制HHT时频图最常见的方式是把时间-频率平面划分为网格，将每个IMF的瞬时幅度填充到对应的网格位置。MATLAB里imagesc搭配accumarray是最快实现。先定义频率轴和时间轴，再把瞬时频率四舍五入到频率网格索引，最后用accumarray累加幅度：

freq_axis = 0:5:1000; % 频率网格，分辨率5 Hz time_axis = t; % 时间轴与原信号一致 [~, freq_bins] = histc(t_freq{1}, freq_axis); hht_spectrum = zeros(length(freq_axis), length(time_axis)); for k = 1:num_imf valid = freq_bins > 0 & freq_bins < length(freq_axis); idx = sub2ind(size(hht_spectrum), freq_bins(valid), round(linspace(1,length(time_axis),sum(valid)))); hht_spectrum(idx) = hht_spectrum(idx) + t_amp{k}(valid).^2; end imagesc(time_axis, freq_axis, hht_spectrum); set(gca,'YDir','normal'); xlabel('时间 (s)'); ylabel('频率 (Hz)'); colorbar;

histc把瞬时频率映射到频率网格标号。sub2ind将二维索引转换为一维线性索引，便于快速累加。这里幅值取平方是因为时频图中我们希望显示能量，而不是幅值的线性值。set(gca,'YDir','normal')修正imagesc默认的倒置纵轴。如果数据量很大，建议预先分配hht_spectrum为sparse矩阵，并用sparse累加，避免密集矩阵占用过多内存。

对于需要矢量输出的论文，可以用pcolor但会特别慢。我常用的替代方案是把时频图绘制为散点图，每个点有瞬时时间和频率，颜色映射到瞬时幅度：

all_freq = []; all_time = []; all_amp = []; for k = 1:num_imf all_freq = [all_freq; t_freq{k}(1:10:end)]; all_time = [all_time; t(1:10:end)]; all_amp = [all_amp; t_amp{k}(1:10:end)]; end scatter(all_time, all_freq, 5, all_amp, 'filled');

注意这里的1:10:end是每隔10个点抽样一次，防止散点太多覆盖噪声区。这个方式的优点是能精确表达瞬时频率的波动，缺点是对频率变化剧烈的区域容易形成密集色块。

4. 时频图失真的三个根源与对应调参对策

4.1 端点效应：HHT时频图最常见的“假频率”

EMD构造包络时，第一个点和最后一个点往往没有完整的极值邻域，三次样条在端部会产生大幅度摆尾，导致IMF在开头和结尾出现异常振荡，进而让瞬时频率在端点附近突变。反映在时频图上，就是图像左右边出现竖直的亮条，频率值远超正常范围。

处理端点效应的常用做法是镜像延拓。在emd调用中，ExtrapolationMethod参数可以设为'mirror'，让函数在端点处对称复制极值。如果使用自实现EMD，可以在数据两端各延长一个周期：

x_ext = [flipud(x(1:200)); x; flipud(x(end-199:end))]; [imf_ext, ~] = emd(x_ext, 'ExtrapolationMethod', 'mirror'); % 截取原始数据对应部分 imf = imf_ext(201:end-200, :);

镜像延拓会人为引入周期假设，对非平稳信号来说不一定完全正确，但至少能消除包络幅值发散。另一种更轻量的做法是在绘制时频图时直接裁剪掉首尾10%的区域，只展示中间稳定段。这个方法比较粗暴，但适用于很长的信号，因为截断后信息损失有限。

4.2 筛分停止准则对频率分辨率的实际影响

希尔伯特变换要求IMF瞬时频率不能有负值，但实际上几乎所有实测数据的瞬时频率都会出现局部负频率，只是持续时间极短。负频率出现的原因通常是IMF局部不满足窄带条件，也就是相邻两个振荡的幅度差距过大。增加筛分次数能让IMF的包络更对称，降低负频率出现概率，但过度筛分会降低信号的时间分辨率。

在MATLAB内置emd中，SiftRelativeTolerance控制的是两次筛分之间候选IMF的能量变化率。设为0.5会提前停止，IMF可能不平滑；设为0.01会接近完全收敛，但耗时显著。对于有冲击特征的数据，我建议设为0.1到0.2，并且统计瞬时频率中负值占比：

negative_ratio = sum(inst_freq < 0) / length(inst_freq);

如果这个比例超过1%，就把SiftRelativeTolerance调小到0.05。如果调小后仍然有很多负频率，问题一般出在数据本身存在不连续点，比如传感器信号截断处。此时需要先对信号做平滑处理或分段分析，不要盲目调小容差。

4.3 频率轴上限与色标动态范围

时频图的纵轴频率上限由采样率决定，理论上最高是奈奎斯特频率。实际显示中如果直接画到奈奎斯特频率，大部分区域颜色会很暗，只有低频部分有亮色。图表的对比度和可读性都很差。我一般将频率轴上限设为关注频带最高频率的1.5倍，比如轴承故障特征频率集中在300到800 Hz，画图时只画到1200 Hz。

色标动态范围使用分位数压缩比线性映射更稳。直接线性映射会把少数高能量点拉高整体色标，导致时频图一片低亮度。用prctile找到第99百分位的幅度值，将上限设置成该值，剔除异常大点：

amp_max = prctile(nonzeros(hht_spectrum(:)), 99); imagesc(time_axis, freq_axis, hht_spectrum, [0 amp_max]);

imagesc的第四个参数直接指定色标上下限，让颜色映射集中在有效能量范围。这样处理后的时频图能清楚看到频带随时间漂移的轨迹，而不是只能看到几个极亮点。

5. 用边际谱验证HHT时频图质量的一个技巧

时频图画出来之后，很多人的下一步是直接读图找特征频率带。但图像容易受色标、噪声和端点效应干扰，一个更客观的验证方法是计算边际谱。边际谱的定义是HHT谱对时间积分：(h(f)=\int_0^T H(t,f),dt)，它表示整个时间长度上每个频率成分积累的总能量。与FFT幅度谱不同，边际谱不要求信号平稳，能更真实反映非平稳信号中“出现过的频率”的能量分布。

在MATLAB中，可以直接对hht_spectrum沿时间轴求和：

marginal_spectrum = sum(hht_spectrum, 2); plot(freq_axis, marginal_spectrum);

拿到边际谱之后，对比FFT谱中的峰值。如果某个频率在FFT谱中有明显峰，但边际谱中几乎没有能量，说明这个频率成分在时间上是极短的瞬态，或者被EMD分解到了残余项里。反过来，如果边际谱有峰而FFT谱没有，说明该频率成分持续时间长且幅度时变，这正是HHT独有的优势。

一个我常用来定位滚动轴承外圈故障的技巧是：先对原始信号做EMD分解，取前两个IMF，计算它们的瞬时频率和瞬时幅度，然后绘制二维平面下瞬时频率随时间的变化曲线，用颜色叠加瞬时能量。由于外圈故障会产生周期性的冲击，每次冲击对应的瞬时频率会有一个先升后降的轨迹，在时频图上表现为一条条短竖线。如果这些短竖线之间的时间间隔等于理论故障频率的倒数，就可以确认故障特征。为了提高信噪比，先对瞬时频率曲线做5点中值滤波，清楚随机尖刺。同时计算瞬时能量大于其均值的两倍的时刻，标记为冲击发生点：

energy = t_amp{1}.^2; threshold = 2 * mean(energy); impact_idx = find(energy > threshold); impact_times = t(impact_idx);

最终将impact_times与故障特征频率在时间上的间隔做对比，如果间隔的标准差小于采样周期的3倍，就能稳定判定故障类型。这个技巧的关键在于瞬时能量的阈值选取，阈值过高会漏掉弱冲击，过低会混入噪声。数据量大时，可以先用边际谱确定故障特征频率范围，再回来细化阈值参数。

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