ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

EMD一维信号去噪实战:MATLAB完整流程与参数调优

EMD一维信号去噪实战:MATLAB完整流程与参数调优 简介MATLAB环境下的EMD去噪实现面向信号处理初学者、科研人员及相关工程开发者用于解决非线性、非平稳一维信号的噪声消除难题。资源包共2个文件均为m格式的MATLAB脚本核心函数负责EMD分解与信号重构示例脚本演示从数据预处理到结果可视化的完整流程压缩包整体仅6KB轻量易学。该资源已有335人学习适合快速入门。代码严格遵循经验模态分解的标准步骤定位局部极值、构造上下包络线、计算平均包络、提取内在模态函数IMF并通过希尔伯特变换分析瞬时频率最终叠加IMF与残余分量实现去噪。运行示例后可直观观察原始信号、各IMF分量及去噪信号的对比理解EMD的自适应分解特性。无论是用于课程实验、科研探索还是工程中的信号清洗这份紧凑的代码都能提供清晰的算法参考与可扩展基础。资源虽小但逻辑严密、注释到位是掌握EMD原理与MATLAB实现的便捷切入点。1. 从傅里叶到emd为什么一维信号去噪选它一维信号去噪的常规思路是先把信号变换到某个频域或小波域把噪声对应的系数压掉再逆变换回来。这套思路的前提是“信号和噪声在变换域可分”但真实场景里大量信号并不满足这个假设旋转机械的振动信号、脑电和心电的生理信号、结构健康监测里的应变数据频率成分随时间变化且噪声往往是宽带的短时傅里叶变换窗口一固定就顾此失彼小波变换则要提前选基函数和分解层数选错一层去噪结果就会把有用瞬态抹平。emd经验模态分解走的是另一条路它不需要任何先验基函数直接从信号自身的时间尺度特征出发把信号分解成一组固有模态函数imf和一个残差项。每个 imf 都满足两个条件——极值点数和过零点数至多相差 1且上下包络的均值处处为 0。这意味着分解是完全数据驱动的对非线性、非平稳的一维信号尤其好用。去噪逻辑也随之变得直观噪声通常落在高频 imf 上把高频分量剔掉或做阈值处理再叠加剩余分量就能得到干净的重构信号。适合的人群很广搞故障诊断的、处理生理信号的、做地质勘探数据预处理的都能直接套用这套流程。2. matlab 里的 emd 实现从内置函数到 eemd/ceemdan 的选型2.1 emd 算法在 matlab 中的三种存在形式先明确一点MATLAB 自身经历了两个阶段。旧版本比如 R2018a 之前里emd不在标准工具箱中很多工程师用的是法国学者 Pierre 等人发布的原始 MATLAB 程序包包含emd.m和若干依赖函数核心就是三次样条包络和筛分迭代。R2018a 之后MATLAB 在 Signal Processing Toolbox 里正式提供了emd函数这也让大量只装了基础工具箱的用户省去了找代码包的麻烦。但对于实际工程信号标准 emd 有个顽固的毛病——模态混叠。一个典型的场景是一段包含时变正弦成分的信号叠加了一个小幅值的高斯白噪声分解结果里某个 imf 可能同时包含了两个不同频率尺度的成分。解决这个问题有两个主流变体在 MATLAB 里也都有成熟的实现路径EEMD集合经验模态分解对原信号多次叠加白噪声后分别做 emd最后对同名 imf 取平均。白噪声的均匀频谱特性会迫使不同尺度的成分被分到正确的 imf 中。CEEMDAN自适应噪声完备集合经验模态分解在每次分解时只加入一次白噪声并且第一个 imf 计算完之后后续分解针对的是残差信号。它比 eemd 迭代次数少重构误差几乎为零。选型建议很直接信号大体平稳、信噪比不低用内置emd就够信号强非平稳、工程实测数据噪声大直接上 CEEMDAN。MATLAB 里 CEEMDAN 没有内置函数我用的是从 File Exchange 获取的第三方实现加载后调用格式和 eemd 几乎一致。2.2 在本地跑通最小可用代码如果你用的是新版 MATLAB打开命令行窗口直接执行下面的代码就能看到 emd 函数的基本行为fs 1000; % 采样率 1000 Hz t (0:999) / fs; % 1 秒时长 x sin(2*pi*50*t) 0.5*sin(2*pi*120*t) 0.3*randn(size(t)); [imf, residual] emd(x);代码里emd的默认停止条件由SiftMaxIterations等参数控制默认值对多数仿真信号都够用。imf是矩阵每一行对应一个固有模态函数第一行通常是最高频分量最后几行对应的是低频趋势residual是残差也就是单调趋势项。运行后你会在工作区看到imf的维度是 5x1000说明原信号被分解成了 4 个 imf 加 1 个残差。如果想看每个 imf 的能量占比可以执行energy sum(imf.^2, 2); plot(energy, o-);按我的习惯在真正做去噪之前会先用这个最小代码跑通一版确认分解行为符合预期再切换到真实信号上。因为 emd 的分解结果对参数选择很敏感先在小样本上把趋势看清楚能省掉后面不少调参的时间。2.3 找不到 emd 函数时的排查顺序很多人在论坛提问“为什么我的 MATLAB 没有 emd 函数”最常见的原因不是安装问题而是信号处理工具箱没有安装或没有激活。排查顺序是先在命令行执行ver(signal)如果返回空说明 Signal Processing Toolbox 没装如果确认装了但还是报Undefined function emd那大概率是版本太旧——R2017b 及更早版本确实没有内置emd。旧版本的一个临时方案是去下载第三方的 emd 工具包解压后执行addpath(genpath(emd_toolbox))把工具箱路径加入当前会话但这样每个新会话都要重新添加。更稳妥的做法是把路径写入startup.m文件这样 MATLAB 每次启动会自动加载。第三方包的质量参差不齐我遇到过同一份信号在不同版本包里分解结果不一致的情况建议优先升级 MATLAB 版本或安装 Signal Processing Toolbox其次是选用被引用较多的稳定版本封装。3. 用 emd 对一维信号去噪的完整流程与核心代码3.1 去噪的核心思路高频剔除与阈值收缩把 emd 用于去噪直接的想法是“前几个 imf 是噪声扔掉”但这种一刀切的做法很容易把真实信号中的瞬态高频成分也丢掉。工程上更稳妥的处理方式是分两步走先通过某种准则把 imf 分成“噪声主导”和“信号主导”两组再对噪声主导的 imf 做阈值收缩而不是彻底清零。我常用的一个准则是基于相关系数突变点计算每个 imf 与原信号的相关系数通常噪声 imf 的相关系数很低而第一个信号主导 imf 的相关系数会突然跳高。找到这个突变位置后这个位置之前的 imf 交给阈值收缩处理之后的 imf 直接保留。阈值收缩参考小波去噪中的软阈值思想对每个噪声主导 imf 计算其标准差再乘一个系数作为阈值。3.2 可直接运行的完整去噪脚本下面这段脚本覆盖了从信号读取到重构输出的完整链路可以直接复制运行我把它作为模板用了很久。% 1. 读入一维信号 % 假设 data 是列向量fs 为采样率 % 如果数据在 CSV 中可用data readmatrix(signal.csv); fs 1000; t (0:length(data)-1) / fs; % 2. 执行 EMD 分解 [imf, residual] emd(data, MaxNumIMF, 8, ... SiftMaxIterations, 100, ... Display, 0); % 参数说明 % MaxNumIMF 限制分解层数防止过分解出低频伪分量一般设 6~10 % SiftMaxIterations 是筛分迭代次数上限默认 100 通常够用 % Display 设为 0 可避免运行时间过长时刷屏 % 3. 计算相关系数并定位突变点 corr_vals zeros(size(imf, 1), 1); for k 1:size(imf, 1) temp corrcoef(data, imf(k, :)); corr_vals(k) temp(1, 2); end % 找相关系数陡增的位置从第二个分量开始找第一个明显跳变点 diff_corr diff(corr_vals); [~, idx] max(diff_corr); % 突变最大处 noise_imf 1:idx; % 前面这些视为噪声主导分量 % 4. 对噪声主导的 IMF 做软阈值收缩 denoised_imf imf; sigma_est median(abs(imf(idx, :))) / 0.6745; threshold sigma_est * sqrt(2 * log(length(data))); for k noise_imf y imf(k, :); y sign(y) .* max(abs(y) - threshold, 0); denoised_imf(k, :) y; end % 5. 重构去噪信号 denoised_signal sum(denoised_imf, 1) residual;参数说明第 3 步里用相关系数突变点做分组依据实际使用时要留意一种情况——如果信号本身能量极低相关系数可能普遍很小突变不明显。这时我会把判别准则换成对每个 imf 做功率谱峰值频率检查峰值频率高于某个先验值的才归为噪声组。第 4 步的标准差估计用的是一维中值估计这是从小波去噪里借来的经验公式0.6745对应高斯噪声的标准差与 MAD 的换算系数。sqrt(2*log(n))是通用阈值如果你想更激进地压制噪声可以把阈值系数提高到 1.2~1.5 倍。3.3 去噪效果的评价方法做完去噪必须量化评估不能只靠眼睛看波形。常见的两个指标是信噪比增益和均方根误差。如果原始干净信号已知直接计算SNR_in 10 * log10(sum(clean.^2) / sum((data-clean).^2)); SNR_out 10 * log10(sum(clean.^2) / sum((denoised_signal-clean).^2)); fprintf(SNR gain: %.2f dB\n, SNR_out - SNR_in);如果干净信号不可得就用平滑度指标或能量占比来间接评估。平滑度的定义是相邻差分平方和的比值去噪后平滑度应当在 1 附近且小于噪声信号。我还会顺手画出重构残差data - denoised_signal的时域图理论上残差里只能看到高频抖动如果残差里出现了明显的低频结构说明阈值收缩把真实信号的一部分也削掉了这时候要先降低阈值系数再检查是不是高频分组分错了。4. imf 筛选的三种实用方法怎么判断哪些分量该留4.1 方法一相关系数排序与突变点判别去噪的核心环节就是区分“信号 imf”和“噪声 imf”。相关系数法是最直观的噪声分量与原信号的相关性低信号分量相关性高。实现方式在第 3 章已经给出但有一个边界情况需要讨论当信号本身是窄带分量时中间的某个 imf 可能相关性和两边都差不多没有明显的突变点。处理这个边界情况我的折中方案是设定一个经验阈值相关系数低于 0.1 的 imf 一律视为噪声主导高于 0.3 的一律视为信号主导中间的根据上下文判断。这个阈值不是固定的信号采样率越低、信噪比越低阈值就应相应下调。4.2 方法二连续均方误差CMSE准则第二种方法和相关系数法完全独立适合做交叉验证。思路是计算相邻重构信号之间的均方误差第一个明显的极小值点对应的位置就是噪声主导和信号主导的分界。具体做法是先把所有 imf 叠加成一个完整重构然后逐个去掉排序靠前的 imf计算每次去掉前后的能量变化n_imf size(imf, 1); cmse zeros(1, n_imf-1); recon sum(imf, 1) residual; for k 1:n_imf-1 recon_k sum(imf(k1:end, :), 1) residual; cmse(k) mean((recon - recon_k).^2); end [~, cut] min(cmse);CMSE 准则背后的逻辑是当去掉的主要是噪声时重构信号变化剧烈CMSE 会有一个明显的峰值而去掉第一个信号分量时CMSE 往往跌到局部最小。cut就是信号分量的起点前面的全是噪声。把相关系数法和 CMSE 法的结果对比如果两者给出的分组不一致优先相信 CMSE因为它基于能量统计受信号波形影响更小。4.3 方法三过零率与瞬时频率联合判据第三种方法更细粒度适合信号中有明确物理含义的场合比如振动分析中轴承故障的特征频率已知。EMD 分解出的每个 imf 都应该对应一个窄带瞬时频率计算每个 imf 的过零率也就是单位时间内穿过零轴的平均次数然后转换成等效频率。如果某个 imf 的等效频率落在远高于关注频段的范围内而且带宽特别宽大概率是噪声残留。实际代码可以通过sum(abs(diff(sign(imf(k,:))))0) / (2*length(t)) * fs来估计。这个方法的额外好处是可以一次覆盖多个频段故障频率集中在某几阶时直接把这些频段对应的 imf 单独挑出来重构去噪的同时还做了特征提取一举两得。三种方法各有适用场景我在实际项目中通常的组合策略是方法适用场景不足相关系数突变点通用信号、信噪比尚可信号分量幅值过小时失效连续均方误差噪声较强、无明显频带特征对低频分量不敏感过零率瞬时频率有先验频段信息的工程信号需要人工设定频段范围5. emd 去噪的关键参数与故障排查从模态混叠到端点效应5.1 筛分停止条件的三个参数在用内置emd函数时几个参数直接决定分解质量。SiftMaxIterations是最大筛分迭代次数默认 100。迭代次数过小时包络均值没收敛imf 不满足正交条件过大时会把有效信号也磨平出现类似过拟合的效应。MaxNumIMF限制最大分解层数默认是log2(n)-1量级但对于短信号实际分解出的层数往往达不到上限。MaxNumExtrema限制了极值点数量阈值当剩余极值点少于该值时分解强制结束。我一般这样设[imf, residual] emd(x, ... SiftMaxIterations, 80, ... MaxNumIMF, 8, ... MaxNumExtrema, 2);这里把迭代上限设到 80 而不是默认的 100原因是工程信号普遍带噪过多的迭代会让包络跟随噪声的局部波动分解结果反而更碎。5.2 模态混叠的三个典型症状与破解方案模态混叠是 emd 去噪中最常遇到的坑症状非常典型一个 imf 的瞬时频率在某个时刻突然跳变或者相邻两个 imf 的频谱大面积重叠。破解方案按优先级排列增加噪声辅助分解——改用 eemd在 MATLAB 中实现时加入强度为原信号标准差 0.1~0.3 倍的白噪声集成次数设 100~200。这招能解决绝大多数轻微混叠。调整筛分迭代次数——迭代次数过少会造成欠筛分imf 中混入相邻尺度的分量适当提高到 150~200 试试。分频段预滤波——如果信号里有明显分离的强频率成分先用带通滤波器把关注频段切出来再对子带信号分别做 emd。这招虽然牺牲了“自适应”这个卖点但工程上极其好用信号成分本来就该先验地分离的场合没必要硬靠分解来分。5.3 端点效应与边界拟合的处理大量实测信号不是周期采样的截断边界处的极值点缺失会让三次样条包络在两端出现幅度渐大的摆动这就是端点效应。结果是最低阶的 imf 和残差项在两端明显偏离真实趋势。有两种可靠解法镜像延拓以边界极值点为对称轴把边界内的信号镜像映射到外部再构造包络。MATLAB 内置的emd已经实现了类似思路但处理的是默认模式如果你用旧版本第三方包需要自己在调用前做一次镜像延拓。极值延拓利用边界处的最后几个极值点通过多项式外推得到虚拟极值再参与包络构造。不过这要求你的采样数据在边界处没有大的突变否则外推极值会发散。验证端点效应是否被压住的方法很简单对比去噪信号和原始信号在两端各取 10% 长度的残差如果残差幅度明显大于中间段说明端点效应还在影响结果需要加强延拓或对边界段做额外平滑。5.4 去噪后信号失真的排查清单症状最可能原因排查动作重构信号幅值变小软阈值收缩把有效成分的峰值也削了降低阈值系数到 0.5~0.8波形整体偏移残差项被误删或阈值处理影响低频分量检查残差是否保留查看最低阶 imf 的频谱高频部分仍有毛刺分组时把第一个信号分量错判为噪声改用 CMSE 准则重新分组局部时段波形异常端点效应或模态混叠加上镜像延拓检查该时段瞬时频率6. 重构质量验证与零均值对称性技巧最后一层功夫放在重构质量的验证上这里有一个容易忽略的细节EMD 分解保证每个 imf 的均值为零或接近零但重构时如果把噪声主导的 imf 用硬阈值清零而不是软阈值收缩重构信号会出现微小的直流偏移。这个偏移单看波形不明显但对后续的积分运算比如振动信号积分到位移会被放大。我的做法是重构完成后计算mean(denoised_signal - data)。正常情况下这个差值均值应该在零附近如果明显偏离零就要检查是不是有 imf 被整体清零导致的。进一步验证分解质量可以用一个对称性检查——把去噪信号反向排列后重新做一次 emd两个方向分解出的 imf 个数和能量分布应该高度一致。如果差异超过 10%说明原始分解对端点敏感需要重新做延拓。最后一个实用技巧如果你处理的是一维信号流比如长时间连续监测数据不要整段做 emd而是用滑动窗口切段每段 1024 或 2048 个采样点窗口重叠 50%每段单独分解去噪后在重叠区做加权平均拼接。这么做既能实时出结果又能避免长序列分解带来的计算量爆炸和端点效应累积。重叠区的加权系数按汉宁窗来取拼接处几乎看不出接缝。把去噪前的信号和重构信号放在同一张图里对比再用第 3 章的 SNR 指标做定量确认这版 emd 去噪管线就算完整落地了。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进