简介这是一份希尔伯特黄变换HHT完整 MATLAB 源码包面向学习非线性、非平稳信号处理的科研人员、工程师及高年级本科生可解决 EMD 分解、瞬时频率提取与希尔伯特谱绘制等核心问题。压缩包总体积 22KB共 17 个 m 文件覆盖基础经验模式分解、针对噪声/局部/在线场景的改进版本、希尔伯特谱分析、可视化与辅助处理等模块目录结构清晰便于按需调用或二次开发。代码中既包含完整 EMD 分解流程与 IMF 提取逻辑也有瞬时频率/幅度计算、时频谱图显示及合成测试信号生成等配套功能可衔接从信号预处理到 HHT 全流程分析便于理解算法原理和拓展实验。已有 3894 人学习下载资料虽小但功能完整适合在 MATLAB 中直接运行、研读算法实现也可作为教学演示、课程设计或项目改造的起点。1. HHT不是又一个FFT它告诉你频率在哪个时刻出现做频谱分析的人第一次接触希尔伯特黄变换HHT时常会有个错觉它像FFT一样输入一段信号输出一张频谱图。实际跑起来才发现它要先经过一次经验模态分解EMD把信号切成若干个本征模态函数IMF再对每个IMF算瞬时频率。这套流程对齿轮点蚀、心电干扰、地震波记录这类频率随时间变化的非线性非平稳信号分辨方式比短时傅里叶更贴近直觉频率不是一排固定的柱状图而是一条随时间连续移动的线。这也是HHT程序的结构特点——它不是单一函数而是由EMD分解、Hilbert变换、时频谱绘制三块串起来的一条计算链。下面我把这条计算链写成可直接运行的MATLAB程序自写EMD分解、希尔伯特谱和边际谱一次给全并补上停止准则、端点效应、插值方法这几个容易把结果带偏的参数让初学的研究生和想从工程上替换短时傅里叶的工程师都能直接照用。2. 从EMD筛分到瞬时频率写MATLAB之前先把算法链默写一遍写代码前我习惯先在一行注释里写下完整的计算链x(t) → EMD → IMF1..IMFn residual → Hilbert变换 → A(t)、f(t) → H(ω,t) → h(ω)。后续程序的所有分支本质上都是这条链的实例化。HHT这个词在实际使用里经常被简写成“EMDHilbert谱”原因就在这里。2.1 EMD筛分筛选IMF的迭代条件不是“一次减均值”对非平稳信号EMD不预设基函数而是靠极值点构造包络。筛分循环的经典步骤如下找出信号所有局部极大值和局部极小值。用三次样条把极大值点连成上包络u(t)极小值点连成下包络l(t)。计算包络均值m(t) (u(t)l(t))/2。从当前信号中减去m(t)得到候选分量。对新候选分量重复第14步直到满足IMF定义保存该IMF。原信号减去该IMF对剩余信号继续下一轮筛分。一个分量满足IMF条件通常要同时看两点极值点个数与过零点个数相等或至多相差1上下包络的平均值趋近0。换句话说IMF是局部窄带的每个瞬时时刻只有一个主导频率。这也是后面求瞬时频率的前提宽带信号直接求相位导数会得到一堆没有物理意义的负频率。筛分过程需要结束条件经典判据是标准差SD。SD的计算式通常是相邻两次筛分结果之差的能量除以当前信号能量。SD越小筛分越精细但并不意味着越好——筛分次数过多会把有物理意义的幅度波动抹平让IMF变成等幅正弦。后面第四章会专门讲这个阈值的设置范围这里先记住筛分是迭代不是一步减法。2.2 Hilbert变换瞬时频率是相位导数不是能量重心单个IMF记作c(t)。对它做Hilbert变换得到H{c(t)}后构造解析信号z(t) c(t) j·H{c(t)} A(t)·e^(jθ(t))其中A(t)是瞬时幅值θ(t)是瞬时相位瞬时频率定义为相位对时间的导数f(t) (1/2π)·dθ(t)/dt。这里必须用unwrap解卷绕否则相位在±π处跳变求导后会得到巨大脉冲。对应到MATLAB核心就是三行z hilbert(imf_k); % 解析信号 inst_amp abs(z); % 瞬时幅值 inst_freq fs/(2*pi) * [0; diff(unwrap(angle(z)))];最后补的0是为了让inst_freq和原信号等长画图或做逐点累加时不用来回截断。瞬时频率出现负值通常说明该IMF不纯混入了其他分量或者端点处样条过冲严重。2.3 一套完整程序最终输出三样东西程序跑完应该同时保留三类结果。这三类结果对应三种不同的工程用途。输出对象怎么得到典型用途IMF和残差EMD筛分循环分离不同时间尺度的分量残差是趋势项希尔伯特谱H(ω,t)每个时刻的瞬时频率和瞬时幅值观察特征频率随时间的变化如转速波动边际谱h(ω)希尔伯特谱对时间积分对特征频率的存在强度做整体判断边际谱不是功率谱密度它和FFT幅值谱的量纲、分辨率都不相同不能直接拿峰值大小去和FFT对比。它回答的问题是“整个时间段里哪些频率成分出现过、累计能量有多强”。理解了这三样东西下面代码每个变量该存什么就清楚了。3. 可运行的HHT完整MATLAB程序EMD、希尔伯特谱和边际谱一次跑通3.1 先构造一个能自检的测试信号不要拿真实数据调试算法。我一般先构造两个非谐波关系的正弦加白噪声比如50Hz和120Hz频率不呈整数倍关系避免分解结果被采样周期巧合干扰。fs 1000; % 采样率 1000 Hz t (0:500-1)/fs; % 0.5 秒数据 x 1.2*sin(2*pi*50*t) 0.8*cos(2*pi*120*t); x x 0.2*randn(size(t)); % 加白噪声观察稳定性50Hz和120Hz落在低频段500个采样点足以让三次样条构造出稳定的包络。白噪声幅度取信号标准差的10%左右太小看不出抗噪效果太大会让高频IMF被噪声主导。这段代码跑通后再换成你自己的传感器数据。3.2 自写EMD函数筛分循环和端点镜像一次拿全先说明一点MATLAB R2018a之后的Signal Processing Toolbox内置了emd和hht工程交付时我通常直接用内置函数因为速度和边界处理更稳定。但调试算法、理解每个IMF怎么筛出来、给论文方法部分画流程图时自写版本不可替代。下面是完整可运行的emd_self.m。function [imf, residual] emd_self(x, Sd, MaxIter) % emd_self 一个可直接读懂的EMD实现 % x: 列向量信号 % Sd: 筛分停止阈值经典取值 0.2 ~ 0.3 % MaxIter: 单个IMF的最大筛分次数防止死循环 x x(:); if nargin 3, MaxIter 200; end if nargin 2, Sd 0.2; end imf zeros(numel(x), 0); % 每列一个IMF r x; while true h r; % 内循环反复减去包络均值直到满足停止条件 for iter 1:MaxIter m envelope_mean(h); h_new h - m; sd sum((h - h_new).^2) / (sum(h.^2) eps); h h_new; if sd Sd break; end end imf [imf, h]; % 保存一个IMF r r - h; % 残差继续分解 % 残差极值点不足2个时认为只剩趋势项 if sum(~isnan(findpeaks(r))) 2 ... sum(~isnan(findpeaks(-r))) 2 break; end end residual r; endfindpeaks来自 Signal Processing Toolbox。如果你只有基础版MATLAB可以用diff(sign(diff(r)))手工找极值点效果等同只是代码长一点。envelope_mean是核心子函数负责构造上下包络并返回均值。function mval envelope_mean(h) % 对信号h分别做上下包络三次样条拟合返回包络均值 [pks_max, loc_max] findpeaks(h); [pks_min, loc_min] findpeaks(-h); if numel(loc_max) 2 || numel(loc_min) 2 mval zeros(size(h)); % 极值点不够不再修正 return; end pks_min -pks_min; % 端点补极值把首尾样本当作极值点减少样条外插发散 loc_max [1; loc_max(:); numel(h)]; pks_max [h(1); pks_max(:); h(end)]; loc_min [1; loc_min(:); numel(h)]; pks_min [h(1); pks_min(:); h(end)]; up spline(loc_max, pks_max, 1:numel(h)); down spline(loc_min, pks_min, 1:numel(h)); mval (up(:) down(:)) / 2; end端点补极值是自写EMD里最值得留意的处理。直接对原始信号做样条拟合首尾附近的包络会因为缺少极值约束而飞出去产生俗称的“飞翼”现象污染物体现在第一个和最后几个IMF上。这里把端点当作极值点参与拟合是一种简单有效的抑制方式第四章会讲更完整的镜像延拓。3.3 对每个IMF做Hilbert变换瞬时幅值和瞬时频率同时拿到主程序里调用完emd_self立刻进入Hilbert变换循环。我习惯用cell数组分别存每个IMF的幅值和频率这样后续画图时按IMF逐条处理不需要再次索引原始矩阵。A cell(1, size(imf, 2)); f_inst cell(1, size(imf, 2)); for k 1:size(imf, 2) z hilbert(imf(:,k)); A{k} abs(z); phase unwrap(angle(z)); f_inst{k} fs/(2*pi) * [diff(phase); 0]; % Hz endhilbert是MATLAB自带函数返回解析信号不需要额外工具箱。瞬时频率数组最后补一个0是为了让f_inst{k}和imf(:,k)长度一致。调试时我会统计每个IMF里负频率的个数如果中间数据段负频率占比超过10%基本可以判断分解不彻底或信号本身不适合直接做HHT。3.4 希尔伯特谱矩阵和边际谱用频率分箱累加希尔伯特谱本质上是把每个时刻的瞬时频率映射到频率轴上把对应瞬时幅值累加进去。分辨率由f_res控制这里取0.5Hz频率轴从0到fs/2。f_res 0.5; f_axis 0:f_res:fs/2; H_spec zeros(numel(t), numel(f_axis)); for k 1:size(imf, 2) idx round(f_inst{k} / f_res) 1; keep idx 1 idx numel(f_axis); for n find(keep(:)) H_spec(n, idx(n)) H_spec(n, idx(n)) A{k}(n); end end figure; imagesc(t, f_axis, H_spec.); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); colorbar; marginal sum(H_spec, 1); % 对时间积分得到边际谱 figure; plot(f_axis, marginal); xlabel(频率 (Hz)); ylabel(边际谱幅值);分箱时丢弃了负频率和超出频率轴范围的点同时把端点飞翼造成的异常瞬时频率挡在了谱图外。注释里的逻辑是每个样本点有一个瞬时频率它落在哪个频段就把该点的幅值加到那个频段上。采样率越高、数据越长这个矩阵就越大处理长数据时可以把f_res放大到1Hz或2Hz牺牲频率分辨率换取内存。4. HHT程序里最该调的3组参数停止准则、样条插值和端点延伸4.1 SD阈值0.2到底改不改emd_self里的Sd默认取0.2这是经典论文里的经验值但不是金科玉律。SD的计算公式在2.1里已经给出它衡量的是相邻两次筛分的差异。实际调参时我按下面这张表来试参数常用范围调参效果什么时候动它SD0.1 ~ 0.3越小IMF越接近纯窄带但筛分次数大幅增加边际谱峰值分不开时降到0.1MaxIter100 ~ 500防止个别数据段不收敛导致死循环遇到冲击信号或大幅值突变时调大IMF个数不直接设置由残差极值点数量决定分解结果太多需要检查是否过筛SD过小的典型特征是某个IMF出现等幅正弦形状完全失去原始信号的幅值波动这就是“过度筛分”。反过来SD太大时第一个IMF还混着另一个分量的残余边际谱上会有拖尾峰。调试时我通常在循环里记录每次迭代的sd值打印出来看它在阈值附近是缓慢逼近还是震荡。4.2 端点镜像延拓样条在端点外插一定会飞翼3.2 的envelope_mean里用的是“端点当作极值点”这是最简处理。更标准的做法是镜像延拓以左右端点为对称轴把信号向外翻转一段长度翻转后的极值点也参与样条拟合拟合完成后再把延拓部分裁掉。function ym mirror_extension(y, n) % 以左右端点为对称轴各镜像n个点 ym [y(end:-1:end-n1); y; y(n:-1:1)]; end调用时把ym送进emd_self分解完只取中间的n1:end-n段。延拓长度n一般取信号长度的5%10%或者直接取一个IMF内平均极值间隔的两倍。需要记住的是无论怎么做端点处理瞬时频率的头尾十几个点仍然不可信画希尔伯特谱时裁掉边界段是常规操作。4.3 spline和pchip各管一段短数据别迷信三次样条spline插值在极值点间距大的地方容易过冲产生幅度夸张的包络进而把IMF的瞬时幅值抬高。对长度只有几百点、且包含冲击成分的数据我经常换成pchip。改动只有一行把envelope_mean里的spline改成pchip。pchip不会产生过冲但包络的平滑度差一些。经验判断是数据干净、极值点密集用spline数据带噪声或采样点稀缺用pchip。如果两种插值的分解结果差异很大说明原始信号本身不适合直接做EMD应该先滤波或重采样。4.4 模态混叠是EMD的短板EEMD加白噪声绕过去模态混叠是EMD被吐槽最多的问题某个IMF可能在一段时间内含有高频分量过一段时间又只剩低频分量频率发生分段跳变。解决思路是在信号里加入白噪声多次分解后取平均利用噪声的统计均匀性把不同尺度的分量分离这就是EEMD。function [imf, residual] eemd_self(x, EnsNum, NoiseAmp) % 集合经验模态分解仅演示核心逻辑 x x(:); imf []; residual 0; for e 1:EnsNum xi x NoiseAmp * randn(size(x)); [imf_i, res_i] emd_self(xi, 0.2, 300); if e 1 imf imf_i; else ncol min(size(imf, 2), size(imf_i, 2)); imf imf(:,1:ncol) imf_i(:,1:ncol); end residual residual res_i; end imf imf / EnsNum; residual residual / EnsNum; end实际使用时噪声幅值取原始信号标准差的0.10.4倍集合次数取50200次。EnsNum越大结果越稳定但计算时间线性增长。如果不想自己维护不同分解次数之间的IMF对齐逻辑工程上可以直接用内置emd配合循环实现或者接现成的EEMD实现做结果对比。噪声幅值过小起不到抑制混叠的作用过大会把低频分量淹没这是最常见的失败原因。5. 对真实数据动手CSV导入、结果校验和自动提取边际谱峰值频率5.1 外部数据导入别把带直流和缺失值的数据直接喂给EMD真实传感器数据的第一步永远是清洗。以最常见的CSV文件为例我通常用readtable读入后先做三个动作再去均值、去趋势、补NaN。不做这些预处理EMD会把直流分量和线性趋势当成一个IMF输出白白消耗分解层数。T readtable(demo_data.csv); t T.Time; x T.Signal; x x - mean(x); % 去直流 x detrend(x); % 去线性趋势 % 缺失值用线性插值补齐 x(isnan(x)) interp1(t(~isnan(x)), x(~isnan(x)), t(isnan(x))); [imf, residual] emd_self(x, 0.2, 300);detrend对非线性趋势无能为力如果数据有明显弯曲先用平滑滤波把趋势估计出来再减掉。缺失值超过总长度5%的数据不建议直接插值因为样条包络会把缺失段的形态猜错导致IMF出现虚假振荡。5.2 结果是否可靠的3个检查分解完成后不要急着画图。我先跑一段检查命令确认三个事实残差的标准差远小于原始信号标准差且残差里不应该还看得出周期波动每个IMF的均值都接近0瞬时频率的负值点应该只集中在头尾边界段中间数据段出现大量负频说明IMF不纯。for k 1:size(imf, 2) fprintf(IMF%d mean%.6f negFreq%d/%d\n, ... k, mean(imf(:,k)), sum(f_inst{k}0), numel(f_inst{k})); end边界段的负频率可以容忍但中间数据段如果出现超过5%的负频点就要回到第四章调SD和插值方式。把检验命令固化在脚本里每次换数据都执行一遍能省去大量排查时间。5.3 自动提取边际谱峰值频率直接作为后续特征边际谱画出来之后从图上读频率是低效的。我习惯用findpeaks直接从marginal里提取峰值频率列表。注意这里marginal的横轴是f_axis峰值的位置要映射回Hz。[mag, loc] findpeaks(marginal, MinPeakHeight, 0.1*max(marginal), ... MinPeakDistance, round(2/f_res)); peakFreq f_axis(loc); % 输出到CSV方便后续特征使用 writetable(table(peakFreq(:), mag(:), VariableNames, {Freq_Hz, Mag}), ... hht_peaks.csv);MinPeakHeight取最大边际谱值的10%用来过滤噪声产生的矮峰MinPeakDistance取2Hz避免同一个真实特征峰被拆成多个相邻峰。如果故障特征频率是边频带结构比如齿轮故障的调制边带可以把MinPeakDistance调小到1Hz然后把peakFreq两两作差差频对应的就是调制频率。这份峰值表后续做阈值告警、聚类或者送给深度学习网络都不用再重跑一次HHT。本文还有配套的精品资源点击获取