简介经验模态分解与希尔伯特-黄变换的MATLAB实现程序包面向信号处理、设备故障诊断、振动分析及非平稳数据分析等领域的研究人员和工程师。该程序利用经验模态分解将原始信号分解为多个本征模态函数进而通过希尔伯特变换获得瞬时频率与瞬时幅值有效解决非平稳信号难以直接处理的问题分解后的分量可叠加重构出平稳信号便于后续特征提取与趋势分析。压缩包结构紧凑仅含三个M脚本文件包括主程序、频谱计算函数与瞬时频率估计函数整体体积只有2KB非常轻量。使用时需配合已安装的经验模态分解工具箱中的分解函数调用适合具备一定MATLAB基础、希望快速掌握该方法或根据需求修改算法的学习者。目前已有一千二百七十七人学习下载代码逻辑清晰、注释完整可作为课程设计、科研实验或工程应用的参考模板帮助读者深入理解希尔伯特-黄变换的全流程。1. 非平稳信号的时间-频率刻度EMD与Hilbert-Huang变换各负责一半Hilbert-Huang变换HHT不是某一种变换公式而是一套流程先用经验模态分解EMD把非平稳信号拆成若干固有模态函数IMF再逐一对IMF做Hilbert变换求瞬时频率和瞬时幅值。拆这一步针对的是“频率随时间变”——傅里叶谱只能告诉你有哪几个频率告诉不了50Hz出现在第几秒直接对原始信号算Hilbert瞬时频率又会因多分量叠加而出现负频率和伪波动。HHT把两步分开恰好补上这两块短板。这套MATLAB程序的主程序是HHT.m压缩包里的hhspectrum.m和instfreq.m负责Hilbert谱参数计算EMD分解依赖已安装的EMD工具箱中的emd函数。适合做振动、故障诊断、脑电和风速分析的工程师拿到后改一下采样率和文件路径就能接进自己的数据。2. EMD筛分原理先有IMF才有可信的瞬时频率2.1 为什么直接算Hilbert瞬时频率会失败Hilbert变换把实信号x(t)映射为解析信号z(t)x(t)jH[x(t)]瞬时频率定义为相位对时间的导数。这个定义对单分量窄带信号成立一个调频正弦的相位导数就是直觉上的频率但多分量信号叠加后解析信号的相位在分量交叠处快速摆动导数频繁出现负值。负频率不是噪声是数学定义对“多分量”失效的信号。EMD的作用就是把x(t)拆成一组单分量IMF。IMF有两个判定条件一是整个数据段内极值点个数与过零点个数相等或最多差1二是任意时刻上包络和下包络的均值接近0即局部对称。实际判断时第二个条件靠迭代逼近不会严格等于0工具箱通过相对容差控制。默认参数下分解出的IMF已经满足工程需求。注意IMF不是正弦它是幅值和频率都被调制的窄带分量正因为允许包络缓慢变化Hilbert谱才能画出随时间变化的能量分布这一点和傅里叶分解有本质区别。2.2 筛分迭代与残差分解的停止条件决定分量个数EMD筛分过程按下面几步循环找出原信号的全部局部极大值和局部极小值点用三次样条分别拟合上包络u(t)与下包络l(t)计算包络均值m(t)(u(t)l(t))/2从信号中减去得到h(t)x(t)-m(t)对h(t)重复第1到第3步直到h(t)满足IMF两个条件得到第一个IMF原信号减去IMF1得到残差r1对r1继续执行第1到第4步得到IMF2如此往复直到残差单调或幅度低于阈值。每次减去包络均值的过程叫筛分sifting它把高频振荡一层层剥出来。筛分次数过多IMF会变成纯调频信号丧失物理意义次数过少IMF又不满足窄带条件。EMD工具箱在这个问题上用标准差阈值做了折中。调用方式统一如下imf emd(x); [npts, ncol] size(imf); num_imf ncol - 1; % 第三方工具箱的最后一列是残差如果用的是MATLAB R2018a之后自带的emd函数返回值是两段[imf, residual] emd(x)需要手动拼成imf [imf, residual]否则后续hhspectrum会少算一个趋势分量。这两种返回格式差异是新手最容易翻车的地方后面的排错表里会再提。2.3 EMD工具箱依赖与路径配置压缩包内没有emd.m主程序HHT.m第一步就会调用它。先确认环境是否就绪which emd返回路径说明工具箱可用。如果返回空把EMD工具箱目录加入MATLAB路径并持久化addpath(D:\Tools\EMD); savepath;这里不推荐在脚本里裸调addpath因为重启MATLAB后设置会丢失。换机器时直接执行一遍这两行或在startup.m里统一维护。确认工具箱可用的判断标准是emd(x)能对任意长度向量返回列数大于等于1的矩阵。与短时傅里叶、小波对比EMD最大的差异在于基函数来自数据本身不预设窗长和母小波方法基函数时频分辨率自适应能力短时傅里叶固定窗正弦窗长决定时间与频率分辨率互相制约无连续小波母小波伸缩平移多分辨仍依赖小波基与尺度步长部分EMDHHT数据驱动IMF瞬时频率逐点定义无窗长约束强这一特性决定了HHT适合频率成分连续变化、且变化速率不固定的信号。代价是分解结果对噪声有一定敏感性轻微扰动可能改变IMF个数所以后面的验证环节必须做。3. HHT.m、hhspectrum.m与instfreq.m的调用链与参数约定3.1 主程序HHT.m先分解再逐列变换主程序把整条流水线收口成一次调用。我习惯把采样频率作为必传参数因为后面画Hilbert谱和边际谱时都要把归一化频率换算成物理频率function [imf, A, f, tt] HHT(x, fs, draw) % HHT 主程序EMD分解 Hilbert谱参数计算 % 输入 % x - 待分析信号行向量或列向量均可 % fs - 采样频率单位Hz % draw - 是否绘制Hilbert谱散点图默认1 % 输出 % imf - IMF矩阵每列一个分量最后一列为残差 % A - 瞬时幅值矩阵时间×分量 % f - 瞬时频率矩阵Hz时间×分量 % tt - 时间轴s if nargin 2 error(至少需要输入信号x和采样频率fs); end if nargin 3 draw 1; end x x(:).; % 统一成行向量 N length(x); t (0:N-1) / fs; % 物理时间轴单位秒 imf emd(x); % 第1步EMD分解 [A, fn, tt] hhspectrum(imf, t); % 第2步瞬时幅值与归一化频率 f fn * fs; % 第3步归一化频率换算为Hz if draw figure; scatter(tt, f(:), 3, A(:), filled); xlabel(时间 (s)); ylabel(频率 (Hz)); colorbar; title(Hilbert-Huang 谱散点密度即能量分布); end这里用scatter做时间-频率-幅值的三维映射比imagesc更适合未经toimage网格化的原始输出因为瞬时频率本身不是均匀网格。emd(x)接受行向量或列向量但hhspectrum内部对行数、列数的判定依赖输入约定所以开头统一成行向量能省掉后面大部分维度报错。3.2 hhspectrum.m逐列Hilbert输出幅值和频率两个矩阵hhspectrum的职责是遍历imf矩阵的每一列对每个IMF做Hilbert变换得到瞬时幅值和瞬时频率function [A, f, tt] hhspectrum(imf, t, l) % HHSPECTRUM 计算每个IMF的瞬时幅值与瞬时频率 % 输入 % imf - EMD分解结果行数为采样点数列数为分量数 % t - 时间轴长度必须等于size(imf,1) % l - 频率平滑点数默认1 % 输出 % A - 瞬时幅值矩阵尺寸(N-1)×M % f - 归一化瞬时频率矩阵单位cycles/sample % tt - 与矩阵行对齐的时间轴长度N-1 if nargin 2 t 1:size(imf, 1); end if nargin 3 l 1; end if size(imf, 1) size(imf, 2) imf imf.; % 统一为每列一个IMF end M size(imf, 2); L size(imf, 1); A zeros(L, M); f zeros(L, M); for k 1:M an hilbert(imf(:, k)); % 解析信号 A(:, k) abs(an); % 瞬时幅值即解析信号模 fn instfreq(imf(:, k), t, l); % 归一化瞬时频率 f(1:length(fn), k) fn(:); % 频率序列比原信号短1点 end A A(1:end-1, :); % 尾部对齐去掉最后一行 tt t(1:end-1);瞬时频率用相邻两点相位差定义天然比原信号少1个点所以这里统一截掉最后一行让三个输出维度对齐。注意instfreq在内部会再做一次Hilbert变换这里直接传原始IMF列即可不需要把an传进去造成重复处理。hilbert是MATLAB自带函数不依赖额外工具箱。3.3 instfreq.m相位差分与平滑参数lfunction [fn, t] instfreq(x, t, l) % INSTFREQ 相邻样本相位差法计算瞬时频率 % 输入 % x - 实信号序列内部自动做Hilbert变换 % t - 时间轴 % l - 滑动平均窗口长度默认1 % 输出 % fn - 归一化瞬时频率单位cycles/sample % t - 对齐后的时间轴 if nargin 2 t 1:length(x); end if nargin 3 l 1; end x x(:).; t t(:).; N length(x); y hilbert(x); % angle(y(n)*conj(y(n-1))) 是相邻样本的相位差 % 除以2*pi换算为周期数即cycles/sample fn angle(y(2:N) .* conj(y(1:N-1))) / (2 * pi); t t(2:N); % 相位差只有N-1个点 if l 1 fn movmean(fn, l); % 滑动平均压制相位噪声 end相邻样本共轭相乘的相位差不会超过±π因此fn天然落在[-0.5, 0.5]区间这是归一化频率的物理边界换算物理频率时直接乘以fs。l参数决定频率曲线的平滑程度l取值效果适用场景1不平滑逐点抖动明显数据长、关注瞬态突变3~5轻度平滑保留主要事件一般振动、语音信号8~15强平滑噪声压制好但事件边缘变钝强噪声、只看整体趋势l取大之后事件起止时刻会被拉宽做故障定位时慎用。实际使用中按信号信噪比调整噪声明显时先从l5起步。4. 完整复现一遍构造非平稳信号、分解、算谱、验证重构误差4.1 构造测试信号调幅调频趋势噪声为了验证HHT的每一个环节构造一个瞬时频率随时间线性增加的调频分量叠加一个2Hz的低频趋势项和轻微白噪声fs 1000; % 采样率1000Hz T 0.5; t (0:T*fs-1) / fs; % 调幅系数0.6载波瞬时频率从50Hz线性增至65Hz x (1 0.6*cos(2*pi*2*t)) .* sin(2*pi*(50*t 15*t.^2)) ... 0.3*sin(2*pi*5*t) 0.05*randn(1, length(t)); figure; plot(t, x); xlabel(时间 (s)); title(测试信号AM 线性调频 低频趋势);线性调频的瞬时频率理论值是φ(t)对t求导的结果f(t)5030t0.5秒内从50Hz扫到65Hz。这个理论值后面直接当标尺对比HHT输出的瞬时频率是否落在正确区间。0.05的噪声幅度相对主分量很小不会破坏分解但足以测试频率曲线的平滑程度。4.2 跑通整条链路并核对三个关键结果[imf, A, f, tt] HHT(x, fs, 0); % 查看分解出几个分量 disp(size(imf)); % 逐列绘制IMF确认低频趋势是否被剥离到残差 figure; for k 1:size(imf, 2) subplot(size(imf, 2), 1, k); plot(t, imf(:, k)); ylabel(sprintf(IMF%d, k)); end xlabel(时间 (s));分解结果至少应该有3列第一列是50~65Hz的高频调频分量第二列是5Hz附近的趋势分量最后一列是残差。把IMF1展开看单个周期包络幅度应该随0.6·cos(2π·2t)缓慢摆动这就是调幅分量被正确分离的标志。重构验证是整条链路是否成立的最直接证据x_re sum(imf, 2); % IMF累加重构 rmse sqrt(mean((x(:) - x_re).^2)); fprintf(重构RMSE %.3e\n, rmse);RMSE在1e-10量级说明EMD分解完备原信号的信息没有丢失。这里也对应摘要里说的“通过将IMF分量累加重构得到平稳信号”——非平稳特性被剥离到各IMF的包络和频率调制里单个IMF在局部时间窗内可以当平稳窄带信号处理。再检查瞬时频率标尺k 1; % 第一个IMF idx tt 0.1 tt 0.4; % 避开端点效应区域 f_theory 50 30*tt(idx); f_est f(idx, k); fprintf(频率估计误差(中位数) %.4f Hz\n, median(abs(f_est - f_theory)));中位数误差应该在1Hz以内。如果误差偏大优先怀疑端点效应和噪声引起的相位抖动而不是去改代码。4.3 高频踩坑维度、返回值和残差混入错误现象常见原因排查与处理未定义函数或变量 emdEMD工具箱未装或未加入路径执行which emd为空则addpath后savepath矩阵维度不符合要求输入x是列向量或t长度与imf行数不一致统一x x(:).t用(0:N-1)/fs频率矩阵出现NaN信号含全0段解析信号相位无定义检查预处理避免零值段进入分解Hilbert谱高频区大片噪点残差被当成IMF参与变换绘图或分析时排除imf最后一列瞬时频率在端点剧烈摆动样条包络在端点缺少极值点约束只取中间段分析或用镜像扩展重构RMSE在1e-3量级输入含NaN/Inf或筛分未收敛用any(isnan(x))检查缩短数据长度特别提醒MATLAB自带emd与第三方工具箱的返回格式不同。内置版本要写[imf, residual] emd(x)再把残差拼回最后一列否则趋势项会丢失。判定当前用的是哪个版本看which emd返回路径里是否包含toolbox/signal字段。5. 验证HHT结果可信度的三个指标以及一个调参技巧5.1 正交性指数IMF之间有没有漏掉的交叉能量EMD理论要求分解完备且正交但有限次筛分会让IMF之间残留少量能量耦合。正交性指数IO统计所有两两IMF内积绝对值之和用信号能量归一化function io orthog_index(imf) % 正交性指数衡量IMF间的能量泄露越小越好 M size(imf, 2); s 0; for i 1:M for j i1:M s s abs(sum(imf(:, i) .* imf(:, j))); end end x_re sum(imf, 2); % 重构信号能量 io s / sum(x_re.^2); end调用io orthog_index(imf); fprintf(IO %.4f\n, io);工程经验上IO小于0.05是合格的。如果大于0.1说明两个相邻IMF在时域上仍有重叠振荡典型原因是模态混叠。此时先检查信号里有没有间歇性冲击——冲击会把高频事件叠加在低频背景上让第一个IMF失真这不是改代码能解决的问题要从信号预处理入手。5.2 边际谱与FFT频谱对比验证频率定位是否可信边际谱是把Hilbert谱沿时间方向累加得到的能量-频率分布。做法是把瞬时频率分箱累加对应幅值nbins 200; fmax fs / 2; fb floor(f(:) / fmax * nbins) 1; fb(fb 1) 1; % 防止噪声造成的负频率越界 fb(fb nbins) nbins; Hm accumarray(fb, A(:), [nbins 1]); fax ((1:nbins) - 0.5) * fmax / nbins; plot(fax, Hm); hold on; % 对比FFT幅值谱包络 X abs(fft(x)); X X(1:floor(end/2)1); X X / max(X) * max(Hm); f_fft (0:length(X)-1) * fs / length(x); plot(f_fft, X, --); xlim([0 200]); legend(HHT边际谱, FFT幅值谱);测试信号的边际谱峰值应该集中在50~65Hz频带和5Hz附近两个区域与FFT谱的峰位置一致。两者差异在于FFT把0.5秒内的能量平铺峰是一片边际谱沿时间累加能保留瞬时能量变化。如果边际谱出现FFT里完全没有的强峰说明某个IMF是伪分量回到5.1的正交性检查。5.3 端点效应对付办法与筛分容差的最后一个参数端点效应是样条包络在数据两端缺少极值点约束导致的包络摆动会污染IMF开头和结尾各几个振荡周期。常用的技巧是镜像扩展后分解再截回原长度nMirror 150; % 约为最低频分量周期的1~2倍 x_ext [fliplr(x(1:nMirror)), x, fliplr(x(end-nMirror1:end))]; imf_ext emd(x_ext); imf_cut imf_ext(nMirror1:end-nMirror, :); % 对比镜像处理前后第一个IMF中间段的相关系数 corr_val corr(imf(100:400, 1), imf_cut(100:400, 1)); fprintf(相关系数 %.4f\n, corr_val);相关系数在0.99以上说明分解稳定两端差异大是正常现象分析时用tt 0.05 tt 0.45这类区间截断即可。若用的是MATLAB自带emd还可以收紧筛分停止容差来减少低频IMF混入伪振荡imf emd(x, SiftRelativeTolerance, 0.01);第三方Rilling工具箱没有这个参数名需要在其全局选项里设置。容差收紧后IMF个数通常不变但每个IMF的包络会更光滑。最后一个实用技巧观察Hilbert谱时把频率轴改成对数刻度低频趋势分量和高频载波分量能同时看清set(gca, YScale, log)一行就够适合振动信号里故障特征频率和主频相差一个数量级的场景。本文还有配套的精品资源点击获取