MATLAB实现EMD+FFT+HHT非平稳信号时频分析完整流程

做信号处理这些年我遇到最多的需求不是“怎么把频谱图画出来”而是“这段信号里到底藏着什么成分它随时间是怎么变的”。FFT能告诉你频率却说不清时间短时傅里叶变换能看时频但窗口一选就定了分辨率而EMD加HHT这套组合天生就是用来处理非平稳、非线性信号的。这篇内容我打算用MATLAB把一个完整流程跑通EMD分解、每个IMF的FFT频谱分析、HHT时频分析再加上后面那套可视化方案。适合正在做故障诊断、生物医学信号分析、振动噪声处理或者被非平稳信号折磨得头大的人参考看完可以直接拿自己的数据套用。1. 整体设计与方案选型1.1 为什么偏偏是这个组合先说清楚三者的分工。EMD经验模态分解负责把复杂信号拆成若干个本征模态函数IMF拆出来的每个IMF都是单分量窄带信号FFT负责告诉我们每个IMF里最主要的频率成分在哪HHT希尔伯特-黄变换则是在EMD基础上做希尔伯特变换求出每个IMF的瞬时频率和瞬时幅值最终画出一张“时间-频率-幅值”三维图谱。为什么不用小波分析小波需要选基函数选db4还是morlet结果差异有时候大到让人怀疑人生。EMD最大的特点是自适应——它不预设基函数完全根据信号本身的时间尺度特征逐层分解。这一点在处理实际采集数据时非常占便宜因为现场信号往往是多成分叠加、频率漂移、甚至有间歇性冲击的你根本不知道什么基函数能完美匹配它。FFT在这个流程里更像一把尺子。EMD分解完之后到底准不准、每个IMF对应什么物理含义光看时域波形很难判断把IMF拉去做FFT主频一测便知。它的作用不是替代HHT而是给后续HHT分析提供“频域验证”。HHT是这一整套的压轴。它给出的是瞬时频率不是平均频率。对变频信号、脉冲信号这类非平稳成分只有瞬时频率才能如实反映频率随时间的变化。最终画的Hilbert谱横轴是时间、纵轴是频率、颜色代表能量强度一张图把信号的“家底”全抖出来。1.2 这套方案适合处理什么类型的数据虽然这套流程看起来很“万金油”但我必须把适用边界说清楚不然你拿去硬套会翻车。适合的典型场景有几个。第一是旋转机械故障诊断比如轴承磨损、齿轮断齿振动信号里往往叠加了转频、啮合频率以及冲击引起的宽带响应EMD能把不同频带拆开HHT能定位冲击发生的时间。第二是生物医学信号像心电ECG、脑电EEG这类信号频率低、非平稳、还夹杂大量噪声和基线漂移EMD可以把呼吸干扰和工频干扰先剥离开。第三是地球物理和声学信号比如地震波、水下声信号这类信号常常是短时瞬态的FFT做出来只是一团模糊的频谱只有HHT能还原瞬态过程。不适合的场景也要注意。如果信号的信噪比太低噪声完全淹没了有效成分EMD会先把噪声拆出来导致分解结果不稳定如果采样率过低频率混叠会直接污染后续所有分析。另外EMD分解对数据长度有一定要求太短的信号段分解不出可靠IMF我个人经验是至少保留5到10个完整振荡周期以上。1.3 可视化在流程里到底起多大作用很多人把可视化当成“最后画个图交差”这是对信号处理最大的误解。在实际工程里可视化的意义是验证和诊断。你在命令窗口里看IMF的数值看不出模态混叠但你把IMF逐条画在同一个图里一眼就能发现相邻两条IMF的波形长得几乎一样——这就是过分解的典型迹象。你打印频域峰值可能被数值误差误导但把FFT幅频谱画出来谱峰是否尖锐、是否有边频带一目了然。更关键的是Hilbert谱如果不经过良好可视化设计再好的分析结果也表达不出来。颜色映射选错了整个能量分布会糊成一片频率轴范围没限制高频噪声会让有效信号被压缩成一条细线时间轴分辨率不够瞬态事件会被平滑掉。所以我在本文里把可视化当做一个正式的环节来讲提供可直接用的绘图方案。2. 核心原理与关键细节2.1 EMD是怎么把信号拆开的EMD的基本假设是任何复杂信号都可以看作若干个“本征模态函数”与一个趋势项残差的叠加。IMF需要满足两个条件一是极值点个数和过零点个数相等或最多相差一个二是上下包络的均值为0即关于时间轴局部对称。分解过程像一个反复剥离的过程。先找出信号的所有局部极大值和极小值用三次样条曲线分别拟合出上包络和下包络取上下包络的平均值得到一个均值曲线用原始信号减去这个均值曲线得到一个“初步分量”判断它是否满足IMF条件。如果不满足就用这个分量重新执行一遍上述操作直到满足停止条件得到第一个IMF。然后用原始信号减去这个IMF对剩余部分继续筛分直到残差是单调函数或幅度很小为止。这个过程看起来简单实际使用中有两个关键参数值得留意。一是筛分停止准则MATLAB的emd函数内部默认有一套迭代停止条件一般不需要动。二是边界效应因为信号端点处无法同时得到左右两侧的极值信息包络拟合经常在端点处发散分解结果的前后两端会明显失真。我在实际处理中一般会先把信号两端做数据延拓预测再把IMF前后各截掉一部分只保留有效区间。这不是可有可无的操作边界效应处理不当会直接影响瞬时频率的计算结果。2.2 FFT验证频域成分的标尺FFT是数字信号处理的基石但用它对IMF做频谱分析时有几个细节必须注意。一是频率分辨率。FFT的频率分辨率是Fs/N其中Fs是采样率N是参与FFT的数据点数。要想让两个频率接近的成分在频谱上分开必须靠增加数据长度而不是提高采样率采样率决定了分析频率范围数据长度决定了频率分辨能力。二是频谱泄漏。如果信号频率不是FFT频率格点的整数倍能量会泄漏到相邻频点表现为谱峰变宽、旁边出现伪峰。解决办法是加窗函数常用的有汉宁窗、汉明窗、布莱克曼窗。加窗也会带来主瓣变宽的代价所以频谱分辨率不是越高越好要根据实际频率间隔取舍。三是幅值归一化。用FFT算出的结果是对称的复数一般不直接看。要画单边幅频谱需要取FFT结果的前一半并对幅值乘以2/N。直流分量则只除以N不需要乘2。很多新手画出来的幅值跟时域信号对不上往往就是归一化没做对。2.3 HHT与Hilbert谱的核心概念HHT的核心是对每个IMF做希尔伯特变换构造解析信号。解析信号有实部和虚部幅值就是瞬时幅值相位的导数就是瞬时频率。这里有个容易混淆的地方瞬时频率的定义不是对信号直接求导而是对解析信号的相位求导。对单分量、窄带信号瞬时频率非常稳定且能精确刻画频率随时间的变化但如果信号本身是多分量混叠的瞬时频率会变成无意义的“摆动”这也是为什么HHT要求先把信号做EMD分解成IMF——只有IMF才能保证希尔伯特变换得出有物理意义的瞬时频率。希尔伯特谱的构造方式是把每个时间点上各IMF的瞬时频率和瞬时幅值记录下来放到“时间-频率”平面上用颜色映射表示幅值大小得到H(t,f)这就是Hilbert谱。把H(t,f)对时间求积分或求和得到的是边际谱表示每个频率成分在整个时间跨度内累计的能量强度。边际谱和FFT频谱看着像但物理含义不同。FFT频谱是正弦分量的稳态叠加边际谱是瞬时频率的累计能量分布。对非平稳信号边际谱更能反映出真实存在的频率成分同时谱峰也更尖锐。3. MATLAB实操从仿真信号到真实数据3.1 环境准备与工具箱检查动手之前先确认环境。本方案需要MATLAB基础环境版本建议R2021a及以上因为自带的emd函数在这个版本之后已经比较稳定语法也很清爽。如果你还在用老版本也可以装G. Riling的EMD工具箱代码调用方式和系统自带的不同但核心思路一样。需要安装的工具箱主要是Signal Processing ToolboxFFT和窗函数都在里面hht函数也依赖它。另外建议装上Statistics and Machine Learning Toolbox后面绘制边际谱做统计时会用到但不是必需。环境确认就两步% 检查emd函数是否存在 if exist(emd, file) disp(emd函数可用); else disp(当前版本没有内置emd请考虑升级MATLAB或下载第三方工具箱); end % 检查hht函数 if exist(hht, file) disp(hht函数可用); else disp(当前版本没有内置hht请确认Signal Processing Toolbox已安装); end实测下来大部分人的问题不是工具箱没有而是工作路径或变量名跟内置函数冲突。尤其注意不要把自己的脚本命名为emd.m或fft.m否则MATLAB会优先找你的脚本导致调用内置函数时直接报错或者递归调用崩溃。3.2 构造一组带“故事”的仿真信号学习这套流程最忌讳一上来就拿真实数据练手因为真实数据里干扰太多你不知道问题出在算法上还是数据上。先用一组仿真信号把流程跑通每个成分都有明确的物理意义对号入座再上真实数据就容易了。我构造的仿真信号包含四个成分直流或低频趋势项模拟传感器漂移或环境缓慢变化一个恒频正弦模拟稳定的周期振动源一个频率随时间线性升高的调频信号模拟转速升高的设备振动一个高频脉冲衰减信号模拟间歇性冲击比如齿轮敲击或轴承局部缺陷。采样率设置为1024 Hz采样时长2秒这样一个周期信号能覆盖40多个周期足够EMD分解和FFT分辨。Fs 1024; % 采样率 1024 Hz T 2; % 采样时长 2 秒 t (0:1/Fs:T-1/Fs); % 时间序列 N length(t); % 成分1低频趋势项用低频正弦模拟 trend 0.8 * sin(2*pi*2*t); % 2 Hz 低频分量 % 成分2恒定频率正弦波 c1 1.5 * sin(2*pi*50*t); % 50 Hz 稳态振动 % 成分3频率随时间变化的调频信号 c2 sin(2*pi*(20*t 8*t.^2)); % 频率从20Hz线性增加到52Hz左右 % 成分4间歇性冲击衰减 impulse zeros(size(t)); impulse_time [0.4, 0.9, 1.4, 1.8]; for k 1:length(impulse_time) idx find(abs(t - impulse_time(k)) 0.05); % 冲击持续0.1秒 tt t(idx) - impulse_time(k); impulse(idx) impulse(idx) 2.0 * sin(2*pi*200*tt) .* exp(-40*abs(tt)); end % 叠加得到最终信号 x trend c1 c2 impulse;这几个成分各有特点2 Hz趋势项会让信号整体上下漂移50 Hz稳态正弦是“背景噪音”调频信号频率一直在变最能体现HHT的价值200 Hz冲击衰减则是典型的非平稳瞬态信号。混合起来后直接用FFT只能看到几个模糊的大峰几乎无法分辨频率随时间的变化这正是引入EMD和HHT的动机。先画一下原始信号的时域波形对整体形态有个感觉figure(Color,w,Position,[100 100 1200 400]); plot(t, x, k); xlabel(时间 (s)); ylabel(幅值); title(原始混合信号时域波形); xlim([0 T]); grid on;3.3 EMD分解与IMF时域可视化调用emd函数得到IMF矩阵和残差代码非常简单imf emd(x);imf的每一列是一个IMF列数由信号特性决定。最后一列通常是残差趋势项。我的经验是在正式分解前先对信号做一次去趋势处理例如用detrend去掉线性趋势可以显著减少EMD把趋势成分拆散的概率让前几个IMF更集中。分解完成后把每个IMF和原始信号画在一起用堆叠图的形式看[numIMF, ~] size(imf); % 注意MATLAB的emd返回的imf维度可能是IMF数量 x N具体以运行结果为准 figure(Color,w,Position,[100 100 900 800]); for k 1:size(imf, 1) subplot(size(imf, 1)1, 1, k); plot(t, imf(k, :), b); ylabel([IMF, num2str(k)]); if k 1 title(EMD分解结果); end xlim([0 T]); grid on; end subplot(size(imf, 1)1, 1, size(imf, 1)1); plot(t, x, k); ylabel(原始信号); xlabel(时间 (s)); xlim([0 T]); grid on;观察这张图的方法不是看它画得美不美而是看几个关键点。第一个IMF应该是信号里频率最高的成分通常对应冲击或高频载波中间的IMF对应中等频率最后的残差是趋势项。如果发现相邻IMF在某一时间段的波形高度相似说明存在模态混叠后面要处理。实际运行中我的经验是第二个IMF往往会捕获50 Hz稳态正弦第三个IMF会捕获调频信号但EMD对调频成分的分解并不总是完美的它可能会把调频信号的频率范围切分成两段分给两个IMF。这种“过分解”现象在信号频率跨度较大时经常出现不能盲目认为每个IMF都有独立物理意义要结合FFT和Hilbert谱验证。3.4 对每个IMF做FFT频谱分析分解完并不代表结束先用FFT给每个IMF做一次“体检”。我习惯写一个自定义函数集中处理去均值、加窗、FFT、幅值修正和频率轴计算。function [f, amp] myfft(x, Fs, winType) x x(:); N length(x); x x - mean(x); % 去直流 if nargin 3 || isempty(winType) winType hann; end win window(winType, N); xw x .* win; % 加窗 X fft(xw, N); P2 abs(X / N); P1 P2(1:floor(N/2)1); P1(2:end-1) 2 * P1(2:end-1); f (0:floor(N/2)) * Fs / N; amp P1; end对前几个非趋势IMF调用这个函数把频谱画在一个图里figure(Color,w,Position,[200 100 1200 700]); nShow min(4, size(imf, 1)); for k 1:nShow [f_cur, amp_cur] myfft(imf(k, :), Fs, hann); subplot(nShow, 1, k); plot(f_cur, amp_cur, b); xlabel(频率 (Hz)); ylabel(幅值); title([IMF, num2str(k), FFT幅频谱]); xlim([0 300]); grid on; end频谱结果通常会给你几个明确信号。第一个IMF的主峰如果在200 Hz附近说明它确实捕获到了高频冲击第二个IMF主峰在50 Hz对应稳态正弦如果第三个IMF的峰值在30到50 Hz之间有一串连续谱峰说明调频成分被它部分吸收。你还可以记录每个IMF的最大峰频率做出一张表格辅助判断IMF序号主峰频率(Hz)可能对应成分IMF1约200冲击衰减IMF2约50稳态正弦IMF3流动变化调频信号残差低频窄带趋势项这里有个常见误区不要把FFT的幅值直接理解为IMF的能量大小。加窗后幅值会比实际值偏低不同窗函数的修正系数还不一样。对比几个IMF时只看相对大小和主峰位置就够真要算能量用bandpower函数更准确。3.5 绘制Hilbert谱与边际谱做完FFT验证接下来进入重头戏——HHT时频分析。MATLAB自带的hht函数可以一步完成希尔伯特变换和谱图计算。% 只对前几个IMF做HHT残差趋势项不需要参与时频分析 imf_used imf(1:nShow, :); % 调用hht函数 [hspec, f_hht, t_hht] hht(imf_used, Fs);hspec就是二维希尔伯特谱矩阵f_hht是频率轴t_hht是时间轴。直接用hht还可以自动画图但自动画图的配色和轴范围经常不适合发表和汇报我建议手动控制绘图参数。手动绘制Hilbert谱核心是用imagesc或surf把二维矩阵显示成伪彩图figure(Color,w,Position,[300 100 1200 500]); imagesc(t_hht, f_hht, hspec); axis xy; xlabel(时间 (s)); ylabel(频率 (Hz)); title(希尔伯特谱 H(t,f)); xlim([0 T]); ylim([0 250]); colorbar; colormap(jet); shading interp;这张图就是整段信号分析的“总览图”。你会在时间轴上看到几条随频率变化的亮线50 Hz处一条水平亮线代表稳态正弦一条从20 Hz斜着爬到50 Hz的亮线代表调频信号在0.4、0.9、1.4、1.8秒位置出现竖条高频亮块代表冲击事件。趋势项在2 Hz附近因为频率轴下限我设到0它会被压缩在最底部这正好验证了分解的正确性。边际谱的求法是对hspec沿时间轴求和因为时间轴采样是等间隔的直接用sum即可marginal_spectrum sum(hspec, 2); figure(Color,w,Position,[300 100 1200 400]); plot(f_hht, marginal_spectrum / max(marginal_spectrum), b); xlabel(频率 (Hz)); ylabel(归一化边际谱幅值); title(边际谱归一化); xlim([0 300]); grid on;边际谱和FFT频谱做对比你会发现边际谱的谱峰更集中各频率成分的界限更分明。原因是FFT对整段信号做平均调频信号的能量被平均到一个频带范围内形成比较宽的峰而边际谱把每个时刻的瞬时频率累积起来能量分布更贴近真实频率变化。画Hilbert谱时有一点要特别提醒hht函数的默认频率分辨率较高如果你的信号比较长二维矩阵会非常大绘图会变慢。可以在调用时通过频率点数和时间点数的设置控制输出规模具体参数请参考MATLAB文档不同版本接口略有差异。3.6 导入CSV数据走完整流程真实工程里信号往往存在CSV文件里第一列是时间第二列或之后几列是传感器数据。用MATLAB导入并走完整流程我会写成下面这样的模板。% 读取CSV % 假设CSV有表头前两列分别是 time 和 signal data readmatrix(sensor_data.csv); % 提取时间和信号列 t_raw data(:, 1); x_raw data(:, 2); % 估计采样率 Fs 1 / mean(diff(t_raw)); % 注意如果时间间隔不均匀建议先重采样到等间隔 % 如果时间列是非数值型或CSV存在空值先做清洗 x_raw fillmissing(x_raw, linear); % 去除异常跳变点和线性趋势 x_detrend detrend(x_raw); % 执行EMD imf emd(x_detrend);导入CSV这一步常见坑很多。第一个坑是时间列单位不统一有些设备导出的时间戳是Unix毫秒有些是相对时间秒直接diff得到的采样周期可能差了好几个数量级导致频率轴完全错误。第二个坑是数据里可能有NaN或Infemd函数遇到NaN会报错所以先fillmissing做插值。第三个坑是采样率不均匀比如某些采集卡在系统繁忙时会丢点时间间隔忽大忽小这会导致FFT分析结果出现虚假频率成分至少要先画一下时间间隔分布确认均匀后再用。我通常处理真实数据的流程是先用findpeaks检查波形质量剔除明显坏道再做去趋势和带通滤波最后才进入EMD。预处理的好坏决定了分解的质量花在预处理上的时间永远值得。4. 常见问题与排查技巧实录4.1 模态混叠模态混叠的表现是某个IMF的波形时密时疏FFT频谱出现几个不相关的大峰相邻IMF的波形相似度极高。成因通常是一个间歇性的高频信号混入低频振荡成分里导致局部极值点分布不均包络均值失去意义。举个最典型的例子一个50 Hz正弦波里叠加了一个只在0.2秒内出现的高频冲击EMD在冲击出现的区间会把高频成分拆到IMF1但没有冲击的区间IMF1又被低频成分主导结果IMF1在时间上出现“频率断裂”频谱上出现宽频乱峰。解决办法有几个方向。一是改用集合经验模态分解EEMD或完全自适应噪声集合经验模态分解CEEMDAN通过在信号上添加有限次数的白噪声让不同尺度的成分被平均化地分离出来再对多次分解结果求平均可以有效抑制模态混叠。缺点是需要多次迭代计算量明显上升而且噪声幅度的初始值需要调参。二是先对信号做带通滤波把频率明显不在同一量级的成分分开再分别做EMD。比如知道50 Hz和200 Hz是无关成分就先做个100 Hz的高通滤波器把高频冲击先滤出来剩下的低频段再做EMD。三是检查采样率是否足够。如果冲击信号很窄但采样率只比奈奎斯特频率高出一点点EMD在极值点插值时无法正确还原冲击波形也会诱发模态混叠。提高采样率或降低分析频带往往立刻缓解。4.2 边界效应在信号的两端没有充分的极值点参与包络拟合三次样条包络会在端点处剧烈弯曲导致IMF端点发散。这在Hilbert谱上表现为时间轴两端出现不合理的超高或超低频“伪影”像两只触角一样向两端伸展。最彻底的办法是数据延拓。常见延拓方法包括镜像延拓、对称延拓、多项式拟合延拓、自回归预测延拓。镜像延拓使用最多它在端点处把信号“反射”出去一段让包络拟合在端点附近有数据支撑之后再把延拓部分截掉。MATLAB的emd内置了边界处理策略但实际表现并不总是理想。我的习惯是如果信号非常长直接放弃两端各5%到10%的数据段。比如一段10秒的信号我只对中间8秒做分析和可视化两端的瞬时频率抖动就不去管它。这样做损失的信息不大但能把谱图的“触角”彻底消除。如果信号本身就短比如只有1秒那就得用镜像延拓这种更精细的方法延拓长度取信号长度的10%到20%延拓太多反而会引入虚假振荡。4.3 FFT频谱读图容易犯的低级错误频率轴算错是我见过最多的错误。f (0:N-1)*Fs/N这个公式看似简单但如果你用了单边谱就要用N/2而不是N-1如果你用linspace(0, Fs, N)生成频率轴又和FFT输出长度对不上。画图之前先在时域上人工确认一下波形里大概有多少个周期的振荡再对照谱峰位置看是否一致这是最直接的验证方法。还有很多人把FFT结果的第一个点也就是0 Hz直流直接当噪声忽略但它在EMD流程里代表趋势项的能量如果你关心信号的整体漂移就要看这个点。另外加窗之后频谱幅值会被压低不同窗函数的幅值修正系数不同比较不同IMF的频谱幅值时最好统一窗函数不然结论没有可比性。4.4 Hilbert谱看不清楚怎么办hht画出来的谱图如果一片模糊通常是三个原因。第一频率轴范围太大。比如采样率是10 kHz信号主要成分在500 Hz以内你不限制ylim整个谱图的有效信息被压缩在底部十分之一的高度里当然看不清楚。处理方法就是将频率轴限制在有效频带内。第二颜色映射不合理。jet色图高亮区域多容易让人觉得大片区域都有强能量换成parula或turbo色图后低幅值部分显示更平滑。如果要突出弱信号可以手动设置caxis的范围不要总是让色标自动缩放。第三瞬时频率在端点或不连续处会跳变画出来是杂乱无章的噪声点。这时候先看IMF的质量如果某个IMF本身就不光滑它的瞬时频率就没有意义应该剔除或重新分解。还要说明一点hht对输入IMF的数量有要求只输入前几个有效的IMF就够了把残差也丢进去会让谱图底部多一条超低频亮带掩盖其他成分。4.5 IMF太多不知道看哪个EMD分解经常产生七八个甚至十几个IMF但并不是每个都有物理意义。判断哪个IMF值得分析我用三个标准。第一是相关系数。计算每个IMF与原始信号的皮尔逊相关系数相关系数高的几个IMF通常是主要能量所在相关性很低的通常包含大量噪声和边界伪影。第二是频谱主峰清晰度。做主峰明显的IMF保留主峰宽平、多个杂峰电平接近的IMF往往是噪声或过度分解的产物。第三是IMF的振荡规律性。如果一个IMF的时域波形频率忽高忽低、毫无章法它的瞬时频率大概率没有物理意义。实操中我常常只保留前3到5个IMF做后续分析。在故障诊断里有效成分通常就集中在头两个IMF剩下的大多是噪声在心电分析里基线漂移会被拆到最后的残差里中间某个IMF对应QRS波群的频带。数据结构不同挑选策略也要跟着变。最后说一个个人小习惯每次跑完这套流程我会把图保存成带时间戳的PNG并在图里直接标注采样率、数据长度、分析区间、窗函数这些关键参数。三个月后再回头看你永远不会记得当年这个数据是用什么窗处理的但图上的标注会替你想起来。这个习惯救过我很多次也建议你从今天开始用上。