★ FEATURED ARTICLE 简介一份基于希尔伯特-黄变换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(end1) 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(Aeps)做对数压缩能把微小冲击从背景里扒出来。这个技巧在轴承外圈故障、齿轮点蚀这类早期微弱信号上非常管用。提示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常见取法为Nstd0.15、NE200。参数改完后IMF数量和形态都会变记得重新做第2.3节的IMF体检不要拿旧表直接套。5.3 负瞬时频率与全图低频轰炸现象时频图低频区出现一大片连续深色能量团把所有细节淹没查看瞬时频率矩阵时发现里面有不少负值。原因可能有两条。一是原始数据包含强趋势趋势项被EMD分离成一个极低频IMF画图时占据低频区二是某个时刻信号幅值很小、噪声占主导相位求导结果在噪声驱动下乱跳负频出现是常态。解决低频能量团先排除残差再看第一个IMF是否本身就是趋势项必要时加深detrend阶数或对信号做高通预滤波。负瞬时频率的推荐处理是第3.3节给出的合法性过滤配合幅值阈值。在这里再强调一遍不要用abs强行把负频率掰成正负频率出现的区域往往对应信号真空区硬掰等于在图上伪造一组不该存在的能量。5.4 与FFT对比总是对不上号现象同一段信号用FFT谱峰是200HzHHT时频图上的频率线却经常飘在240Hz附近图里颜色最亮的位置和谱峰对不上让人怀疑算法有误。原因FFT是全场平均频谱对时长内的能量做加权平均谱峰是时间的统计结果HHT时频图展示的是瞬时频率信号频率本身在漂移图上每个时刻的线都是瞬时值。若信号从150Hz扫到250HzFFT谱峰落在频率轨迹的能量重心处不会正好等于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测试信号跑一遍确认时频图与理论频率轨迹对得上再放心处理真实信号。时间久了这个动作从锦上添花变成了例行自检。多花几秒钟却能把一堆时频图玄学挡在门外。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
