做非平稳信号分析的朋友应该都有过这种体验拿一段实测的振动或生物电信号想看看里面到底有哪些频率成分FFT一出来谱上一堆峰看着挺全但哪个峰对应哪个时间段的什么变化完全说不清楚。原因很简单FFT把时间信息丢光了频率随时间变化的信号在频谱上只表现为一片模糊的宽带突起。我之前做旋转机械故障诊断时也卡在这个问题上后来把EMD分解、FFT和HHT这套组合流程跑通之后时变频率的轨迹才终于能看得一清二楚。这篇文章就把我实际跑通的这套组合方案完整拆开来讲。核心思路是EMD负责把复杂信号分解成有限个本征模态函数IMFFFT负责在频域里做全局“摸底”和分解前后的验证HHT负责给出时间和频率的联合分布。三者配合再加上MATLAB的可视化输出基本能解决绝大多数非平稳信号的分析需求。不管你是在读研究生处理实验数据还是工程师分析现场振动信号只要手里有一份MATLAB这套流程都可以直接抄作业。1. 方案选型为什么是EMDFFTHHT的组合1.1 一句话讲清三者的分工很多初学者会把EMD、FFT、HHT当成三个互相竞争的算法好像选了HHT就可以不用FFT。这是最大的误区。实际工程场景里三者是典型的“接力关系”各自承担不同环节的任务。FFT是全局频域分析工具擅长回答“信号里有哪些频率成分、幅度多大”缺点是一旦信号频率随时间变化频谱就会变得模糊因为FFT默认信号是平稳的。EMD是数据驱动的自适应分解工具不需要预设基函数能把信号按时间尺度从高频到低频逐层剥开得到一组IMF和一个残差这组IMF自带局部特征天然适合非平稳信号。HHT是在EMD基础上对每个IMF做Hilbert变换得到瞬时频率和瞬时幅值最终生成一张时间-频率-幅度的三维时频谱。一句话总结就是EMD负责把人打散FFT负责看静态构成HHT负责看动态演变。三者缺一不可EMD分解质量直接决定HHT谱图质量FFT则是验证每步结果是否合理的标尺。1.2 MATLAB环境准备与工具箱自查这套流程对MATLAB版本有一定要求我在代码里也踩过不少版本坑这里先帮你排雷。从R2018a开始MATLAB内置了emd函数属于Signal Processing Toolbox不需要额外安装工具箱。hht函数也是从R2018a开始内置的早期版本里这个功能分散在若干子函数里用起来比较复杂。如果你用的是R2018a以下的版本就需要下载第三方的EMD工具箱比如G. Rilling那套经典实现通过addpath手动加载。动手之前先自查一遍环境% 检查Signal Processing Toolbox是否存在 ver(signal) % 检查emd和hht是否可用 which emd which hht如果which返回空字符串说明当前版本没有内置函数。which emd能找到内置函数路径说明可以直接用这点很多教程不提示新手经常跑完报错才发现版本不对。我个人建议尽量升级到R2021b以后的版本新版的emd函数在效率和稳定性上有明显改进尤其是对长时间序列的支持好很多。1.3 一套可复用的分析流程总览我把这套流程固定成了一个标准管线每次拿到新数据都按这个顺序走数据清洗与预处理去直流、去趋势、检查采样率对原始信号做FFT获取全局频谱特征建立“频率先验”执行EMD分解得到IMF矩阵和残差逐个IMF做FFT验证各分量的频率范围是否合理对IMF执行Hilbert变换计算瞬时频率和瞬时幅值绘制HHT时频谱观察频率随时间的变化轨迹计算HHT边际谱与FFT频谱对照完成自检这套流程看起来简单但实际上每个环节都有细节要把握。后面我会用一个仿真信号和一个真实CSV数据走完整条管线把中间的关键参数和容易出错的地方全部摊开讲。2. 数据准备从仿真信号到真实CSV导入2.1 构造一个含非平稳成分的仿真信号为了把原理讲透我先构造一个贴近工程实际的仿真信号。这个信号模仿旋转机械升速过程的振动基频从20Hz逐渐升到80Hz同时叠加一个高频共振成分和工频干扰。fs 2000; % 采样率 2000 Hz T 2; % 信号时长 2 秒 t 0:1/fs:T-1/fs; % 时间向量 N length(t); % 采样点数 % 升速过程中的基频振动频率从20Hz线性增大到80Hz f_inst 20 30*t; % 瞬时频率 20Hz - 80Hz phi 2*pi*cumsum(f_inst)/fs; % 相位积分保证瞬时频率严格等于f_inst x1 1.5 * sin(phi); % 非平稳分量关键研究对象 % 高频共振成分356Hz固定频率 x2 0.8 * sin(2*pi*356*t); % 工频干扰50Hz x3 0.2 * sin(2*pi*50*t); % 构造最终信号并加入轻微噪声 x x1 x2 x3 0.03 * randn(size(t));这里有个非常关键的细节构造频率随时间变化的信号时不能直接写sin(2*pi*f_inst.*t)。如果直接乘t求解出的瞬时频率会变成f_inst t*f_inst多出一项实际频率就偏了。正确做法是先用cumsum对瞬时频率做相位积分再取sin这样瞬时频率才是严格等于f_inst的。这个坑我在早期写仿真时栽过当时HHT谱图里扫频轨迹总是比预设偏高排查了半天才发现是信号构造方式错了。2.2 CSV真实数据导入与预处理细节仿真信号只是为了原理验证工程上更多时候是导入采集器输出的CSV或TXT文件。CSV导入看起来简单但有几个细节直接影响后续分析。% 方式一readmatrixR2019a及以上推荐 data readmatrix(vibration_data.csv); % 方式二如果第一行是表头用readmatrix自动跳过 data readmatrix(vibration_data.csv, NumHeaderLines, 1); % 方式三老版本兼容方案 data csvread(vibration_data.csv);导入之后先看变量尺寸和数据类型别急着分析whos data常见的CSV有两种格式一种只有一列信号幅值没有时间列另一种有两列第一列是时间戳第二列是幅值。如果是第一种时间轴需要根据采样率重建fs 1024; % 采样率必须从采集设备那里拿到不能猜 t_data (0:length(data)-1) / fs; x_data data(:, 1);如果是第二种需要注意时间戳的起始点是否为零以及时间间隔是否均匀。很多采集卡导出的CSV时间戳是相对设备启动的毫秒数而且可能存在丢帧导致时间轴不均匀。碰到不均匀时间轴EMD函数会直接报错或者结果不可靠这时候不要骗自己老老实实先做插值重采样。% 如果时间轴略微不均匀用线性插值重采样到均匀网格 t_uniform 0:1/fs:data(end, 1); x_uniform interp1(data(:, 1), data(:, 2), t_uniform, linear); x_uniform x_uniform(~isnan(x_uniform));2.3 采样率、去趋势和去直流容易忽略的三个坑这三个坑几乎每个新手都会踩而且都是后知后觉的类型。第一是采样率。采样率过低高频成分会混叠到低频段HHT时频谱里出现假频率采样率过高数据量巨大EMD分解耗时成倍增加。选采样率时保证最高关心频率不超过奈奎斯特频率的一半也就是目标频率上限不要超过fs/4留足余量。比如关心500Hz以内的信号采样率至少设2000Hz。第二是去趋势。实测信号经常有基线漂移表现为整体趋势项。EMD会把趋势项当成一个低频IMF分解出来这本身没问题但趋势项会对Hilbert变换产生干扰让低频段的瞬时频率出现异常的大幅波动。% 去除线性趋势 x_detrend detrend(x_data, linear); % 去直流去均值 x_clean x_detrend - mean(x_detrend);第三是去直流。直流分量说白了就是信号的均值EMD分解第一个IMF时对极值点的分布影响巨大。实测数据如果均值不为零分解后IMF1可能会被直流偏移拉出不必要的畸变。我把去直流这一步放在所有处理之前已经成了肌肉记忆。3. EMD分解实战用法、参数与IMF解读3.1 emd函数的基本调用与返回结果现在进入核心环节。先用仿真信号跑一遍EMD分解。% 执行EMD分解 [imf, residual] emd(x, Display, 1);输出imf是一个矩阵每一行是一个IMF分量从高频排到低频residual是残差向量代表信号的趋势项或均值。矩阵的列数等于信号长度行数由信号复杂度决定。如果是在R2018a之前的老版本环境或者想用第三方工具箱调用方式略有区别一般是用imf emd(x)返回结果里最后一列是残差。为了兼容新版MATLAB统一用[imf, residual]这种语法在新版下最稳妥。Display, 1会在命令窗口打印分解过程中的筛选迭代次数。这个参数看起来不起眼实际是排查问题的利器。正常分解时每个IMF的筛选迭代次数大概在10到50次左右。如果你发现某个IMF迭代了200次以上才收敛说明这个IMF对应的频率成分不稳定或者端点效应很严重需要对信号做预处理。3.2 IMF可视化与物理意义判读分解完成后第一件事不是急着画HHT谱而是先把IMF曲线全部画出来看一眼。这一步花不了几秒钟但对判断分解质量至关重要。n_imf size(imf, 1); figure(Color, w, Position, [100 100 850 900]); tiledlayout(n_imf 2, 1, TileSpacing, compact); nexttile; plot(t, x, k); ylabel(原始); title(EMD分解结果原始信号 IMF 残差); for k 1:n_imf nexttile; plot(t, imf(k, :), b); ylabel([IMF, num2str(k)]); % 统一y轴范围避免视觉误导 ylim([-max(abs(imf(k, :)))-0.1, max(abs(imf(k, :)))0.1]); end nexttile; plot(t, residual, r); ylabel(残差); xlabel(时间/s);画完之后怎么判读以我构造的仿真信号为例理想情况下应该看到IMF1是356Hz的共振成分频率稳定IMF2是50Hz工频干扰IMF3是20到80Hz扫频分量波形能看到频率逐渐变密的趋势IMF4或更低频段是趋势项。如果分解结果和预期对不上就要回头检查信号预处理。我习惯在判读时做一个“IMF频率分层”的检查把各IMF曲线的过零点密度看一眼高频分量过零密度大低频分量过零密度小正常情况下相邻IMF之间应该有明显的频率层级差异。如果两个相邻IMF过零密度差不多说明存在模态混叠。3.3 端点效应和模态混叠的现场应对EMD最大的两个原生毛病就是端点效应和模态混叠。端点效应表现为信号两端发散曲线像被风吹起来一样剧烈摆动。原因很简单EMD的包络拟合依赖极值点而信号端点处只有一侧有极值点拟合出来的包络会在端点失稳。处理手法有几招一是修改插值方式。EDM函数支持Interpolation参数改成pchip分段三次Hermite插值后端点处的包络稳定性比默认的spline好一些。[imf, residual] emd(x, Interpolation, pchip);二是截取边界。如果端点效应只影响两端一小段可以在分析时把头尾各截掉5%到10%只分析中间稳定段。工程信号里我经常直接这么做损失一点点时间长度换来结果可靠。三是镜像延拓。用信号左右两端各一段数据做镜像延拓扩展极值点覆盖范围然后再分解。这个手法在第三方工具箱里有现成函数新版内置emd没有直接开放这个参数需要自己实现。模态混叠则是另一个顽固问题。当信号里有两个频率成分比较接近时EMD可能无法把它们准确分离到两个不同IMF里而是混在一个IMF中。经典应对方案是EEMD集成经验模态分解思路是给原始信号加入多组小幅白噪声分别做EMD再对结果求平均用噪声扰动破坏极值点分布的不规则性。% 用新版MATLAB自带的EEMD做对比如果环境支持 [imf_eemd, residual_eemd] eemd(x, NumEEMD, 100, NoiseAmplitude, 0.02);补充一句EEMD也不是银弹加噪幅值选太小没效果选太大又会污染信号本身。经验法则是噪声幅值设为信号标准差的0.1到0.3倍集成次数50到200次。我一般先跑100次如果分解结果不稳定再往上加。4. FFT在整套流程中的三个作用4.1 分解前用FFT摸清全局频谱EMD分解之前先用FFT给信号做个“全身CT”非常有必要。原因很简单FFT结果可以告诉你信号里大概有哪些主要频率成分这些成分在哪个频段幅度量级如何。有了这个先验信息你才能判断EMD分解出的IMF是否合理。% 对原始信号做FFT N length(x); f_axis_single (0:floor(N/2)-1) * fs / N; X fft(x); X_mag abs(X(1:floor(N/2))) / N; % 单边谱幅值校正直流分量不乘2 X_mag(2:end) X_mag(2:end) * 2; figure(Color, w); plot(f_axis_single, X_mag, LineWidth, 1.2); xlabel(频率/Hz); ylabel(幅值); title(原始信号FFT频谱); xlim([0 500]);按我构造的仿真信号这个频谱图上应该能看到三处明显峰值50Hz附近工频、大约80Hz附近的宽带突起扫频分量的能量堆叠、356Hz附近高频共振。需要注意的是扫频分量在FFT频谱上不是一条线而是一段宽带突起。因为信号频率随时间从20Hz扫到80HzFFT把所有时刻的频率能量累计在整段频率范围内无法告诉你频率变化过程。这正是后面HHT要补上的信息。4.2 分解后对IMF做FFT验证分离效果EMD分解完对每个IMF单独做FFT就能验证分离效果到底好不好。figure(Color, w, Position, [100 100 850 900]); tiledlayout(n_imf, 1, TileSpacing, compact); for k 1:n_imf nexttile; imf_fft abs(fft(imf(k, :))); imf_fft imf_fft(1:floor(N/2)) / N; imf_fft(2:end) imf_fft(2:end) * 2; plot(f_axis_single, imf_fft, LineWidth, 1.1); ylabel([IMF, num2str(k)]); xlim([0 500]); if k 1 title(各IMF分量的FFT频谱验证); end if k n_imf xlabel(频率/Hz); end end判读标准是理想情况下每个IMF的频谱应该是“干净”的单峰或窄带谱峰值频率清晰可辨。如果某个IMF的频谱出现多个分隔明显的峰值说明这个IMF还没有彻底分离干净如果相邻IMF的频谱严重重叠说明分解出了问题。在实机调试中这个验证步骤帮我省了大量回头路。有一次我处理齿轮箱振动信号EMD结果看起来有6个IMF但逐个FFT后发现IMF3和IMF4的频谱几乎一样都是500Hz附近。后来排查发现是齿轮啮合频率的两个边带实际应该归为一个IMF通过这步验证我及时调整了参数避免拿错误结果往HHT里喂。4.3 FFT频谱做幅值校正的小细节这一节算是个补充细节但恰恰是很多人画完频谱后觉得幅度对不上的原因。FFT直接计算结果除以N之后得到的是单频成分的实际幅值的一半因为能量分散在正负频率两边。所以要得到真实幅值单边谱需要除N后把非零频段的幅值乘以2。如果不做这个校正356Hz共振分量的幅值会显示成0.4而不是0.8更容易让你产生“信号变小了”的错觉。代码上注意索引不要越界X_mag abs(fft(x)) / N; X_mag_single X_mag(1:floor(N/2)); % 取单边 X_mag_single(2:end) X_mag_single(2:end) * 2; % 2:end不包含直流加到窗函数也是另一个常见优化点。直接fft相当于加了矩形窗频谱泄漏明显。如果想要更精准的幅度估计可以在FFT前乘一个汉宁窗win hann(N, periodic); X_win abs(fft(x .* win)) / sum(win) * 2; % 幅值恢复系数用窗求和注意用了窗函数之后幅值恢复公式不再是2/N而是2/sum(win)因为窗函数已经改变了信号能量分布。这个细节很多教程不提但我实测下来用对窗函数后356Hz处的幅值误差能从15%降到2%以内。5. HHT时频谱构建与可视化输出5.1 从IMF到瞬时频率Hilbert变换原理简述HHT的关键在于从每个IMF里提取“瞬时频率”。瞬时频率的定义是相位对时间的导数但直接对一个信号求瞬时频率是没有意义的因为工程信号往往是多分量叠加瞬时频率概念无法直接套用。EMD解决了这个问题通过把信号分解成单分量IMF每个IMF在任意时刻只有一个主导频率此时求瞬时频率才有物理意义。求瞬时频率的标准做法是对IMF做Hilbert变换构造解析信号z hilbert(imf(k, :)); % 解析信号 inst_amp abs(z); % 瞬时幅值 inst_phase unwrap(angle(z)); % 瞬时相位 inst_freq diff(inst_phase) / (2*pi) * fs; % 瞬时频率瞬时频率的计算也可以用instfreq函数直接完成。hht函数内部就是封装了这套流程但对参数的控制不如手动精细。这里的关键点是unwrap如果不做相位解卷绕相位跳跃会导致瞬时频率出现巨大的尖峰伪影这是新手经常遇到的问题。5.2 hht函数绘制时频谱的关键参数新版MATLAB直接用hht函数就能出时频谱但默认参数画出来的图经常不理想。我把它长期固定成以下调用方式figure(Color, w, Position, [100 100 900 500]); [hs, f_hh, t_hh] hht(imf, fs, ... FrequencyLimits, [0 500], ... FrequencyResolution, 0.5); imagesc(t_hh, f_hh, hs); axis xy; xlabel(时间/s); ylabel(频率/Hz); title(Hilbert-Huang时频谱); colorbar;解释几个关键参数。FrequencyLimits限制显示频率范围默认是0到奈奎斯特频率实际使用时经常只关心某一段频带限制范围能让目标区域的细节更清楚。FrequencyResolution是频率分辨率单位Hz数值越小频率轴越细但计算量也会增大。对于2000Hz采样率、2秒数据0.5Hz的频率分辨率已经足够细腻。还有一个常被忽略的细节当IMF数量很多且幅度差异大时时频谱上的高幅值成分会压得低幅值成分几乎看不见。解决方法是把颜色轴改成对数刻度set(gca, ColorScale, log);这一步效果非常明显。我处理振动信号时356Hz共振成分的幅值是扫频分量的好几倍线性颜色下低频段几乎一片深蓝什么都看不清切换到对数颜色后20到80Hz的扫频轨迹清晰浮出水面。还有一点时频谱末端经常出现颜色特别亮的竖条这是端点效应的典型表现说明IMF在信号末尾处瞬时频率失稳。不用惊慌要么截掉两端要么在解读时跳过边界区域。5.3 边际谱与FFT频谱的对比验证HHT时频谱是时间-频率-幅值的三维展示但很多时候我们还想看一个一维结果来跟FFT对照那就是边际谱。边际谱就是把时频谱沿时间轴积分得到“每个频率上总能量贡献”的分布。% 边际谱 对时频谱在时间维上求和 marginal_spectrum sum(hs, 2); figure(Color, w); plot(f_hh, marginal_spectrum, b, LineWidth, 1.2); xlabel(频率/Hz); ylabel(幅度); title(HHT边际谱与FFT频谱对比); xlim([0 500]); hold on; % 叠加FFT频谱做对照做归一化量纲一致才可对比 fft_norm X_mag * max(marginal_spectrum) / max(X_mag); plot(f_axis_single, fft_norm, r--, LineWidth, 1.1); legend(HHT边际谱, FFT频谱(归一化), Location, northeast);这个对比是整套流程里我认为最有价值的自检手段。如果HHT边际谱和FFT频谱的主峰位置能对上说明EMD分解、Hilbert变换、时频谱计算整条链路都是正常的如果两者差异非常大比如HHT里多出一个明显的频率峰而FFT里没有说明分解或变换过程有问题需要回头检查参数。这里要注意HHT边际谱和FFT频谱的量纲不一样直接叠画没有意义必须先做归一化。归一化之后两者的主峰位置和相对形状应该基本一致但幅值没法直接对应。6. 常见问题排查与经验速查6.1 五个高频问题的定位与解决我在实际使用和帮朋友调试过程中遇到过不少反复出现的问题。整理成了速查表方便你直接对照排查。问题现象可能原因解决方案未定义函数或变量emdMATLAB版本低于R2018a或未安装Signal Processing Toolbox升级版本或安装工具箱或下载Rilling第三方工具箱并addpathIMF两端曲线剧烈发散端点效应导致包络拟合失稳改用Interpolation,pchip分析时截掉两端5%-10%数据相邻IMF频谱重叠严重模态混叠频率成分过近或噪声干扰极值点分布改用EEMD/CEEMDAN适当增加集成次数检查信号信噪比HHT时频谱一片糊看不到清晰频率轨迹颜色动态范围不合适或频率分辨率设置太高/太低尝试set(gca,ColorScale,log)调整FrequencyLimits和FrequencyResolutionFFT峰值幅度远小于真实信号幅值频谱泄漏或没有做幅值校正加汉宁窗用2/N无窗或2/sum(win)加窗做幅值恢复6.2 几点能少走弯路的实操心得这套流程跑多了以后我沉淀了几个固定习惯每次处理新数据都先照做省了很多重复排错的时间。第一条先看数据再跑算法。拿到任何一段信号先plot看一眼原始波形确认没有明显异常值、断点、台阶再跑预处理。别嫌这一步啰嗦信号里有NaN或Inf时emd会直接报错或者静默输出错误结果到时候排查的功夫远大于先看一遍波形的时间。第二条EMD分解前务必先做FFT哪怕只是粗略看一眼频谱。这样你能预判信号里有几个主要分量、大致在什么频段分解完之后对照判断是否合理。没有先验信息去解读IMF和盲人摸象差不多。第三条HHT时频谱的显示范围宁可窄一点。把FrequencyLimits设得太宽目标频段会被压缩成几像素宽什么都看不清。先设宽范围跑一遍确定主要频率范围后再收窄到关心频段精细调参。第四条不要迷信IMFs越多越好。有时候EMD会继续把噪声拆成高频IMF这些分量幅值小、频率杂乱对分析主信号没有帮助。用MaxNumIMF参数或用相关系数筛选有效IMF只保留与原始信号相关性显著的分量。第五条我个人最深的体会是这套方案的价值不在于某个算法多么“高级”而在于EMD的时频定位能力和FFT的全局检验能力刚好互补。单纯用FFT我们只能看到有哪些频率但看不到变化的时刻单纯用HHT会忽略分量分离的质量验证。两者结合才形成了一条数据驱动的非平稳信号分析闭环。另外补充一点后续可做的扩展当EMD处理强噪声信号力不从心时可以尝试CEEMDAN当需要自动识别故障特征时可以结合能量熵、样本熵等特征提取方法如果数据规模很大还可以把IMF提取的特征向量喂给分类器做模式识别。这些都是在这套EMDFFTHHT基线上可以自然延伸的方向等你把今天这套流程跑熟了再往这几个方向深入会顺畅很多。