简介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的包络振幅是否还有明显波动如果包络变成一条平滑直线就说明筛分过头了。参数内置默认值典型调整范围调整目的SiftRelativeTolerance0.20.05~0.5控制IMF收敛精度越小越严格MaxNumSifting10050~500限制单次筛分迭代次数防过筛MaxNumIMF103~15限制分解出的IMF数量防止过度分解Display00 或 1输出分解过程便于调试调整这些参数时建议先用内置默认值分解一次画出所有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倍就能稳定判定故障类型。这个技巧的关键在于瞬时能量的阈值选取阈值过高会漏掉弱冲击过低会混入噪声。数据量大时可以先用边际谱确定故障特征频率范围再回来细化阈值参数。本文还有配套的精品资源点击获取