处理信号数据时经常遇到一堆波形叠在一起的情况有用信息埋在噪声下面趋势项又盖住细节直接看频谱图就像一锅粥。我用的环境是Matlab但手头没有信号处理工具箱也不想为个分解算法专门装一堆依赖所以自己用纯Matlab基础函数写了一套EMD经验模态分解工具。这个工具不用工具箱、不预设分量个数把数据丢进去就能自动出原始信号图、分解效果图和频谱图替换成自己的数据也很简单。这篇我把实现思路、代码细节和踩过的坑一次讲清楚。这套代码最开始是为了处理一段旋转机械的振动信号。振动传感器采集到的数据往往是转频、倍频、故障特征频率外加随机噪声全部叠在一起直接用FFT看频谱有些靠得很近的峰根本分不清。我需要的不是再做一个滤波器而是先把信号拆成不同频率层次的波形再逐个分析。EMD非常适合干这个活。用惯了工具箱的人可能不理解为什么非要自己写。几个原因老版本Matlab不一定有内置emd函数同事之间传代码容易缺工具箱写论文或做课程设计时评委经常会问“你的核心算法怎么实现的”自写版本可以直接把流程摆出来更重要的是自写代码可以二次开发改成EEMD、CEEMDAN都很顺手对比实验也更公平。这个代码我在好几个项目里用过稳定可靠下面完整拆解。1. EMD到底在做什么先想明白再写代码很多人第一次接触EMD最困惑的是它和傅里叶变换、小波变换到底有什么区别。我常打一个比方你面前有一锅汤里面浮着花椒、辣椒段和葱段。傅里叶变换是去问“这锅汤里红色占比多少、绿色占比多少”只告诉你颜色比例不告诉你这些东西在锅的哪个位置。小波变换像用漏勺在不同位置去捞能留下“某时刻某层”的信息但捞的过程依赖你选什么尺寸的漏勺。EMD就不一样它直接徒手一层一层捞先把最表面的辣椒捞走再捞葱段最后剩下锅底沉淀。捞出来的每一层就是IMF固有模态函数。EMD是自适应的它不提前假设信号由哪些频率组成也不用你设置滤波器通带。算法会自己去寻找信号的局部极大值点和极小值点用它们构造上下包络再通过反复筛分把高频波动一层层剥离。整个过程是数据驱动的对不同形态的信号都能自动适应这正是它作为通用信号分解方法的底气。1.1 一次筛分是怎么完成的EMD的核心动作叫筛分sifting。假设有一段信号x(t)算法要做这几件事找到所有局部极大值点用三次样条插值连成上包络线。找到所有局部极小值点同样用三次样条插值连成下包络线。计算上包络和下包络的平均值m(t)。用原始信号减去m(t)得到第一版候选IMFh1(t) x(t) - m(t)。检查h1(t)是否满足IMF条件如果不满足对h1(t)重复上述过程直到满足为止。IMF条件用专业点的话说有两个第一在整个数据段内极大值点数量和极小值点数量之差最多为1意思是波形上下起伏比较规整第二任意时刻上包络和下包络的均值接近0意思是波形相对零线基本对称。设置这两个条件是为了保证分解出来的每个分量有明确的瞬时频率意义。“接近0”需要一个量化标准。实际操作中我用的是SD准则这也是经典论文里的做法SD sum((h_{k-1} - h_k).^2) / sum(h_{k-1}.^2)当相邻两次筛分结果的SD小于某个阈值时就认为筛分收敛了。阈值一般取0.2到0.3我常用0.25。阈值越小筛分越精细但计算量也越大。1.2 为什么要逐层分解而不是一次分完得到IMF1之后用原始信号减去IMF1得到残差r1(t)。对r1(t)再做同样的筛分得到IMF2再减得到r2一直重复下去。这个过程像剥洋葱一层一层往里剥越到后面频率越低最后的残差通常是单调趋势或者幅度极小的剩余项。逐层分解是EMD能自适应确定分量个数的关键。如果信号里只有两个主频加噪声通常分解出三到五个IMF就停了如果信号非常复杂可能分解出十几个。你不用提前告诉它“给我分解成6个”这也是它比小波、VMD好上手的地方。VMD需要预设模态数K和惩罚因子调参调得人头疼EMD把这一步省了。停止条件需要认真处理。我在代码里设了三种情况任何一个满足就退出剩余信号的极值点少于2个没法再做包络插值剩余信号近似单调再分下去也没有意义剩余信号幅度已经小于原始信号最大幅值的百万分之一能量基本耗尽。这三个条件必须都写上。只写一个不是卡死就是分过头。2. 不靠工具箱的Matlab实现核心函数怎么设计我写的核心分解函数叫emd_self输入是一段一维信号输出是IMF矩阵和残差向量。整个实现只用Matlab基础函数不依赖任何工具箱。下面把代码思路一步步拆开讲。2.1 主函数分解入口主函数负责初始化变量、循环调用筛分函数、判断停止条件。考虑到不同版本Matlab兼容性参数用varargin解析不用较新的arguments块。function [imfs, residual] emd_self(x, varargin) % EMD_SELF 自写经验模态分解 % 输入 % x - 一维信号行向量或列向量均可 % 可选参数 % MaxSift - 最大筛分次数默认200 % SDTol - 筛分停止阈值默认0.25 % MaxIMF - 最大IMF数量默认20 % 输出 % imfs - IMF分量矩阵每行一个IMF % residual - 残差向量 % 参数解析 MaxSift 200; SDTol 0.25; MaxIMF 20; for i 1:2:length(varargin) switch varargin{i} case MaxSift MaxSift varargin{i1}; case SDTol SDTol varargin{i1}; case MaxIMF MaxIMF varargin{i1}; end end x x(:); N length(x); imf_cell {}; residual x; k 0; while k MaxIMF [imf_new, success] siftOnce(residual, MaxSift, SDTol); if ~success break; end k k 1; imf_cell{k} imf_new; residual residual - imf_new; % 终止条件1极值点不足 if length(findLocalMax(residual)) 2 || length(findLocalMin(residual)) 2 break; end % 终止条件2残差能量太小 if max(abs(residual)) 1e-6 * max(abs(x)) break; end end imfs zeros(k, N); for i 1:k imfs(i,:) imf_cell{i}; end end2.2 找极值点自己动手替代findpeaks很多号称“不用工具箱”的教程偷偷用了findpeaks函数这个函数其实属于Signal Processing Toolbox不是Matlab基础函数。为了让代码真正零依赖我写了一个简单的循环版局部极值搜索函数。function [pks, locs] findLocalMax(x) n length(x); pks []; locs []; for i 2:n-1 if x(i) x(i-1) x(i) x(i1) locs(end1) i; pks(end1) x(i); end end end function [pks, locs] findLocalMin(x) n length(x); pks []; locs []; for i 2:n-1 if x(i) x(i-1) x(i) x(i1) locs(end1) i; pks(end1) x(i); end end end这个写法非常直白对大多数信号足够用。如果数据量特别大可以用diff和find的组合优化性能function [pks, locs] findLocalMaxFast(x) d diff(x); signChange diff(sign(d)); locs find(signChange 0) 1; pks x(locs); enddiff(x)得到相邻两点的差值对差值取符号符号从正变负的位置就是局部极大值。这个版本比循环快很多但如果信号有平台段也就是连续几个点数值相等diff会产生0sign会返回0符号变化位置就不准了。所以稳健性优先的话先用循环版本保正确真有性能瓶颈再换快速版。2.3 包络插值三次样条与端点处理包络生成是整个EMD最容易出bug的地方。我用的interp1配合spline三次样条。但三次样条在首尾两端没有极值点约束会把包络甩出很大一个弯这就是著名的端点效应。缓解办法是手工补点把信号首尾值加入极值点序列强制包络从端点出发。function [upper, lower] envelopeXY(x) N length(x); [pksU, locU] findLocalMax(x); [pksD, locD] findLocalMin(x); if length(locU) 2 || length(locD) 2 upper []; lower []; return; end % 首尾补点缓解端点效应 locU [1, locU, N]; pksU [x(1), pksU, x(N)]; locD [1, locD, N]; pksD [x(1), pksD, x(N)]; upper interp1(locU, pksU, 1:N, spline); lower interp1(locD, pksD, 1:N, spline); end这里有个细节如果极值点数量少于2interp1会直接报错所以必须先做数量判断返回空数组让上层函数决定是终止还是跳过。很多网上的代码没处理这个边界情况数据一长就崩。2.4 筛分函数循环迭代的判断逻辑筛分函数是算法的发动机它负责把一段信号反复减去包络均值直到满足IMF条件。退出条件有两个一个是SD小于阈值一个是迭代次数达到上限。function [imf, success] siftOnce(x, MaxSift, SDTol) prev x; sd Inf; count 0; success true; while sd SDTol count MaxSift [up, down] envelopeXY(prev); if isempty(up) || isempty(down) success false; break; end m (up down) / 2; h prev - m; denom sum(prev.^2); if denom eps success false; break; end sd sum((prev - h).^2) / denom; prev h; count count 1; end imf prev; end我调试中发现有些信号筛分到后期SD值会在某个水平来回震荡不会继续下降。这种情况就靠MaxSift强制截断避免死循环。如果数据里带了NaNsum运算会全出问题建议调用前先清洗数据。3. 三类图一次性全部画出来很多人要的其实不只是分解结果而是能直接用在报告里的图。我的绘图函数plotEMDResult一次生成三张图原始信号图、分解效果图、频谱图。三张图各司其职从波形、分量组成、频率分布三个角度展示信号。3.1 原始信号图与时间轴画图的第一步是把时间轴算对。很多新手plot(x)之后横轴是采样点数如果采样率是1000Hz横轴刻度完全不对应秒。正确做法是先用采样率构造时间向量。fs 1000; t (0:length(x)-1) / fs;如果数据里没有采样率信息就把fs设为1横轴代表“每个采样点”频谱图横轴是“每个采样点周数”单位不是Hz。这点心里要有数别分析到一半搞混。原始信号图通常用细线、浅色背景加网格标题写清楚直接能贴报告。代码实现不复杂就是plot加一点格式化。3.2 分解效果图每个IMF占一行分解效果图是整个可视化里信息量最大的一张。行数等于IMF数量加1个残差每一行画一个IMF按顺序从上到下排列。rows size(imfs,1) 1; figure(Name,EMD分解结果,Color,w); for i 1:size(imfs,1) subplot(rows,1,i); plot(t, imfs(i,:), LineWidth, 0.8); ylabel([IMF num2str(i)]); xlim([t(1) t(end)]); set(gca,XTickLabel,[]); grid on; end subplot(rows,1,rows); plot(t, residual, r, LineWidth, 1.0); xlabel(时间/s); ylabel(残差); xlim([t(1) t(end)]); grid on;子图横轴范围统一方便对比同一时间段内各IMF的波动情况。y轴不统一因为不同IMF的幅值可能差几个量级统一了反而看不清低幅值分量。3.3 频谱图每个IMF的主频一目了然频谱图是整个可视化里最值钱的。每个IMF做完FFT后单独放在一个子图里横轴是实际频率纵轴是幅值能直观看到每个分量集中在哪个频段。这样分解结果和频率特征就能对上号。function plotSingleSpectrum(y, fs) N length(y); nfft 2^nextpow2(N); f (0:nfft/2-1) * fs / nfft; Y abs(fft(y, nfft)); Y Y(1:nfft/2) * 2 / N; plot(f, Y); xlim([0 fs/2]); end几个要点必须强调fft得到的是双边谱频率从0到fs我们通常只看单边所以只取前一半幅值要乘以2/N才是真实幅值否则plot出来的数值不是信号实际幅度补零到2的幂次会平滑曲线但不会提高频率分辨率。真实分辨率只取决于数据长度N和采样率fs的比值这个不要搞错。3.4 一键出图封装把三个图的代码封装成一个函数调用时非常清爽。我在实际项目里就是一行代码出全部图。function plotEMDResult(t, x, imfs, residual, fs) N length(x); % 图1原始信号 figure(Name,原始信号,Color,w); plot(t, x, LineWidth, 1.0); grid on; xlabel(时间/s); ylabel(幅值); title(原始信号); % 图2分解效果 rows size(imfs,1) 1; figure(Name,EMD分解结果,Color,w); for i 1:size(imfs,1) subplot(rows,1,i); plot(t, imfs(i,:), LineWidth, 0.8); ylabel([IMF num2str(i)]); xlim([t(1) t(end)]); set(gca,XTickLabel,[]); grid on; end subplot(rows,1,rows); plot(t, residual, r, LineWidth, 1.0); xlabel(时间/s); ylabel(残差); xlim([t(1) t(end)]); grid on; % 图3频谱 figure(Name,EMD频谱,Color,w); for i 1:size(imfs,1) subplot(rows,1,i); plotSingleSpectrum(imfs(i,:), fs); ylabel([IMF num2str(i)]); set(gca,XTickLabel,[]); grid on; end subplot(rows,1,rows); plotSingleSpectrum(residual, fs); xlabel(频率/Hz); ylabel(残差); grid on; end4. 直接替换把自己的数据跑起来代码写完最重要的一步就是让读者能用在自己的数据上。我设计时就坚持一个原则拿过来替换数据源改一个采样率剩下什么都别动直接出图。4.1 跑通一个模拟混合信号先造一个混合信号验证整个流程。我用的是50Hz正弦加10Hz正弦加随机噪声加趋势项这种信号非常接近真实场景里的“多频率叠加加背景噪声加趋势漂移”。fs 1000; t (0:1999) / fs; x sin(2*pi*50*t) 0.6*sin(2*pi*10*t) 0.2*randn(size(t)) 0.5*t; [imfs, residual] emd_self(x, MaxSift, 200, SDTol, 0.25); plotEMDResult(t, x, imfs, residual, fs);跑完之后分解效果图应该能看到IMF1基本是50Hz正弦IMF2是10Hz正弦残差接近趋势项。频谱图里IMF1的主峰落在50HzIMF2的主峰落在10Hz对应关系一目了然。这个例子也演示了“自动确定分量个数”的特点。你不需要告诉它“给我分解成3个IMF”它自己根据信号复杂度找到了合适的层数。4.2 换自己的数据三步完成替换“直接替换”具体操作只有三步第一读入数据。Matlab里读数据的函数很多readmatrix、csvread、load等都可以关键是最后得到一个数值向量x。如果数据是多列的选中信号所在的那一列。第二设置采样率。找到采集设备的参数把fs改成实际值。如果数据文件里没写先按照时间戳算一下两个采样点时间间隔的倒数就是采样率。第三调用函数。把x和fs传给emd_self和plotEMDResult。以读取csv为例data readmatrix(my_signal.csv); x data(:, 2); % 假设第二列是信号 fs 1000; % 根据采集设备参数修改 t (0:length(x)-1) / fs; [imfs, residual] emd_self(x); plotEMDResult(t, x, imfs, residual, fs);就这三步没有任何多余配置。我自己在项目里换新数据基本上改一行读取路径和一个采样率就够了。4.3 多通道信号怎么处理脑电、肌电、振动监测这类场景数据往往是多通道的每个通道一路信号。EMD本身是单通道算法不能把多通道直接揉在一起分解那样会丢掉通道间的空间信息。正确做法是每个通道单独分解把结果存在cell数组里。channels size(X, 2); % X是 N x CC个通道 all_imfs cell(1, channels); for ch 1:channels [imfs, residual] emd_self(X(:, ch)); all_imfs{ch} imfs; end如果想看同一阶IMF在不同通道的分布可以把第k个IMF取出来按通道排列画在一个大图里。这样能直观看出不同通道在同一频段的联动关系。4.4 作为对比方法的典型用法EMD经常在论文里当对比方法这个场景有两个容易踩的坑。第一个坑是控制变量。如果对比方法A用了内置工具箱EMD是自写代码别人可能会质疑实现标准不一致。我在写论文时会明确写清楚EMD采用经典sifting实现参数设置如下代码附在附录。这样审稿人能复现就不会揪着实现细节不放。第二个坑是量化指标。EMD输出的是IMF分量不是最终的分类或回归结果中间还要补特征提取。常见的特征有能量占比每个IMF能量占总能量的比例反映该频段贡献大小样本熵或排列熵衡量IMF序列复杂度峭度衡量冲击特征轴承故障诊断里很有用瞬时频率均值反映IMF的中心频率。以能量占比为例计算逻辑很直观energy_total sum(sum(imfs.^2, 2)); energy_ratio sum(imfs.^2, 2) / energy_total;把这些特征拼成一行向量后面接分类器或者回归模型就是经典的“EMD特征机器学习”方案。我在齿轮和轴承故障诊断里用过这个套路作为baseline很稳。5. 换数据实测踩坑记录这部分是干货。我把自己实际调试中遇到的典型问题、排查过程和解决办法整理出来每个问题都是真实踩过的。5.1 代码运行了很久不结束遇到过两次。一次是数据里有NaN一次是筛分SD值不收敛。NaN会让sum操作变成NaN循环条件判断永远不满足直接卡死。排查时先看数据简况isnan(sum(x))就能暴露问题。SD值不收敛通常发生在信号存在明显突变点或者强脉冲的场景。解决办法是设MaxSift上限我默认200次但还是卡住的话可以在筛分循环里加一个SD连续多次上升就提前终止的判断if count 3 sd sd_prev early_stop_count early_stop_count 1; if early_stop_count 20 break; end else early_stop_count 0; end这样避免了无意义的迭代分解结果也不会差太多。5.2 分解出来的IMF1跟原始信号长得差不多这种情况通常说明信号本身频带很窄比如一个比较干净的正弦波EMD把最高频成分一次就摘走了剩下残差几乎为零。这不是bug是EMD的正常特性。如果你是想把两个频率很近的分量分开比如48Hz和50HzEMD会有点吃力这是它的天然限制。想缓解的话把SDTol调小到0.1筛分更精细有一定帮助。更彻底的办法是用EEMD在信号里注入白噪声辅助分解能明显改善模态混叠。5.3 首尾两端出现明显发散端点效应是EMD的老毛病。三次样条插值在首尾两端缺少极值点约束包络会甩出去一个大弯。我处理过三种方案首尾补点把信号端点加入极值点集合。简单有效大多数情况下够用。镜像延拓把端点看成镜面将极值点对称复制到信号外侧再插值包络。效果更好但代码复杂度高。分析时截掉两端不要使用首尾各10%的数据只取中间80%做后续分析。最省事适合数据量充足的情况。做Hilbert谱分析时端点效应会影响瞬时频率计算这时候我会用镜像延拓。普通分解做特征提取首尾补点就够了。5.4 频谱图横轴和实际频率对不上绝大多数情况是采样率传错了。采样率1000Hz但画图时fs传成100所有频率都缩小10倍。还有一种情况是忘了取单边谱横轴从0画到1000看起来频率翻倍。检查这两处基本能解决。还有一个容易忽略的是频率分辨率。fs/N是分辨率如果数据长度N太短两个靠得很近的频率峰可能连成一个。这时候要增加采样时长不是靠FFT补零。补零只是让曲线变平滑不改变分辨率这个常识经常被误解。5.5 自写结果和内置工具箱结果不一致新版Matlab自带emd函数我拿自己的结果和它对比过主体IMF基本一致边界细节有差异。原因在于内置实现用了更精细的端点处理和终止标准内部参数不一定公开。这是正常现象不同软件包的EMD结果本来就有差异受极值插值方式影响很大。做复现实验时论文里写清楚“经典EMD实现默认参数”就足够。别人如果对不上第一件事应该检查数据长度、采样率、SDTol和端点处理方法是不是一致。6. 自写EMD的边界与变体扩展自写代码最大的好处是可以随意改造成各种变体。这里提两个最常见的扩展我已经在别的项目里验证过。EEMD的思路简单说在原始信号上多次叠加白噪声每次做一次EMD最后把同阶IMF取平均。白噪声的作用是让不同尺度的信号自动映射到合适的尺度上显著缓解模态混叠。代价是计算量成倍增加。如果只是把IMF当特征用EEMD通常够用。CEEMDAN是EEMD的改进版它在每次分解的残差里加入特定噪声收敛更快重构误差更小。核心逻辑和自写代码不冲突就是在筛分主循环外面再套一层噪声副本的管理。如果你后面要做信号重构CEEMDAN是更稳的选择。还有VMD变分模态分解思路跟EMD完全不同需要预设模态数K和惩罚因子调参很麻烦。EMD是“自动分层”VMD是“指定分几层再优化每层中心频率”。我做实验通常把两者都放进去对比各有优劣。EMD的优点是参数少缺点是模态混叠VMD恰恰相反参数多但分离效果更可控。另外EMD分解完之后的IMF不只是拿来看的还能做很多事去趋势直接把残差丢掉用所有IMF重构出去趋势的信号去噪把能量占比很低、又集中在高频的IMF丢掉剩下的重构出平滑信号特征输入把各IMF的熵、能量比、峭度等指标拼成特征向量喂给分类器或者回归模型。我自己做过“EMD去噪加LSTM预测”的组合先对原始序列做EMD把高频噪声IMF剔除剩余分量重构后作为LSTM的输入预测精度比直接用原始数据好一些。原因不难理解EMD把不同时间尺度的波动分离出来模型更容易捕捉到各自的规律。做SOC估计或者电价预测时这个套路值得一试。最后分享一个个人习惯。每次拿到新数据我不会直接上全套分析而是先跑一遍EMD看分解效果和频谱图心里对信号成分有个底。这就跟看病先量体温一样是诊断的第一步。这套自写代码已经在振动信号分析、电池数据预处理和课程设计里反复用过稳定可靠遇到问题也知道去哪排查。你拿过去跑你自己的数据如果遇到我这个清单里没写过的新问题欢迎根据自己的场景继续加判断条件。EMD这东西代码看着简单但每个细节都有讲究。把底层逻辑吃透了不管是改EEMD还是接深度学习你都能举一反三。希望这篇把MATLAB下自写EMD的完整思路讲清楚了能帮你在信号分解和对比实验里少走点弯路。