EMD+FFT+HHT组合实战：用MATLAB实现非平稳信号时频分析

做非平稳信号分析的朋友应该都有过这种体验：拿一段实测的振动或生物电信号，想看看里面到底有哪些频率成分，FFT一出来，谱上一堆峰，看着挺全，但哪个峰对应哪个时间段的什么变化，完全说不清楚。原因很简单，FFT把时间信息丢光了，频率随时间变化的信号在频谱上只表现为一片模糊的宽带突起。我之前做旋转机械故障诊断时也卡在这个问题上，后来把EMD分解、FFT和HHT这套组合流程跑通之后，时变频率的轨迹才终于能看得一清二楚。

这篇文章就把我实际跑通的这套组合方案完整拆开来讲。核心思路是：EMD负责把复杂信号分解成有限个本征模态函数（IMF），FFT负责在频域里做全局“摸底”和分解前后的验证，HHT负责给出时间和频率的联合分布。三者配合再加上MATLAB的可视化输出，基本能解决绝大多数非平稳信号的分析需求。不管你是在读研究生处理实验数据，还是工程师分析现场振动信号，只要手里有一份MATLAB，这套流程都可以直接抄作业。

1. 方案选型：为什么是EMD+FFT+HHT的组合

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 一套可复用的分析流程总览

我把这套流程固定成了一个标准管线，每次拿到新数据都按这个顺序走：

1. 数据清洗与预处理（去直流、去趋势、检查采样率）
2. 对原始信号做FFT，获取全局频谱特征，建立“频率先验”
3. 执行EMD分解，得到IMF矩阵和残差
4. 逐个IMF做FFT，验证各分量的频率范围是否合理
5. 对IMF执行Hilbert变换，计算瞬时频率和瞬时幅值
6. 绘制HHT时频谱，观察频率随时间的变化轨迹
7. 计算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导入看起来简单，但有几个细节直接影响后续分析。

% 方式一：readmatrix（R2019a及以上推荐） 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扫到80Hz，FFT把所有时刻的频率能量累计在整段频率范围内，无法告诉你频率变化过程。这正是后面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 五个高频问题的定位与解决

我在实际使用和帮朋友调试过程中，遇到过不少反复出现的问题。整理成了速查表，方便你直接对照排查。

6.2 几点能少走弯路的实操心得

这套流程跑多了以后，我沉淀了几个固定习惯，每次处理新数据都先照做，省了很多重复排错的时间。

第一条，先看数据再跑算法。拿到任何一段信号，先plot看一眼原始波形，确认没有明显异常值、断点、台阶，再跑预处理。别嫌这一步啰嗦，信号里有NaN或Inf时，emd会直接报错或者静默输出错误结果，到时候排查的功夫远大于先看一遍波形的时间。

第二条，EMD分解前务必先做FFT，哪怕只是粗略看一眼频谱。这样你能预判信号里有几个主要分量、大致在什么频段，分解完之后对照判断是否合理。没有先验信息去解读IMF，和盲人摸象差不多。

第三条，HHT时频谱的显示范围宁可窄一点。把FrequencyLimits设得太宽，目标频段会被压缩成几像素宽，什么都看不清。先设宽范围跑一遍，确定主要频率范围后，再收窄到关心频段精细调参。

第四条，不要迷信IMFs越多越好。有时候EMD会继续把噪声拆成高频IMF，这些分量幅值小、频率杂乱，对分析主信号没有帮助。用'MaxNumIMF'参数或用相关系数筛选有效IMF，只保留与原始信号相关性显著的分量。

第五条，我个人最深的体会是：这套方案的价值不在于某个算法多么“高级”，而在于EMD的时频定位能力和FFT的全局检验能力刚好互补。单纯用FFT，我们只能看到有哪些频率，但看不到变化的时刻；单纯用HHT，会忽略分量分离的质量验证。两者结合，才形成了一条数据驱动的非平稳信号分析闭环。

另外补充一点后续可做的扩展：当EMD处理强噪声信号力不从心时，可以尝试CEEMDAN；当需要自动识别故障特征时，可以结合能量熵、样本熵等特征提取方法；如果数据规模很大，还可以把IMF提取的特征向量喂给分类器做模式识别。这些都是在这套EMD+FFT+HHT基线上可以自然延伸的方向，等你把今天这套流程跑熟了，再往这几个方向深入会顺畅很多。