简介面向需要在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 ToleranceSD sum((h_{k-1} - h_k)^2) / sum(h_{k-1}^2)当SD小于阈值时认为筛分收敛。典型阈值在0.2到0.3之间。MATLAB内置的emd函数也提供类似参数常见设置如下参数典型值作用SiftRelativeTolerance0.2控制筛分收敛越小越严格MaxNumIMF空默认全分解限制最多输出IMF个数防止过分解MaxSiftIterations100单次筛分的最大迭代次数Interpolationpchip包络插值方式spline更平滑但可能过冲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*(t1); 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 dBRMSE从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(k1)的相关系数高频噪声之间的相关系数很低而真实成分的IMF在相邻层之间往往存在显著相关性。当相关系数从低变高的转折点出现时转折点后第一个IMF记为信号起始层adjCorr zeros(size(imfs,1)-1, 1); for k 1:size(imfs,1)-1 a imfs(k,:) - mean(imfs(k,:)); b imfs(k1,:) - mean(imfs(k1,:)); 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上死磕改用子带滤波器和小波硬阈值组合速度会快得多。本文还有配套的精品资源点击获取