简介：中央大学数据研究中心的EMD/HHT Matlab工具箱由黄老师推荐，面向信号时频分析、机械故障诊断、地球物理及生物医学等领域的Matlab研究者与工程师。核心部分基于经典EMD算法，并附EEMD扩展实现，便于对比研究模态混叠抑制效果。压缩包共45个文件，以41个.m源码为主，覆盖IMF分解、希尔伯特变换、包络/极值处理、数据规范化等完整流程；2个csv文件提供测试数据，另有1个.p文件保存预计算函数，1个txt说明辅助配置，整体仅188KB，轻量无需额外安装。目前已有2199人学习下载，属于高校实验室常备参考资料。仔细研读代码可掌握从局部极值检测到包络拟合、从模态分解到瞬时频率计算的每一步细节，还可借鉴blocknormalize、splinenormalize等多种归一化策略与边界处理方法，方便二次开发或与G-Rilling版本相互校验。
1. 中央大学 EMD HHT Matlab 源代码：为什么这个工具箱至今没人能绕过
做旋转机械故障诊断的工程师大概都有过这种经历：测到一段振动信号，FFT 一看峰值杂乱，特征频率被噪声和边带糊成一片，小波换了三四个基函数也分不清哪个是冲击哪个是转频。我当年遇到这类信号，第一个想到的不是自己写算法，而是去找一套能直接跑出结果的现成实现——中央大学的 EMD/HHT Matlab 源代码就是这么被捞出来的。它做的事情一句话能说清：把一段非平稳、非线性的信号自适应地拆成若干个本征模态函数（IMF），再对每个 IMF 做 Hilbert 变换，得到瞬时频率和时频谱。这套工具箱解决的痛点很具体：信号频率随时间变化时，FFT 只给你一张静态频谱，HHT 却能告诉你“哪个时刻出现了哪个频率”。适合做故障诊断、地震信号分析、生物医学信号处理和语音研究的从业者，也是新手入门 HHT 时最容易跑通的一条路。
2. EMD 和 HHT 到底是什么：先弄懂分解逻辑，再谈参数调优
2.1 EMD 的筛分过程：信号是怎么被一层层剥开的
经验模态分解的核心假设是：任何复杂信号都可以看成若干个 IMF 叠加的结果，而每个 IMF 必须满足两个条件——极值点个数与过零点个数相等或最多相差一个；上下包络关于时间轴局部对称。满足这两个条件的分量才有资格做 Hilbert 变换，因为只有窄带信号算出来的瞬时频率才有物理意义。
具体实现时，算法对原始信号反复执行“筛分”。每一轮迭代做四件事：找出当前信号的局部极大值点和极小值点；用三次样条曲线分别拟合出上包络和下包络；计算上下包络的均值；用原信号减去这个均值。减完之后检查剩余信号是否满足 IMF 条件，不满足就继续重复这四步，满足就得到一个 IMF，然后用原信号减去这个 IMF，对剩余部分再开始新一轮筛分。这个过程很像剥洋葱，一层层去掉波动最快的分量，留下的残余项通常是趋势项或者常数偏移。
这套工具箱里的核心函数emd就是把上面的筛分过程封装成一句话调用。源码里用了三次样条插值做包络拟合，停止条件用的是两个相邻筛分结果的标准差阈值，默认值一般取 0.05 到 0.3 之间，具体看你要留多少细节。需要提醒的是，EMD 筛分过程没有任何解析公式，每一步都依赖数值计算，所以同一个信号用不同实现跑出来结果会有细微差别——这是正常的，不要觉得是自己代码写错了。
2.2 为什么不能直接用 FFT 或小波：HHT 的两步走优势
FFT 把信号映射到频域，等于假设信号是平稳的、频率不随时间变化。真实世界的振动信号、心电信号、地震波，哪个不是频率在漂的？小波变换虽然加了时间定位，但前提是你得先选一个小波基函数，选错了分辨率就不对，而且小波基和信号本身没有必然联系。HHT 的思路完全不同：IMF 是从信号自身的数据特征里“长”出来的，不需要预设基函数，这就是它常被叫做自适应时频分析方法的由来。
HHT 完整流程分两步：第一步用 EMD 把信号拆成 IMF；第二步对每个 IMF 做 Hilbert 变换，构造解析信号，算出瞬时幅值和瞬时频率。把瞬时频率画在时间轴上就得到 Hilbert 谱，再对时间积分就得到边际谱。边际谱的物理含义和 FFT 频谱不一样——FFT 的幅值代表“这个频率成分在整个时长里的总强度”，边际谱代表“这个频率出现的累计概率密度”，因此边际谱对瞬态冲击的检测往往更敏感。我见过不少故障诊断的案例，滚动轴承的内圈故障特征频率在 FFT 谱里被调制边带淹没，但在边际谱里却是一条清晰的峰，这就是 HHT 的实用价值。
2.3 这套工具箱在同类实现里的定位
Matlab 里不是没有别的 EMD 实现，但中央大学这套用得最广，主要原因是它对标准 EMD 算法的还原度足够高，而且提供了完整的后处理链条：分解用emd，频谱分析用hht，画三维时频谱用hhspectrum配合toimage。单看函数名你可能觉得没什么，实际动手就会发现，一个能直接从 IMF 算到边际谱的工具箱能省掉你两三天写胶水代码的时间。
另一类常见选择是那些把 EMD、EEMD、CEEMDAN 整合在一起的第三方扩展包，功能更全，但代价是依赖关系复杂、版本迭代频繁，跑通环境有时候比调算法还费劲。中央大学这套相对克制，核心就一条线：EMD 分解 + Hilbert 谱分析。对于只想把 HHT 用到自己数据上、不想折腾额外依赖的工程师来说，这个取舍恰恰是优点。下面三个章节的步骤和参数说明都以这套工具箱为基准。
3. 安装与首次跑通：从下载目录到画出第一张 IMF 图只需要十几分钟
3.1 把工具箱路径加入 Matlab 搜索路径
这套工具箱是以源代码文件形式散发的，没有安装包，也不存在install.m脚本。你下载解压后拿到的是一个文件夹，里面有核心函数的.m文件。第一步要做的是把这个文件夹连同子目录一起加入 Matlab 搜索路径，否则调用emd会直接报“未定义函数或变量”。
常见做法是写一行命令：
addpath(genpath('D:\emdsrc')); % 改成你解压后实际的路径这里genpath的作用是把指定目录下所有子目录一并加进来，比单独addpath一个根目录稳妥。执行完可以用which emd验证：
which emd如果返回了完整路径，说明搜索路径没问题，工具箱已经可以被正常调用。需要提醒的是addpath只在当前 Matlab 会话里有效，关掉软件重新打开就失效了。要么每次启动后手动执行，要么把这行命令写进startup.m，要么在“预设 -> 常规 -> 路径”里永久添加。我个人推荐用startup.m，因为写死在 Matlab 配置里一劳永逸，换电脑重装也不会忘。
3.2 构造一个测试信号跑通 EMD 分解
上手阶段别直接用真实数据，先用一个频率随时间变化的合成信号验证工具箱是好的。最经典的测试信号是线性调频信号叠加正弦波：
fs = 1000; % 采样率 1000 Hz t = 0:1/fs:1; % 时长 1 秒 x = sin(2*pi*(100 + 50*t).*t) + 0.5*sin(2*pi*30*t); % 第一项是频率从 100 Hz 线性扫到 150 Hz 的调频分量 % 第二项是 30 Hz 的固定频率分量 imf = emd(x); % 执行 EMD 分解这段代码里的x必须是行向量，这是这套工具箱的硬性要求。很多人在这一步翻车——数据读进来默认是列向量，直接传给emd会报维度错误。如果拿到的是列向量，先做一次转置：x = x(:)';。emd返回的imf是一个矩阵，每一行对应一个 IMF 分量，最后一行通常是残余项。分解完成后用size(imf)查看行数，行数多说明信号复杂度高，这是正常现象，不用慌。
接下来把分解结果画出来，直观验证每个 IMF 是否符合预期：
figure; for k = 1:size(imf, 1) subplot(size(imf, 1), 1, k); plot(t, imf(k, :)); title(['IMF ', num2str(k)]); end画完你应当看到：第一个 IMF 对应扫频信号的高频成分，第二个 IMF 对应 30 Hz 的正弦分量，最后一个残余项接近一条平线。这个结果说明工具箱工作正常，分解逻辑符合 EMD 的基本预期。如果你画出来的图出现第一行和第二行明显互相混叠、形态异常，先检查信号本身是不是太短——EMD 对短序列的支持很差，低于 20 个完整振荡周期的信号分解质量基本不可用。
3.3 确认分解质量的两个快速检查
一个信号分解完，光看波形还不够，你得确认拆出来的 IMF 是不是真的满足窄带条件。最简单的检查是算每个 IMF 的过零点数和极值点数，看它们是否相等或至多差一个：
for k = 1:size(imf, 1) zc = sum(abs(diff(sign(imf(k, :)))) > 0); % 过零点数 [~, pk] = findpeaks(imf(k, :)); % 局部极大值点 [~, tr] = findpeaks(-imf(k, :)); % 局部极小值点 fprintf('IMF%d: 过零点=%d, 极值点=%d\n', k, zc, length(pk)+length(tr)); end注意这里findpeaks是信号处理工具箱的函数，如果你的 Matlab 没有装该工具箱就换成islocalmax和islocalmin，效果等价。检查结果里前几个 IMF 会满足条件，越往后的 IMF 偏差越大——这不算算法错误，因为残留分量常常只包含单调趋势，极值点定义本身就失去意义。这个检查的意义在于帮你判断分解是否过度：如果第一行 IMF 的过零点数就明显违反条件，说明停止准则设得太宽松，需要调emd的阈值参数。
4. 核心函数逐个拆解：emd、hht、hhspectrum、toimage 到底怎么配合
4.1 主函数 emd 的输出结构和可选参数
emd是整个工具箱的入口，最基础的调用形式是imf = emd(x)。完整的调用形式是：
[imf, residual, info] = emd(x, 'SiftMaxIteration', 100, 'MaxIMF', 5);返回值有三个：imf是 IMF 矩阵，每行一个分量；residual是最终的残余趋势项；info是一个结构体，里面记录了每层分解的筛分迭代次数、停止准则的值等信息。info是一个很实用的调试工具——当某层分解迭代次数长时间不收敛，就要怀疑停止准则对当前数据过于严苛。
可选参数里，SiftMaxIteration控制单层筛分的最大迭代次数，默认值通常是 100。设得越大，每个 IMF 越光滑，但计算时间线性增长；设得太小，IMF 可能还混着相邻尺度的成分。MaxIMF限制最多分解出几层 IMF，适合你明确知道信号里有几个主要分量时使用。另有一个隐藏参数MAXITERATIONS控制整体分解外循环次数，一般不用动。还有SiftThreshold系列参数控制停止准则的具体阈值，后面第五章会展开讲。
参数调整有一个总原则：先跑默认参数，看结果再改。不要一上来就追求参数最优，因为 EMD 没有全局最优解一说，不同的参数组合对应不同的分解精细度，关键是匹配你后续分析的需求。做故障诊断，IMF 分量稍微混一点也许没关系；做瞬时频率精细估计，对停止准则就要严格一些。
4.2 从 IMF 到希尔伯特谱：hht 和 hhspectrum 的分工
拿到 IMF 之后，下一步是对每个分量做 Hilbert 变换。这套工具箱里，直接用hht函数就能算瞬时频率和瞬时幅值：
[inst_freq, inst_amp] = hht(imf, fs);inst_freq和inst_amp都是矩阵，行对应 IMF，列对应时间点。inst_freq的单位是 Hz，这和你传入的采样率fs保持一致；inst_amp是解析信号的包络幅值。注意这里的频率是瞬时频率，会随每个时刻的相位变化率波动，不要期望它像 FFT 频谱那样平滑。
如果你需要进一步画出 Hilbert 谱的三维图，就要用hhspectrum加toimage的组合：
[hs, f, t] = hhspectrum(imf, 1, fs); % 参数 1 表示对 imf 逐行做 Hilbert 变换 [t_img, f_img] = toimage(hs, f, t, length(t));这里hhspectrum返回的hs是复矩阵，实部是瞬时幅值，虚部是瞬时频率经过某种映射后的结果；f和t是对应的频率轴和时间轴向量。toimage再把稀疏的瞬时频率-时间点插值到规则网格上，生成一个矩阵t_img，这个矩阵就是可以直接imagesc或pcolor画出来的时频谱图像——横轴是时间，纵轴是频率，颜色深浅表示该时刻该频率上的能量强度。f_img是对应的频率坐标。
4.3 边际谱与瞬时频率的提取示例
画时频谱直观，但很多工程场景下你需要的是一条“谱线”用于特征提取，这时候算边际谱更实用。边际谱就是把 Hilbert 谱对时间积分：
t_img = abs(t_img); % 取幅值 marginal_spectrum = sum(t_img, 2); % 对时间轴求和 f_axis = f_img; % 频率坐标 plot(f_axis, marginal_spectrum);跑完这段代码你就得到了一条和 FFT 频谱形态类似的曲线，但每个峰的物理意义不同。拿轴承故障信号举例，FFT 谱里外圈故障特征频率的峰值可能被转频谐波盖住，边际谱里同一个频率的峰往往因为“该频率出现的累计能量”而显现得更清楚。这不是玄学，而是因为 HHT 把每一个时刻的瞬时频率都统计进去了，瞬态冲击成分不会被长时平稳成分稀释掉。
需要特别说明：marginal_spectrum的纵轴量纲取决于t_img取绝对值的方式。如果你用abs取模，得到的是能量幅值积分；如果你对复数谱取实部，得到的是有正负的瞬时幅值积分，习惯上取模的做法更常见。我在实际项目里通常两种都算一遍，对比看哪个对特征频率更敏感。
5. 避坑指南：这套工具箱藏得最深的 5 个坑，每一个我都踩过
5.1 输入数据必须是 double 型行向量，否则报错莫名其妙
现象：数据从 CSV 读进来，直接调用emd(x)，Matlab 提示“维度不匹配”或者干脆输出全 NaN。原因：工具箱内部用size(x, 2)推断信号长度，如果输入是列向量，size(x)返回的第一个维度是数据点数，函数就把整个逻辑搞乱了。解决：进入emd之前强制转换类型和形状。
if ~isa(x, 'double') x = double(x); end if size(x, 1) > 1 x = x(:)'; % 强制转成行向量 end这段转换代码应该形成肌肉记忆。我见过有人在这个坑上折腾了半天，以为自己工具箱没装对，实际上只是 Excel 导出时把行数据存成了列。
5.2 端点效应：信号两端出现巨大幅值飞翼
现象：分解结果的第一个 IMF 两端出现明显的发散抖动，振幅远大于信号本身；Hilbert 谱的边缘出现竖直亮条。原因：EMD 做三次样条包络拟合时，端点处找不出极值点，样条外推不可靠，导致包络失真。解决：最省事的方法是分解前对信号做端点延拓——镜像延拓是常见做法之一，把边界值对称折叠出去再分解，分解完裁掉两端延拓部分。
ext_len = round(length(x) * 0.1); % 延拓 10% 长度 x_ext = [x(ext_len:-1:1), x, x(end:-1:end-ext_len+1)]; % 左右各做一次镜像反射 imf_ext = emd(x_ext); imf = imf_ext(:, ext_len+1:end-ext_len); % 裁掉延拓段镜像延拓并不能完全消除端点效应，但能把飞翼推向延拓段，使有效区间内的分量更可信。分析时还有一个习惯：只看信号的中间 80% 区间，两端各舍弃 10%，这是很多论文里实际上在用的做法，只是通常不会写出来。
5.3 模态混叠：间歇性高频信号把分解搞乱
现象：信号里有一串周期性的冲击（比如轴承内圈故障的冲击序列），分解后这些冲击的能量没有集中在一个 IMF 里，而是散落在相邻两三个分量中，导致每个 IMF 都长得不伦不类。原因：EMD 的筛分过程按时间尺度逐层剥除成分，如果某一成分的出现是间歇性的，它和邻近频率成分的交织会让包络均值失真。解决：常见做法是用掩膜信号（mask signal）把间歇成分“勾”出来——选一个高于目标成分频率的正弦波叠加到原信号上，一起做 EMD，分解完再从对应 IMF 里减掉这个正弦波的贡献。还有一类做法是改走集合经验模态分解（EEMD）路线，给原信号多次加白噪声再平均，缺点是计算量倍增。聊到 EEMD 要提醒一句：这个方案不在中央大学这套工具箱的原生支持范围内，你得上手改源码或换第三方包。
5.4 过度分解：IMF 数量多到离谱，低频段毫无物理意义
现象：一段简单正弦信号，分解输出 8 个 IMF，后面几个频率很低、振幅也很小的分量怎么看都是噪声。原因：停止准则设得过严，算法不认为残余项已经“足够单调”，继续往下剥离出数值噪声级别的小波动。解决：调大停止阈值，或者限制最大分解层数。直接看参数怎么动：
[imf, residual, info] = emd(x, 'MaxIMF', 4, 'SiftThreshold', 1e-3);MaxIMF设为一个有限值是最直接的控制手段，你心理预期这段信号有几个物理成分，就设多少，额外多出来的层大概率都是分解过度产物。SiftThreshold是停止准则的关键参数，默认值在 1e-4 到 1e-5 量级，改到 1e-3 左右可以显著减少无效 IMF，代价是前几层 IMF 的分辨率有所下降。
5.5 瞬时频率出现负值：Hilbert 变换的物理性失灵
现象：计算出的瞬时频率在部分时刻为负，画在时频谱上出现向下的毛刺。原因：瞬时频率的数学定义是解析信号相位的导数，当 IMF 不满足 Bedrosian 定理（即信号带宽过宽、包络变化过快时），Hilbert 变换得到的相位信息失真，导数就会出现负值。解决：检查对应 IMF 的波形是否对称，如果某个 IMF 的对数上升沿和下降沿存在明显不对称，先对 IMF 做一次窄带滤波再求瞬时频率；另外在使用hhspectrum时把负频率直接置零是工程界的普遍做法。更进一步的做法是改用解调技术分析瞬时频率，用 Hilbert 变换得到的是近似解，精度不够时别硬撑。
6. 进阶技巧：三个验证步骤，把分解质量从“能出图”提升到“能信”
用合成信号验证分解正确性的方法值得这么做。构造一段已知瞬时频率解析式的信号，比如频率随时间线性扫的信号x(t) = sin(2π(100+50t)t)，分解后用hht求瞬时频率，和理论值 100+100t 对比。偏差在 1% 以内说明分解质量良好，偏差超过 5% 说明参数要调整或端点效应干扰太大。这个验证流程我每次都跑，因为它能把“算法内部黑匣子”变成“可量化校准的仪表”。
如果验证发现瞬时频率偏差偏大，常见修复路径是缩短分析区间，先把端点影响排除在外，看区间中段的误差是否收窄。中段误差也大，问题就不在端点上，而是分解本身混入了相邻成分。这时优先调SiftThreshold，把它从默认值向 1e-3 方向调整，降低筛分精度，反而能让 IMF 更符合窄带条件。两种方向都试过仍不理想，换掩膜信号往往能解决问题。
关于正交性指标，EMD 理论认为各个 IMF 之间应该近似正交，实际数值上不可能完全为 0。计算正交性指数的方法是：对 IMF 矩阵做两两内积归一化，取非对角元的绝对值均值。这个值小于 0.1 说明分解质量可接受，超过 0.3 就说明模态混叠严重，前几章的避坑策略该用上了。我习惯把这段指标计算代码存成一个独立脚本，每次换数据集都会跑一遍，作为筛选是否进入后续时频分析的“准入条件”。
我在这套工具箱上交过不少学费：最开始迷信默认参数，拿到故障信号就直接分解，画出来的时频谱花里胡哨，功率谱看什么都像特征频率；后来学会先用合成信号标定参数，再处理实际数据，效果立刻有质的提升。做 HHT 分析没有捷径，多花二十分钟做验证，能帮你省下后面整理结论时的几天返工。这些习惯希望能帮到你。
本文还有配套的精品资源，点击获取
