HHT变换MATLAB例程：EMD分解与希尔伯特时频分析实战详解

简介：这是一份面向MATLAB信号处理学习者的希尔伯特黄变换（HHT）例程包，围绕希尔伯特变换与经验模态分解（EMD）详细展示时频分析实现思路。压缩包共3个文件，全部为.m脚本，体积仅10KB。三个脚本分别承担主程序、核心EMD分解与特征提取功能，代码模块划分清晰，便于按功能拆分学习和调试。已有199人浏览学习，适合需要快速上手HHT原理与代码的初学者。通过这套例程，用户可以了解利用hilbert()函数构建瞬时频率、编写EMD迭代分解流程，并基于HHT谱进行信号特征提取；代码体量小、注释直观，适合直接运行修改，有助于理解非平稳信号在生物医学、机械故障诊断等场景中的时间频率分析基本方法。

1. 拿到振动数据先别急着画频谱，这组 HHT 例程其实更值得跑一遍

一段轴承振动信号放到 MATLAB 里，多数人第一步就是pwelch求功率谱。可当信号频率随时间漂移、又叠着调幅冲击时，固定窗长的频谱会把不同时刻的频点搅成一个平均意义上的峰，包络里真正有用的突变反而被平滑掉了。HHT.zip这组 matlab 例程走的是另一条路：emd.m负责把信号自适应拆成若干固有模态函数，hht1204.m把分解和希尔伯特变换串成完整时频分析流程，hht_feature_extraction.m再把时频信息压缩成能喂给分类器的特征向量。整套流程不依赖第三方工具箱，对做机械故障诊断、生物医学信号分析或波浪数据处理的人来说，属于值得完整过一遍源码的 matlab 例程。下面按真正会去用它的顺序拆解。

2. EMD 分解原理与 emd.m 调用约定

2.1 包络均值减掉的是"信号里的本征波动"

经验模态分解的核心思想很朴素：任意复杂信号可以看成若干个本征波动叠加。所谓本征模态函数 IMF，需要满足两个条件：极值点个数与过零点个数相差不超过 1；上下包络线在任意时刻的均值是 0。实际的做法不是直接求解方程，而是靠迭代筛分。

筛分过程通常是这样：先用findpeaks找出局部极大值和极小值，分别插值生成上包络、下包络，取两者的均值 m，再用原始信号减去这个均值得到细节分量 h。如果 h 不满足 IMF 条件，就把 h 当成新的输入继续筛，直到满足停止准则。这中间有两个常见坑：插值方式选'spline'还是'pchip'影响包络形态；筛分次数太多会把 IMF 洗成纯调频信号，振幅波动被吸干，失去物理意义。

包里的emd.m一般会暴露几个可调参数，调用格式通常是[imf, residual, info] = emd(x)。下面这段代码可以用合成信号先验证分解行为：

clear; close all; addpath('./hht'); % HHT.zip 解压后的目录名按实际改 fs = 1000; % 采样率，单位 Hz t = (0:999) / fs; x = sin(2*pi*50*t) + 0.5*sin(2*pi*125*t) + 0.1*randn(size(t)); [imf, residual] = emd(x); % 使用本例程包里的 emd.m figure; for k = 1:min(3, size(imf,1)) subplot(4,1,k); plot(t, imf(k,:)); ylabel(['IMF', num2str(k)]); end subplot(4,1,4); plot(t, residual, 'k'); ylabel('残差');

imf是模态矩阵，每一行对应一个 IMF，按频率从高到低排列，列数等于信号长度；residual是最后剩下的趋势项。第一次跑通时先不要上真实数据，用正弦叠加白噪声确认每个 IMF 的频率成分是否清晰。若某个 IMF 同时包含 50 Hz 和 125 Hz，说明筛分没有收敛或者信号长度不够。

2.2 这类 emd.m 内部常用的停止准则参数

这类例程里的emd.m表面只接收信号，内部却留着几个可调挡位。常见参数如下：

一般建议先用默认参数跑，再重点观察后几个 IMF 是否出现明显的频率混叠。如果末尾几个 IMF 形似噪声而不是光滑波形，多半是MAXMODES没设上限，或者停止阈值太小。

3. 用 hilbert() 把 IMF 变成瞬时幅度与瞬时频率

3.1 从窄带 IMF 到解析信号

得到 IMF 之后，希尔伯特变换只是配角却至关重要。对每个 IMF 做hilbert(imf)，相当于构造复信号 z(t) = x(t) + j·H[x(t)]，这样就能同时得到瞬时幅度、瞬时相位和瞬时频率。瞬时幅度是解析信号的模，瞬时相位是解析信号的角度，瞬时频率是相位对角度的导数。关键是 HHT 里的频谱不再像 FFT 那样把信号投影到固定基函数上，而是逐时刻给出频率，所以能够描述频率随时间连续变化的信号。

使用 EMT 分解后的 IMF 再做希尔伯特变换，前提是 IMF 足够窄带。如果直接把原始信号丢给hilbert，瞬时频率会在多个频率分量之间剧烈跳变，得到的瞬时频率没有任何物理意义。这也是为什么必须先跑emd.m再跑hilbert()，顺序不能反过来。

3.2 组装时频平面：希尔伯特谱的计算代码

把每个 IMF 的瞬时频率和瞬时幅度按时间对齐，累加到时间-频率平面上，就得到希尔伯特谱。下面这个函数可以直接存成脚本使用：

function [hs, faxis, taxis] = hilbert_spectrum_from_imfs(imf, fs) % 输入 imf: emd() 输出的模态矩阵，行是模态，列是采样点 % 输出 hs: 频率×时间 的希尔伯特谱矩阵 n = size(imf, 1); ncol = size(imf, 2); taxis = (0:ncol-1) / fs; fmax = fs / 2; nbin = 256; % 频带划分数，可按需提高 fres = fmax / nbin; faxis = linspace(0, fmax, nbin); hs = zeros(nbin, ncol); for k = 1:n z = hilbert(imf(k,:)); % 解析信号 a = abs(z); % 瞬时幅度 phase = unwrap(angle(z)); % 相位展开，避免跳变 f = [diff(phase), 0] / (2*pi) * fs; % 瞬时频率，单位 Hz for i = 1:ncol if f(i) < 0 || f(i) > fmax continue; % 丢弃异常频点 end idx = floor(f(i) / fres) + 1; idx = min(idx, nbin); hs(idx, i) = hs(idx, i) + a(i); % 幅度能量累加到该频带 end end end

这个循环里最关键的是[diff(phase), 0]：相位差分后长度比原信号少 1，补 0 是为了让频率序列和幅度序列保持同长度，否则累加时索引会错位。负频率直接丢弃是因为负的瞬时频率通常由端点效应或相位跳变造成，没有物理含义。

绘制时频图用imagesc就够了：

figure; imagesc(taxis, faxis, hs); set(gca, 'YDir', 'normal'); xlabel('时间/s'); ylabel('频率/Hz');

set(gca, 'YDir', 'normal')一定要加上，否则频谱会被翻转成上大下小的图像，容易误读。

3.3 频率分辨率由数据本身决定，不是由窗函数决定

希尔伯特谱的频率分辨率取决于相位差分后的频带累计，和傅里叶频谱不同的是它没有固定的窗长概念。nbin=256只是一个显示分档，实际瞬时频率是连续的，把nbin提高并不会带来虚高的分辨率，只是让谱图更细腻。这是 HHT 和短时傅里叶最大的区别：STFT 时间分辨率和频率分辨率相互拖累，HHT 的时间分辨率等于采样周期，频率精度取决于相位差分质量。

4. hht_feature_extraction.m 特征提取：把时频图变成可量化指标

4.1 从希尔伯特谱到特征向量的常用维度

时频图直观，但故障诊断和模式识别需要的是定长特征向量。hht_feature_extraction.m这类文件在这个包里承担的就是最后一公里：把emd.m和希尔伯特谱的结果压缩成一组数字。常见做法是围绕边际谱和瞬时能量做统计。边际谱就是把希尔伯特谱沿时间轴求和：margin = sum(hs, 2)，它表示每个频率上累计的能量密度。相比 FFT 频谱，边际谱的峰更集中，对调幅冲击更敏感。

常用的特征项有这些：

4.2 按工程惯例补全一个可用的特征提取函数

打包文件里的hht_feature_extraction.m通常只负责把流程串起来，参数细节需要自己补。我一般用它写成如下结构：

function feat = hht_feature_extraction(x, fs) % x: 单段信号，列向量 % 输出 feat: 8 维特征向量 imf = emd(x); % 第 1 步：分解 % 第 2 步：组装希尔伯特谱 [hs, faxis] = hilbert_spectrum_from_imfs(imf, fs); % 第 3 步：边际谱统计特征 margin = sum(hs, 2); % 沿时间求和 margin = margin / (sum(margin) + eps); % 归一化 [~, idx_peak] = max(margin); feat(1) = faxis(idx_peak); % 峰值频率 feat(2) = margin(idx_peak); % 峰值能量 % 低频段能量比：0~200Hz 之间的能量占比 low_band = faxis <= 200; feat(3) = sum(margin(low_band)); % 瞬时能量序列的均值和方差 inst_energy = sum(hs, 1); feat(4) = mean(inst_energy); feat(5) = std(inst_energy); % 边际谱熵 margin_safe = margin(margin > 0); feat(6) = -sum(margin_safe .* log(margin_safe)); % 前 3 个 IMF 的瞬时幅度峭度 z1 = abs(hilbert(imf(1,:))); z2 = abs(hilbert(imf(2,:))); z3 = abs(hilbert(imf(3,:))); feat(7) = kurtosis(z1); feat(8) = kurtosis(z2) + kurtosis(z3); feat = feat(:); end

这里要注意两点：边际谱归一化放在熵计算之前，保证熵值不受绝对能量影响；峭度作为无量纲指标，对冲击型故障非常敏感，但只用第一个 IMF 容易被噪声主导，所以把前三个 IMF 的峭度合并。如果你的信号采样率不是 1000 Hz，低频段阈值 200 需要按实际特征频率调整。

4.3 特征稳定性与归一化问题

EMD 对样本扰动比较敏感，轻微噪声可能导致分解出的 IMF 顺序或形态变化。为了特征稳定，常见做法是每个样本重复提取若干次，加不同随机种子的小噪声扰动后取平均，把最后的结果作为该样本的特征向量。送入分类器前建议再对每个特征维度做 z-score 归一化，否则边际谱熵和峰值频率量纲差异太大，树模型之外的大部分分类器都会偏向量纲大的特征。

5. 滚动轴承故障诊断场景下把例程整体跑通

5.1 分段与重叠率设置

真实工程里的故障信号不是无限延长的，特征提取前通常要加窗分段。滚动轴承数据常见的采样率是 12 kHz 到 48 kHz，每段选择 1024 或 2048 点比较常见。分段过短，低频特征频率在谱里难以分辨；分段过长，故障冲击会被平均掉。

读取和分段的代码片段如下：

load bearing_vibration.mat; % 假设变量名为 x，长度 N fs = 12000; win = 2048; % 每段样本数 step = 1024; % 重叠 50%，保留更多样本 nseg = fix((length(x) - win) / step) + 1; seg = zeros(win, nseg); start = 1; for k = 1:nseg seg(:,k) = x(start:start+win-1); start = start + step; end

重叠率 50% 是常用做法：既保证相邻段之间有连续性，又不至于让重复样本过多而导致分类器过拟合。若样本量不足，可以把重叠率提高到 75%，代价是计算时间上升。

5.2 批量特征矩阵构建与标签对齐

分段完成后，对每段调用特征提取函数：

X = zeros(8, nseg); % 8 维特征 for k = 1:nseg X(:,k) = hht_feature_extraction(seg(:,k), fs); end % 假设 seires 是逻辑向量，1 表示故障 labels = repmat(0, nseg, 1); labels(idx_fault) = 1; % 特征矩阵转置为 样本×特征，供分类器使用 X = X'; save('hht_feat.mat', 'X', 'labels', '-mat');

这里常见的问题是特征提取函数内部会重新调用emd，计算量较大。批量跑上万段时建议用parfor替换for，同时每段长度固定，避免数组动态增长。对后续分类任务，用随机森林或 XGBoost 就够了，HHT 特征已经做了非线性压缩，不需要再用深度网络从头学。

5.3 实际调参建议表

故障诊断场景中，优先看边际谱峰值频率和第一个 IMF 的峭度。轴承内圈故障时峰值频率会落在特征频率及其倍频附近，峭度随冲击加剧而增大；如果峰值频率落在固有频率附近而不是故障特征频率上，需要考虑是否发生了共振。

6. 三个容易让 HHT 例程输出"漂亮但错误"的坑

6.1 端点效应比想象中严重

EMD 插值包络时，信号首尾只有一侧有极值点，包络会在这里出现明显摆动，导致首尾 IMF 变形、瞬时频率向极端值漂移。常见做法是镜像延拓法扩展端点，再进行分解：

% 数据延伸一倍后再分解，结束后截掉扩展部分 x_ext = [flipud(x(1:200)); x; flipud(x(end-199:end))]; imf_ext = emd(x_ext); imf = imf_ext(:, winning);

延拓长度取几百个点即可，太长会改变整个分解结果。如果没有延拓，至少要把结果的首尾各 5%~10% 丢弃，不纳入特征统计。

6.2 瞬时频率出现负值时不要强行归零

hilbert对近似线性调频信号做相位展开通常没问题，但遇到剧烈冲击时会算出负瞬时频率。直接把这些点设 0 在谱图上会产生一条虚假的低频线。更合理的处理是记录这些点的比例，当负频点占比超过 10% 时，说明该段信号不适合做 HHT，应该回到 EMD 参数去检查筛分质量。负频本身是很好的质量指示器，不要只想着掩盖它。

6.3 用线性调频信号验证整套例程的正确性

跑真实数据前，先用一个已知瞬时频率的合成信号验证emd.m和谱图代码是否正常：

fs = 2000; t = 0:1/fs:1-1/fs; f0 = 50; f1 = 300; % 频率从 50Hz 线性扫到 300Hz 的 chirp 信号 x = sin(2*pi*(f0*t + (f1-f0)*t.^2/2)); imf = emd(x); [hs, faxis, taxis] = hilbert_spectrum_from_imfs(imf, fs); % 提取谱峰作为实测瞬时频率 [~, idx] = max(hs, [], 1); f_est = faxis(idx); f_theory = f0 + (f1-f0)*t; % 理论瞬时频率 err = mean(abs(f_est(200:end-200) - f_theory(200:end-200))); fprintf('平均瞬时频率误差: %.3f Hz\n', err);

把首尾各 200 点去掉是为了避开端点效应。误差小于 1 Hz 说明分解、希尔伯特变换和谱图组装这段链路基本正确；误差偏大时，优先查unwrap是否漏掉相位跳变，再查 EMD 停止准则是否导致分解不彻底。验证通过后，再把这套流程交给真实数据做故障诊断或生物医学信号分析，输出才有信心。

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