企业官网建设流程全解析

处理信号数据时，经常遇到一堆波形叠在一起的情况：有用信息埋在噪声下面，趋势项又盖住细节，直接看频谱图就像一锅粥。我用的环境是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)，得到第一版候选IMF：h1(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{i+1}; case 'SDTol' SDTol = varargin{i+1}; case 'MaxIMF' MaxIMF = varargin{i+1}; 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 end

2.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(i+1) locs(end+1) = i; pks(end+1) = 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(i+1) locs(end+1) = i; pks(end+1) = 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); end

diff(x)得到相邻两点的差值，对差值取符号，符号从正变负的位置就是局部极大值。这个版本比循环快很多，但如果信号有平台段，也就是连续几个点数值相等，diff会产生0，sign会返回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

这里有个细节：如果极值点数量少于2，interp1会直接报错，所以必须先做数量判断，返回空数组，让上层函数决定是终止还是跳过。很多网上的代码没处理这个边界情况，数据一长就崩。

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强制截断，避免死循环。如果数据里带了NaN，sum运算会全出问题，建议调用前先清洗数据。

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; end

4. 直接替换：把自己的数据跑起来

代码写完，最重要的一步就是让读者能用在自己的数据上。我设计时就坚持一个原则：拿过来，替换数据源，改一个采样率，剩下什么都别动，直接出图。

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的主峰落在50Hz，IMF2的主峰落在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 C，C个通道 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和50Hz，EMD会有点吃力，这是它的天然限制。

想缓解的话，把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的完整思路讲清楚了，能帮你在信号分解和对比实验里少走点弯路。