简介：一份基于希尔伯特-黄变换（HHT）的MATLAB时频图绘制脚本，针对非线性、非平稳信号分析场景而设计，面向信号处理、故障诊断、生物医学工程等领域的入门与进阶学习者，可帮助快速得到信号的时频分布。资源包仅含1个m文件，压缩后大小约1KB，代码量精简但处理流程完整，主要涵盖经验模态分解（EMD）、希尔伯特变换、瞬时频率与瞬时幅度计算，以及时频图的绘制与基础后处理。该资源已有441人学习下载。借助这个示例，读者可以掌握HHT时频图的完整生成链路：先通过EMD将原始信号自适应分解为若干个本征模态函数（IMF），再对每个IMF分量进行希尔伯特变换，计算其瞬时频率与幅度，然后组合成一张时间-频率-幅度的二维图像，直观观察频率随时间的动态变化。脚本结构清晰，适合作为配套代码模板学习算法原理，也便于替换输入数据，迁移应用到轴承故障诊断、心电信号分析、地震信号处理等场景，同时可结合参数调整与去噪步骤，为后续深入研究提供基础。

1. HHT时频图：拆开Untitled.m，看MATLAB里时频图产出全流程

实际场景是这样的：车间里测到一段滚动轴承振动信号，长度几十秒、采样率25600 Hz，你拿FFT做频谱，峰值一大堆，根本分不清哪一撮频率成分是启动时出现的、哪一撮是停机前消失的。频率成分随时间变化，这就是非平稳信号的痛点。HHT时频图（希尔伯特-黄变换时频图）把信号先做经验模态分解，把原始波拆成若干个瞬时频率和振幅都清晰的分量，再通过Hilbert变换映射到时间-频率平面，形成一幅动态谱图。这份.rar里的Untitled.m正是一个完整的MATLAB实现，从加载信号、EMD分解、Hilbert变换到绘图一步不缺，跑通一遍，你就有了一副自己的HHT时频图画法脚手架。

适合谁？做机械故障诊断的、碰过地震波和脑电、心电这类瞬态信号的人，都用得上。它教的不是调一个现成函数交差，而是让你理解每个环节的边界条件——哪个参数直接决定图能不能用，哪个步骤偷懒会画出漂亮但没物理意义的图。下文按真实复现顺序把代码拆开讲，并在第5章把最常见的翻车地点列出来，先有预期，后面操作会顺手很多。

2. 先跑通EMD分解：IMF不选对，后面的时频图全是五彩斑斓的噪声

2.1 EMD为什么是数据驱动的解耦器

第一次接触HHT的工程师会问：我要时频分析，为什么不用短时傅里叶变换（STFT）或连续小波变换（CWT）？这两种确实成熟，但前提是先选基函数和窗函数。STFT窗口一固定，频率分辨率就受窗长限制，窗开得窄就换不来时间与频率的双重精度；小波可以缩放，但小波函数本身仍是人为规定的。

EMD的特别之处在于不预设基函数，它根据信号自身的包络极值特性，一层层剥出有物理意义的分量，叫本征模态函数（IMF）。IMF满足两个约束：极值点数量与过零点数量相等或最多相差一个；上下包络均值在任意位置近似为零。实际执行时，EMD用三次样条把局部极大值点连成上包络、局部极小值点连成下包络，求均值后从原信号减掉，反复迭代直到满足停止条件。

听起来不复杂，但迭代次数和极值点选取对IMF形态影响很大，直接决定时频图质量：IMF残留了不该有的成分，瞬时频率就不干净，时频图中间会出现互相纠缠的伪能量带。

2.2 在MATLAB里跑EMD：一份能跑通的最小路径

% ===== 信号加载与预处理 ===== fs = 25600; % 采样率，单位Hz t = (0:N-1)/fs; % 时间轴 x = signal_data(:); % 确保列向量，经典踩坑点：行向量容易把 % emd结果矩阵的维度搞反 % 去均值和去趋势：EMD对直流偏移和线性趋势都敏感 x = x - mean(x); x = detrend(x); % ===== EMD分解 ===== % 新版MATLAB自带emd函数；旧版或这份Untitled.m里往往是自定义实现 imf = emd(x, 'MaxNumIMF', 8, 'Display', 0); % imf：每一列对应一个IMF，最后一列通常是残差 [nSamp, nImf] = size(imf); disp(['样本点数: ', num2str(nSamp), '，分解出IMF数量: ', num2str(nImf)]);

逻辑说明：emd返回的矩阵里，列数就是分解出的分量个数，最后一列是残差，它本质上是一条趋势项，不能当作振荡分量处理。MaxNumIMF限制最多分解出8个分量，对多数振动信号够用；Display设为0关闭日志，避免控制台刷屏。如果MATLAB版本没有内置emd，就把Untitled.m里自带的自定义函数路径用addpath加到搜索路径，接口保持一致即可。

detrend这一行值得单独说。原始信号若有直流偏移或缓慢漂移，EMD会把大趋势优先剥成第一个IMF，真正关注的瞬态分量反而被挤到后面。预处理不干净，后面所有IMF的物理含义都变成一张糊涂账。

2.3 丢弃残差与低频拖尾：选IMF的三个经验指标

EMD跑完，会得到若干个IMF，但画图时不是全都要用。我一般按三条指标过滤一遍：

% ===== IMF体检 ===== for k = 1:nImf r = corrcoef(x, imf(:,k)); varRatio = sum(imf(:,k).^2) / sum(x.^2); zeroCross = sum(abs(diff(sign(imf(:,k)))) > 0); fprintf('IMF%d: 相关系数%.3f, 方差占比%.3f, 过零点数%d\n', ... k, r(1,2), varRatio, zeroCross); end

参数含义与经验阈值：相关系数看IMF与原信号的重合程度，一般取0.1以上才有保留价值，低于这个值说明分量主要是噪声；方差占比看能量贡献，占比较低但相关性也低的IMF不要画进图里；过零点数代表振荡快慢，越往后过零点越少，末尾几个IMF的振荡周期已接近数据长度，物理意义趋近于趋势而非振动。

最后一列残差几乎肯定过零点极少，方差占比却可能不低。如果强行画进时频图，会在频率接近0的位置拉出横贯整张图的能量条，把有效信息全部淹没。经验做法是把残差排除，再看前几个高频IMF的频率轨迹是否清晰，逐个确认后再进入Hilbert变换环节。

提示：IMF体检不是走流程，每换一组数据都必须重跑一次。同一个算法在不同信号上分解结果差异极大。

3. 从IMF到瞬时频率：Hilbert变换的调用姿势与频率合法性验证

3.1 为什么不能对IMF直接做FFT当瞬时频率

IMF本身有一定带宽，直接FFT会得到一片频谱，得不到某个时间点上的频率值。Hilbert变换的思路是把实信号补成复数解析信号：解析信号z(t) = x(t) + j·H[x(t)]，其中H是希尔伯特变换，瞬时频率定义为相位φ(t)对时间的一阶导数再除以2π。

Hilbert变换是一种线性操作，对单分量信号效果最好，正好和EMD切分出的IMF组合成一个完整链路：EMD负责把多分量信号拆成单分量，Hilbert负责从单分量提取瞬时频率。两个步骤缺一个都立不住。直接拿原始信号做Hilbert变换，多组分频率互相调制，瞬时频率结果没有稳定含义。

3.2 逐IMF计算瞬时频率的代码实现

把上一章的IMF矩阵拿过来，逐列处理：

% ===== 对每个IMF做希尔伯特变换、求瞬时频率 ===== nUse = nImf - 1; % 排除残差，只处理前nUse个IMF H = zeros(nUse, nSamp); % 瞬时频率矩阵 A = zeros(nUse, nSamp); % 瞬时幅值矩阵 for k = 1:nUse xk = imf(:,k).'; zk = hilbert(xk); % 解析信号 phase = unwrap(angle(zk)); % 相位展开，避免±π跳变 freq = diff(phase) / (2*pi) * fs; % 数值微分使矩阵长度少1，把末尾值复制补齐 freq(end+1) = freq(end); H(k,:) = freq; A(k,:) = abs(zk); end

这里有三个细节必须注意。第一，相位展开：angle返回的相位被限制在实时间轴(-π,π]内，假如瞬时频率在采样间隔内跨过π边界，不展开就会在断点处产生跳跃，求导后得到一个极大的假频率尖峰。unwrap能自动加上2π的倍数让相位连续，但它要求相邻采样点间的相位差不能超过π，否则会判断不了方向，这时需要提高采样率或先对IMF做平滑。

第二，数值微分：diff结果是nSamp-1个点，直接赋值会破坏矩阵维度，所以要补一位。严格说末尾补值法会引入一个小误差，但只影响最后一个样本点，整图看不到；想要严谨可以改用gradient(phase)/(2*pi)*fs，它对端点用单边差分，中间用中心差分，整条曲线更平滑。

第三，平滑窗：直接对解析信号求导得到的频率曲线通常有高频毛刺，尤其在噪声段。我一般会用movmean(freq, 11)做11点滑动平均，窗参数根据采样率调整，采样率越高窗可以越大。但窗别开过宽，否则短时的真实频率跳变会被磨平，故障冲击特征就丢了。

3.3 频率合法性与负值毛刺处理

瞬时频率物理限制是必须在0到fs/2之间。噪声较大或端点震荡较强的区段，求导可能得到负频率，也可能超过奈奎斯特频率。画图前做一次合法性过滤：

% 非法频率统一置NaN，渲染时自动跳过 valid = (freq >= 0) & (freq <= fs/2); freq(~valid) = NaN;

这里必须讲透一个处理原则：负瞬时频率不代表真实存在的反向旋转，更多是相位噪声在低幅值段落里主导了求导结果。单纯用abs(频率)把负值掰成正值，等于在时频图上伪造能量，噪声区的假频率会被抬进有效频段。正确顺序是先看瞬时幅值，把幅值低于某个阈值的区域整体置NaN，再过滤频率合法性。哪个阈值合适？一般取该IMF幅值均值的10%，低于这个水平时信号局部接近消失，瞬时频率没有分析价值。

实际操作中还可以在hilbert之前对IMF做一次带通滤波，滤掉噪声成分，Hilbert变换会稳定很多。但滤波本身又会改变相位，所以要控制滤波器阶数和通带宽度，尽量避免滤波引入相移。这个平衡算是HHT实践里的玄学部分，只能根据具体信号反复试。

4. 把瞬时频率画成HHT时频图：imagesc、坐标轴与颜色映射的四个参数细节

4.1 从逐点频率行到可渲染的时频矩阵

到这一步，H矩阵是nUse×nSamp，第k行对应第k个IMF，每一列对应同一个时间点。可视化的本质是把这些点按频率拼到时间-频率平面。这里最容易困惑的是：每个IMF在某个时刻只能给出一个瞬时频率和一个瞬时幅值，所以HHT时频图不是传统意义上的二维谱图，更像是一张三维分布投影。

如果用散点一个个画，几千个点叠在一起，图面会非常乱，也看不出能量强度。所以我通常用imagesc渲染：行索引映射为频率，列索引映射为时间，颜色深浅代表该时刻该频率上的能量强弱。

4.2 渲染代码与参数说明

% ===== 构建可渲染矩阵并绘图 ===== figure('Color', 'w'); freqAxis = linspace(0, fs/2, nUse); imagesc(t, freqAxis, A); axis xy; % 让纵轴从0向上增长，否则频率轴倒置 colormap(parula(256)); c = colorbar; c.Label.String = '瞬时幅值'; xlabel('时间 (s)'); ylabel('频率 (Hz)'); ylim([0 2000]); % 只显示0~2000 Hz观察窗，按信号实际调

参数说明：imagesc的Y轴是把矩阵行索引等分映射到第二个参数给的范围，我这里用linspace(0, fs/2, nUse)构造频率轴，保证行号与赫兹值一一对应。如果直接传[0 fs/2]，MATLAB只把它当成坐标范围的上下限，内部仍按行数线性均分，效果一样，但后者在你想做频率轴非线性映射时会失灵，建议用前者。

axis xy几乎必写。MATLAB里imagesc默认按图像行列方式渲染，行索引向下增长，也就是频率轴从高到低排列，不纠正就是一张倒谱图。parula比jet更适合做物理量渲染，jet在颜色阶跃处会产生视觉上的伪轮廓，容易把连续频率误读成边界。幅值动态范围很大时，建议改成imagesc(t, freqAxis, log10(A + eps))，这样小幅值故障特征也能看见。

4.3 后处理：裁剪边缘和动态范围

端点效应是HHT的经典坑，三次样条包络在数据首尾缺少足够的极值点支撑，包络均值会外飘，对应IMF在端点变形，时频图两端出现明显云团。最省力的办法是绘图前裁剪首尾5%的时间点：

% ===== 裁剪边界，削弱端点效应 ===== cut = round(nSamp * 0.05); t_plot = t(cut:nSamp-cut); A_plot = A(:, cut:nSamp-cut);

动态范围方面，幅值矩阵A里最大值和最小值可能差几个数量级。机械故障诊断里，早期故障的特征幅值很小，和正常振动混在一起，线性色标下几乎看不见。用log10(A+eps)做对数压缩，能把微小冲击从背景里扒出来。这个技巧在轴承外圈故障、齿轮点蚀这类早期微弱信号上非常管用。

提示：eps的作用是防止对数里出现0值。也可以改成log10(max(A(:))/1000 + A)，动态范围压缩得更深。

5. HHT时频图画完对不上直觉：四个常见坑的排查与修复

5.1 端点效应：图两头出现钩形色块

现象：时频图在时间轴开头和末尾出现明显钩子，频率在端点异常爬高或跌落，色块显示时两端各弯出一条能量带，中段反而干净。

原因：EMD的包络拟合在边界处缺少有效极值约束，三次样条外推让包络均值在首末两端产生较大误差，对应的IMF在端点变形，Hilbert变换求出的瞬时频率随之异常。

解决：第一，做EMD前采用镜像延拓把数据两端补一段，分解完成再截断；第二，绘图时统一裁剪首尾5%的时间窗，这是最省力的做法，上一章代码可直接用；第三，如果端点区域恰是你关心的阶段，比如设备启动过程，就不能靠裁剪，要回到原始信号检查端点处是否真有突变，必要时先做带通滤波再分解。

5.2 模态混叠：时频图里频率线拧成麻花

现象：本该在时频图上表现为一条平滑连续线的频率成分，断断续续折叠跳到别的频段，图上两条线像拧麻花，无法分辨多个分量。

原因：信号中存在间歇性成分，比如周期性冲击夹杂在里面，EMD的筛分过程常常把某个高频成分的一部分分给另一个IMF，导致分解结果不符合单分量假设。强噪声也会干扰极值点定位，加重混叠。

解决：第一步先对原始信号做带通滤波，把关注频段之外的成分压掉，再进行EMD。第二步改用EEMD或CEEMDAN，通过多次加白噪声集成平均抑制模态混叠。MATLAB新版本里eemd或ceemdan的调用方式和emd基本一致，主要参数是噪声标准差Nstd和集成次数NE，常见取法为Nstd=0.15、NE=200。参数改完后IMF数量和形态都会变，记得重新做第2.3节的IMF体检，不要拿旧表直接套。

5.3 负瞬时频率与全图低频轰炸

现象：时频图低频区出现一大片连续深色能量团，把所有细节淹没，查看瞬时频率矩阵时发现里面有不少负值。

原因：可能有两条。一是原始数据包含强趋势，趋势项被EMD分离成一个极低频IMF，画图时占据低频区；二是某个时刻信号幅值很小、噪声占主导，相位求导结果在噪声驱动下乱跳，负频出现是常态。

解决：低频能量团先排除残差，再看第一个IMF是否本身就是趋势项，必要时加深detrend阶数或对信号做高通预滤波。负瞬时频率的推荐处理是第3.3节给出的合法性过滤，配合幅值阈值。在这里再强调一遍：不要用abs强行把负频率掰成正，负频率出现的区域往往对应信号真空区，硬掰等于在图上伪造一组不该存在的能量。

5.4 与FFT对比总是对不上号

现象：同一段信号用FFT谱峰是200Hz，HHT时频图上的频率线却经常飘在240Hz附近，图里颜色最亮的位置和谱峰对不上，让人怀疑算法有误。

原因：FFT是全场平均频谱，对时长内的能量做加权平均，谱峰是时间的统计结果；HHT时频图展示的是瞬时频率，信号频率本身在漂移，图上每个时刻的线都是瞬时值。若信号从150Hz扫到250Hz，FFT谱峰落在频率轨迹的能量重心处，不会正好等于200Hz。

解决：不要拿整段FFT去对整张时频图。截取一小段，比如信号中段200个采样点单独做FFT，与时频图对应时段内的平均瞬时频率对比，误差在5%以内说明流程没问题。若偏差显著，回到5.2节排查模态混叠。

6. 拿合成chirp做一次自检：验证时频图不是伪影的有效方法

最终步骤强烈建议你用一段已知答案的信号，从加载到出图把整个HHT流程端到端验一遍。我的习惯动作是生成一个从低频扫到高频的chirp，叠加一个固定频率的正弦分量，两者同时送入管线。chirp的瞬时频率有精确解析解，理论上时频图上应出现一条从起点频率线性爬到终点频率的连续亮线，另有一条水平穿过整个时间轴的恒定频率亮线。两条线交叠处会有轻微干涉，但不会互相吞没。

% ===== 合成信号自检 ===== fs = 1000; t = 0:1/fs:5; x = 0.6 * chirp(t, 20, 5, 200) + 0.4 * sin(2*pi*85*t); % 之后把第2章到第4章的EMD->Hilbert->绘图流程原样套用

跑出来的图，正常情况这两条轨迹清晰独立。如果只有一条线亮，大概率是分解环节把另一路能量弄丢了，回IMF体检表查幅值占比；如果两条线的频率与理论值差超过百分之几，去查Hilbert变换和采样率换算；如果曲线头尾有异常弯曲，说明端点效应清理得不够彻底。这一测，等于给整套流水线做了个冒烟测试，哪一环掉链子立刻暴露。

从那以后我每次拿到新数据，会先把同一个合成chirp测试信号跑一遍，确认时频图与理论频率轨迹对得上，再放心处理真实信号。时间久了，这个动作从锦上添花变成了例行自检。多花几秒钟，却能把一堆时频图玄学挡在门外。希望帮到你。

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