EMD信号去噪实战：MATLAB实现与IMF筛选策略

简介：面向需要在MATLAB中对一维信号进行去噪的开发者与研究人员，这里提供基于经验模态分解（EMD）的完整示例代码。资源压缩包共2个文件、均为m脚本，体积仅6KB，包含一个核心去噪函数和一个可直接运行的演示脚本，便于从零开始运行与调试。已有335人下载学习，适合信号处理初学者作为入门参考。代码覆盖EMD去噪的主要环节：定位信号的局部极值并构造上下包络，计算平均包络后从原始信号中分离出内在模态函数（IMF），反复迭代至残余分量不再满足IMF条件；随后对各IMF应用希尔伯特变换获得瞬时频率与幅度，并将有效分量重构为去噪后的信号。演示脚本还展示原始信号、各IMF、去噪结果及瞬时频谱的可视化，帮助读者理解EMD如何在非线性、非平稳信号中保留主要特征并去除噪声，从而轻松迁移到自己的数据中。

1. 为什么一维信号去噪该先试EMD，而不是直接上滤波器

做过振动分析或生理信号处理的人都有过这种体验：信号里既有缓慢的趋势，又有突然的冲击，用带通滤波器一滤，冲击被抹平了，趋势也被改动了。原因是传统滤波器依赖事先设定的通带频率，而真实一维信号往往是非线性、非平稳的，瞬时频率本身就在变。经验模态分解（EMD）的做法完全不同：它不预设基函数，直接根据信号的局部极值把信号自适应地拆成若干本征模态函数（IMF），然后通过选择性重构实现去噪。这套做法对冲击特征、趋势漂移和随机噪声混杂的信号特别有效。如果你手里只有一段一维数据，想在MATLAB里快速看到分解效果，同时理解每一步在做什么，那么基于EMDdenoise.m和example.m这套代码往下拆是最直接的路径。适合信号处理入门者，也适合想用EMD做冷启动分析的工程师。

2. EMD分解的数学直觉与MATLAB内置实现

2.1 什么是IMF，为什么它比“频率带”更适合去噪

IMF需要满足两个条件：整个信号上极值点数目与过零点数目相等或最多差一；在任意局部，由局部极大值定义的上包络和局部极小值定义的下包络的均值为零。第一条保证了IMF是一段窄带振荡，第二条保证了瞬时频率有物理意义。去噪场景里，噪声往往对应振荡密集的高频IMF，而真实趋势和低速变化落在后面的低频IMF和残余项里。与传统滤波器按固定频带切分不同，IMF的划分是数据驱动的，因此当信号频率随时间漂移时，EMD依然能把不同时刻的振荡归入合适的IMF，不会出现滤波后频率成分被“切歪”的问题。

2.2 筛分（sifting）迭代：包络均值与停止条件

EMD的核心是筛分过程。对输入信号x，先找到所有局部极大值和局部极小值，用三次样条分别拟合上包络和下包络，求平均包络m，然后让h = x - m。如果h满足IMF条件就作为本征模态，否则对h重复这个过程。每提取一个IMF后，用r = r - IMF得到残余，再对残余继续筛分，直到残余单调或极点数不足两个。这个迭代必须有一个停止条件，否则会无限筛下去。最常用的判据是标准偏差（Sifting Tolerance）：

SD = sum((h_{k-1} - h_k)^2) / sum(h_{k-1}^2)

当SD小于阈值时，认为筛分收敛。典型阈值在0.2到0.3之间。MATLAB内置的emd函数也提供类似参数，常见设置如下：

2.3 在MATLAB里用内置emd函数快速看分解结果

如果你只是想先看看EMD对一段信号做了什么事，不需要立刻进入手写实现。MATLAB从R2018a开始把emd作为主工具箱函数提供，调用方式很直接：

x = load('sensor_signal.mat').x; % 读取一维信号 imf = emd(x); % 默认参数分解

这段代码里，x必须是列向量。如果x是行向量，转置一下再用。emd返回一个二维矩阵，行对应IMF序号，第一行通常是高频分量，最后一行是残余趋势。你可以直接plot(imf')查看各个IMF的形态。内置函数的好处是经过优化，稳定性好，但问题在于它对去噪场景没有做“哪些IMF保留”的判断，只会把分解结果全给你。真正要去噪，还需要自己写筛选逻辑。

提示：不同MATLAB版本对内置emd的参数名和默认值有差异，实际使用前用doc emd查看当前版本文档，别拿旧命令硬套。

3. EMDdenoise.m手写实现：从筛分到重构的完整拆解

3.1 输入输出与预处理

手写EMD的意义在于理解边界条件和参数控制。EMDdenoise.m的输入通常是一个含噪一维信号，输出是去噪信号和分解出的IMF矩阵。函数开头先做预处理：保证输入是列向量，检查长度是否足够，处理NaN，然后确定筛分阈值和最大迭代次数。这样后续循环更安全。

function [x_denoised, imfs] = EMDdenoise(x, varargin) % EMDdenoise 对一维信号执行EMD去噪 % 输入： % x - 一维信号，列向量 % varargin - 'SiftTolerance', 阈值；'MaxNumIMF', 最大IMF数 % 输出： % x_denoised - 去噪后信号 % imfs - 分解出的IMF矩阵，每一行是一个IMF p = inputParser; addRequired(p, 'x'); addParameter(p, 'SiftTolerance', 0.2); addParameter(p, 'MaxNumIMF', []); parse(p, x, varargin{:}); x = p.Results.x(:); % 强制列向量 siftTol = p.Results.SiftTolerance; maxIMF = p.Results.MaxNumIMF; % 简单NaN处理：线性插值填补 if any(isnan(x)) t = (1:length(x))'; x = interp1(t(~isnan(x)), x(~isnan(x)), t, 'linear'); end

这里用inputParser统一参数管理。SiftTolerance默认0.2是工程上比较常用的起始值，如果你发现分解出的IMF没有意义，可以调小到0.1，但迭代次数会明显增加。MaxNumIMF为空表示不限制层数，实际使用时建议限制在4~8层，因为真实信号分解过头会把噪声也变成“伪IMF”。

3.2 核心筛分循环：包络线构造与停止判断

EMD主循环需要反复做“找极值→拟合包络→做差→判断”。这里的包络线可以用MATLAB的findpeaks找极值，再用spline做三次样条插值。下面是提取单个IMF的辅助函数：

function imf = extractIMF(signal, siftTol, maxSift) imf = signal; for iter = 1:maxSift % 找局部极大值和极小值 [pks, locMax] = findpeaks(imf); [valls, locMin] = findpeaks(-imf); locMin = locMin; vals = -valls; % 还原真实极小值 % 如果极值点太少，无法继续筛分 if length(locMax) < 2 || length(locMin) < 2 break; end % 构造上下包络 upEnv = spline([1; locMax; length(imf)], [imf(1); pks; imf(end)], 1:length(imf)); lowEnv = spline([1; locMin; length(imf)], [imf(1); vals; imf(end)], 1:length(imf)); meanEnv = (upEnv + lowEnv) / 2; % 原始IMF候选减去包络均值 h = imf - meanEnv; % 计算筛分收敛标准偏差 denom = sum(imf.^2); if denom == 0 break; end SD = sum((imf - h).^2) / denom; imf = h; if SD < siftTol break; end end end

这段代码里需要注意：spline端点处理很关键，直接在信号两端补了原始端点值，可以缓解端点效应，但无法完全消除。findpeaks默认要求极值点比其他点大，如果信号是纯噪声，可能会找到很多伪极值，所以后面在调用层会先用较小的MaxNumIMF限制。SD的分母用的是上一轮imf的能量，这比用原始信号能量更灵敏，不容易过早停止。

3.3 去噪策略：用相关系数或能量占比决定哪些IMF保留

分解得到IMF矩阵后，怎么区分哪些是噪声？常用的方法有两个：一是看IMF与原始信号的相关系数，真实成分通常相关性高；二是看IMF的过零率，高频噪声过零率很高。更稳定的是用相关系数做阈值。先计算所有IMF与去均值原始信号的相关系数，然后只保留相关系数超过某个阈值的IMF，并剔除落在噪声主导区间的最前面几个IMF。

imfCorr = zeros(size(imfs, 1), 1); xCentered = x - mean(x); for k = 1:size(imfs, 1) tmp = imfs(k, :) - mean(imfs(k, :)); denom = sqrt(sum(tmp.^2) * sum(xCentered.^2)); if denom == 0 imfCorr(k) = 0; else imfCorr(k) = sum(tmp .* xCentered) / denom; end end % 找出相关系数显著高于噪声层的IMF位置 threshold = 0.3; % 阈值需要根据实际信号调整 denoisedIMFIdx = find(imfCorr > threshold); x_denoised = sum(imfs(denoisedIMFIdx, :), 1);

这里的阈值0.3是经验值。如果信号本身比较干净，阈值可以提高到0.5；如果噪声很重，0.2可能更合适。还有一种做法：把所有IMF按相关系数从大到小排序，取前50%作为保留IMF，但这样容易把小幅度真实成分丢掉。我更倾向于先用相关系数排序，然后观察前几个IMF的时域波形，确认没有明显噪声振荡后再自动筛选。

3.4 完整函数框架与调用

把上面几段拼装成EMDdenoise.m，主流程就是：预处理、循环提取IMF、计算相关系数、重构信号。在example.m里，你可以这样调用：

x = load('signal.mat').x; fs = 1000; t = (0:length(x)-1) / fs; [x_clean, imfs] = EMDdenoise(x, 'SiftTolerance', 0.15, 'MaxNumIMF', 6); figure; subplot(2,1,1); plot(t, x); title('原始含噪信号'); subplot(2,1,2); plot(t, x_clean); title('EMD去噪信号');

这种方式适合离线分析，因为EMD本身是批量处理，不太适合实时流式去噪。实时场景通常用滑动窗口，每来一段时间窗做一次EMD，只取窗口后半部分的重构结果，以缓解边缘效应。

4. 运行example.m：把去噪流程跑通并量化效果

4.1 构造一段含噪的非平稳信号

没有现成采集数据时，先造一个已知真值的合成信号，这样能定量评估去噪效果。典型的非平稳信号是频率调制中频信号叠加趋势和冲击，再加白噪声：

fs = 1000; t = (0:1.5*fs-1)' / fs; trueSignal = sin(2*pi*(10 + 5*t).*t) + 0.5*sin(2*pi*3*t) + 0.3*(t>1); noise = 0.4 * randn(size(t)); x = trueSignal + noise;

这段代码生成一个频率随时间增加的正弦分量、一个3Hz低频分量以及一个在1秒后出现的阶跃趋势，噪声标准差是0.4。用这样的信号跑EMDdenoise，能同时检验EMD对频率漂移、趋势突变和随机噪声的分离能力。如果你用固定通带的低通滤波器，频率漂移部分会被削掉一部分，但EMD不会。

4.2 调用EMDdenoise并绘制原始/噪声/重构信号对比

在example.m中，我们按3.4节方式调用，并额外画出每个IMF：

[x_clean, imfs] = EMDdenoise(x, 'MaxNumIMF', 7); figure; subplot(4,1,1); plot(t, x); title('含噪信号'); subplot(4,1,2); plot(t, imfs(1,:)); title('IMF1 (高频)'); subplot(4,1,3); plot(t, imfs(end-1,:)); title('IMF4 (中低频)'); subplot(4,1,4); plot(t, x_clean); title('重构去噪信号');

观察重点：IMF1是否呈现明显的随机振荡，IMF2或IMF3是否能捕捉到频率调制的正弦部分，残余里是否有阶跃。如果IMF1和IMF2都像噪声，说明阈值设得太低，把噪声也保留了；如果去噪信号比真值平滑很多，说明阈值太高，丢掉了有效振荡成分。

4.3 用SNR、RMSE评估去噪效果

纯粹用眼睛看图不够，需要量化。常用指标是信噪比（SNR）和均方根误差（RMSE）。注意幅值单位要一致：

denoiseErr = trueSignal - x_clean'; snr_d = 20 * log10(norm(trueSignal) / norm(denoiseErr)); rmse_d = sqrt(mean(denoiseErr.^2)); fprintf('EMD去噪后 SNR = %.2f dB, RMSE = %.4f\n', snr_d, rmse_d);

有一次合成信号实验中，原始含噪信号的SNR约5.8 dB，去噪后能到13~15 dB，RMSE从0.28降到0.12左右。具体数值和多次平均有关，但趋势很稳定。如果去噪后SNR反而下降，先检查是不是信号里存在强冲击被EMD当成噪声去掉了，或者端点效应污染了低频IMF。

4.4 最容易踩的坑：端点效应与模态混叠

端点效应是EMD去噪中最常见的问题。信号两端的极值缺失，样条包络在端点处容易发散，导致首尾的IMF出现异常抖动。缓解办法有三种：端点处用镜像延伸法延长数据；只取信号中间80%的重构结果，丢弃两端10%；或者用小波包做预滤波后再做EMD。模态混叠则表现为一个IMF里同时包含明显不同时间尺度的成分，通常是因为信号中存在间歇性大幅扰动。遇到这种情况，最简单的是把MaxNumIMF调大一点，给高频部分更多的分解空间；更可靠的是换用EEMD或CEEMDAN，后续章节会讲。

5. 再往深走：让EMD去噪更可靠的几个实用调整

5.1 限制分解层数，避免过分解

不是IMF越多越好。当信号长度只有几百点，分解出七八个IMF时，后面几个IMF往往波形畸变。我一般把MaxNumIMF设为min(6, floor(log2(length(x)))-1)，这个经验值在多数振动和生物信号上表现不错。你也可以在分解后检查最后一个IMF的极值点数量，如果明显少于前面IMF，说明分解已经到趋势项了，后面的层可以丢弃。

5.2 用相邻IMF相关系数自动选阈值

固定阈值不好用，可以用相邻IMF相关系数的变化来自动判断噪声与信号的分界。计算imf(k)与imf(k+1)的相关系数，高频噪声之间的相关系数很低，而真实成分的IMF在相邻层之间往往存在显著相关性。当相关系数从低变高的转折点出现时，转折点后第一个IMF记为信号起始层：

adjCorr = zeros(size(imfs,1)-1, 1); for k = 1:size(imfs,1)-1 a = imfs(k,:) - mean(imfs(k,:)); b = imfs(k+1,:) - mean(imfs(k+1,:)); adjCorr(k) = sum(a.*b) / sqrt(sum(a.^2)*sum(b.^2)); end startIdx = find(adjCorr > 0.5, 1, 'first') + 1; x_denoised = sum(imfs(startIdx:end,:), 1);

这里0.5是参考值。如果信号中趋势分量很强，相邻IMF相关系数可能一直较高，那么可以从最后一层往前找首次低于0.3的位置作为噪声层边界。

5.3 当EMD不够用：CEEMDAN与EEMD的选择

模态混叠严重时，可以先用集合平均的方法，比如EEMD，给原信号加多次白噪声做EMD再平均。CEEMDAN在MATLAB里有第三方实现，要自己放入路径。它比EEMD收敛快，残余噪声更小。换用时只需要把EMDdenoise内部的分解循环替换为ceemdan(x, 'NumEnsembles', 200)，之后的相关系数筛选逻辑完全通用。代价是计算时间增加一个数量级，适合对离线的、特征微弱的信号做精细去噪。如果你的目标是嵌入式实时处理，建议不要在EMD上死磕，改用子带滤波器和小波硬阈值组合，速度会快得多。

本文还有配套的精品资源，点击获取