简介：一份基于灰狼算法优化特征模态分解的Matlab完整源码与数据包，面向机械故障诊断、振动信号处理与模式识别领域的研究人员、工程师及高年级本科生，旨在解决传统FMD算法中滤波器长度、模态数等关键参数依赖人工经验试凑、难以快速收敛到最优解的问题。代码将GWO全局搜索与FMD分解相结合，内置包络熵、信息熵等6种可选目标函数，可根据信号特性自动寻优；完整实现信号自适应分解，并绘制原始信号、优化迭代曲线、参数变化过程及分解后的各模态分量，便于直观分析算法效果。压缩包共8个文件，其中6个.m脚本覆盖主程序、FMD核心算法、GWO优化器、初始化配置、目标函数计算和模态绘图模块，另含2个xlsx表格分别存放测试数据与分解结果，整包仅61KB，轻量易用、开箱即跑。目前已有216人学习使用，尤其适合作为FMD类算法复现、超参数对比研究及毕业论文实验部分的参考样例。

1. GWO-FMD：为什么特征模态分解要把参数交给灰狼算法去搜

做旋转机械故障诊断的人，大概率都遇到过这个场景：手里有一截轴承振动信号，噪底很大，转频和故障特征频率挨得又近，直接用FFT看频谱，谱线上全是毛刺，特征峰埋在里面。于是想到用模态分解把信号剥开。用过VMD的都知道，VMD的惩罚因子α和模态数K基本靠试，试错成本高，换一组数据又要重新调。FMD（Feature Mode Decomposition，特征模态分解）比VMD更聚焦——它通过频带划分和迭代滤波把特征分量分离出来，但FMD同样不是开箱即用，滤波器长度、步长、频带划分方式都会明显影响分解质量。GWO-FMD这套Matlab源码的思路，就是用灰狼优化算法去自动搜索FMD的关键参数，把原本靠经验的调参过程变成一次有约束的寻优。这篇文章我会从FMD和GWO的原理讲起，再落到Matlab源码包的目录结构、运行步骤、参数设置和避坑记录。适合两类人：一类是刚接触特征模态分解、想找一个可直接跑通的Matlab例程的研究生；另一类是已经在用VMD/EEMD做故障诊断、想对比FMD效果并省去手动调参的工程师。

2. 先看懂FMD与GWO的组合逻辑：从频带划分到狼群寻优的完整链路

2.1 特征模态分解的核心机制：为什么它能分离混叠分量

FMD的核心思想是把信号分解看成一次“滤波—重构—再滤波”的迭代过程。每一轮迭代，算法根据当前信号的频谱特征设计一组带通滤波器，把信号划分成若干个频带区间，然后通过滤波和重构得到一组模态分量。关键点在于：FMD的滤波器初始中心频率不是随便定的，而是基于原始信号自相关谱或功率谱的峰值位置来初始化，这让它在处理密集频谱分量时比VMD更有针对性。VMD在频域里强制各模态中心频率分离，遇到两个分量间隔小于带宽时容易把一个真实分量劈成两半；FMD用迭代估频的方式去逼近真实分量位置，因此对模态混叠的容忍度更高。

FMD的迭代过程里有几个核心参数：滤波器长度L、步长（也叫带宽步长或滤波步进）、惩罚系数、最大模态数。其中滤波器长度L直接影响频带分辨能力，L太小则频带分割粗糙，L太大则计算量成倍上涨且容易过拟合噪声。步长决定每次更新滤波器中心频率时的移动距离，步长越小收敛越慢但精度更高，步长大了则可能跳过最优频带。惩罚系数控制模态分量的平滑度，偏向于抑制噪声还是保留细节。这几个参数相互耦合，手动调参很难一次到位，这正是引入GWO的理由。

2.2 灰狼优化算法在做什么：它替人试了哪些参数

灰狼优化（Grey Wolf Optimizer，GWO）模拟灰狼群体的等级制度与狩猎行为。种群里有四个等级：α狼是头狼，β狼和δ狼是次级领导，ω狼是底层个体。优化过程分为包围、追捕、攻击三个阶段。每一轮迭代，ω狼根据α、β、δ三只头狼的位置来更新自己的位置，这样种群不至于只围绕一个局部最优解打转，三个头狼的加权指引相当于一个“小范围并行搜索”，不容易在早期就陷入局部收敛。

放在GWO-FMD这个源码包的具体语境里，GWO要搜索的参数通常是FMD的滤波器长度L、步长、惩罚系数，以及FMD的频带划分控制参数。每一只狼的位置就对应一组FMD参数组合。适应度函数的选择直接决定优化方向。

提示：GWO-FMD里最常见的适应度函数是包络熵。包络熵越小，说明分解出的模态分量包络波形越稀疏、冲击特征越明显，正好符合故障诊断中“提取冲击成分”的需求。

用伪代码表达GWO-FMD的嵌套逻辑如下：

%% GWO-FMD 两层嵌套结构 % 外层：GWO优化器 for iter = 1 : MaxIter for wolf = 1 : NPop params = Pos(wolf, :); % 一组FMD参数 fitness(wolf) = FMD_Fitness(params, x); % 包络熵评估 end % 更新 alpha/beta/delta 三只头狼 % 更新所有狼位置 end % 内层：FMD完整分解流程 function fitness = FMD_Fitness(params, x) modes = FMD(x, params); % 完整FMD分解 env = abs(hilbert(modes)); % Hilbert解调取包络 fitness = EnvelopeEntropy(env); % 包络熵越小越优 end

参数对应关系：NPop是灰狼种群数（15～30），MaxIter是最大迭代次数（30～80），Pos(wolf,:)是每个个体的参数向量。实际项目中，GWO搜索的参数上限和下限要根据信号采样率和目标频带提前卡好，不能直接用默认范围，否则狼群会跑到不合理频段去。

2.3 为什么不用网格寻优或贝叶斯优化

FMD的参数搜索空间是连续且非凸的，如果只用网格搜索，滤波器长度L取20～200，步长取0.01～1，每种组合再跑一次完整分解，需要几百次迭代，而每次FMD分解在长信号上要数十秒，总耗时不可接受。贝叶斯优化虽然采样效率高，但它依赖代理模型的拟合质量，对FMD这种参数与结果之间关系不光滑、噪声较大的问题，代理模型容易失真。GWO属于无梯度优化，不假设目标函数的解析形式，只需要能算出适应度值就能跑，实现也简单。更现实的理由是——源码包里已经配好了完整的GWO优化循环和适应度函数，拿到手直接改参数就能跑，不用额外装优化工具箱，这对Matlab版本比较老或没有Global Optimization Toolbox的环境尤其友好。

3. Matlab源码包结构：核心函数、测试数据与第一遍跑通的步骤

3.1 下载后第一件事：先看清目录里有什么

拿到源码包后先把压缩包完整解压，不要直接在压缩包内双击.m文件运行。目录结构大致是：

我建议先从main_GWO_FMD.m开始读，它决定了整条流程的骨架。用Matlab打开后，不要急着点“运行”，先把前四十行读完，确认信号和数据文件路径没有问题。这套源码的数据文件通常是.mat格式，文件名类似bearing_signal.mat，加载方式和数据字段名要看代码里写的是load还是matfile，版本差异可能导致字段引用失败。

3.2 第一遍运行：保持默认参数跑通全流程

大多数情况下，源码包自带的仿真信号是叠加了高斯白噪声的故障信号，FMD默认参数和GWO默认参数组合是能直接跑出结果的。第一遍运行的目的只有一个：确认整个链路没有报错，能画出图来。这里给出常见的启动流程。

%% 第一遍跑通全流程 clc; clear; close all; % 加载测试数据 load('data/bearing_signal.mat'); % 得到变量 x, fs fs = 12000; % 采样率，按实际数据修改 N = length(x); % 信号长度 t = (0:N-1)/fs; % 时间轴 % 调用GWO-FMD主函数 % Pop=20, MaxIter=30 是低成本的快速验证组合 [BestPos, BestFitness, ConvergenceCurve] = main_GWO_FMD(x, fs, 20, 30); % 输出最优参数 fprintf('最优滤波器长度 L = %.2f\n', BestPos(1)); fprintf('最优步长 lambda = %.4f\n', BestPos(2)); fprintf('最优惩罚系数 tau = %.2f\n', BestPos(3));

这里main_GWO_FMD接收四个输入：信号x、采样率fs、种群数20和迭代次数30。输出BestPos是最优参数向量，BestFitness是最小包络熵，ConvergenceCurve是每轮迭代的最佳适应度值，用来画收敛曲线判断优化是否有效。关于参数初值：种数取20、迭代取30是先验证流程的低成本组合，完整实验建议改到Pop=30, MaxIter=50，效果会有明显提升，但耗时也从几分钟级别跳到十几分钟级别。

3.3 用GWO返回的最优参数做正式分解

GWO优化完成后，得到的只是一组最优参数，真正要用的还是FMD函数本身。需要把最优参数传给FMD重新跑一次完整分解，并把各模态分量画出来。这部分代码在demo_plot.m里已经写好了，但很多新手会直接在GWO.m内部调用FMD，然后试图在GWO的函数内部绘图——结果要么图很多很乱，要么变量作用域不对报错。正确做法是把寻优和分解解耦：

%% 用最优参数执行最终分解 params.L = round(BestPos(1)); % 滤波器长度取整 params.lambda = BestPos(2); % 带宽步长 params.tau = BestPos(3); % 惩罚系数 params.K = 5; % 目标模态数，按信号先验设置 % 调用FMD分解 modes = FMD(x, fs, params); % 绘制前四个模态的时域波形 figure; for i = 1 : 4 subplot(4,1,i); plot(t, modes(i,:), 'b'); ylabel(['IMF', num2str(i)]); xlim([0, 0.1]); % 截取前0.1秒，观察冲击细节 end xlabel('时间/s');

逻辑说明：GWO搜索时为了效率，通常在适应度函数里限制FMD的迭代次数并使用简化参数，因此最优参数下正式分解时，要把params.K、params.MaxIter_FMD这些分量参数重新设置。modes是一个K×N矩阵，每一行是一个模态分量，绘图时截取前0.1秒是因为早期冲击特征最明显，完整信号太长会淹没在缩略图的压缩里。

注意：params.K不要设置得过大。FMD最突出的问题是模态数高于真实分量数时，会把单个物理分量拆成多段，尤其在噪声较强时。建议先跑一次K=3～6的区间，对比各模态的中心频率和波形特征，再锁定最终K。

4. 五个核心参数的实战取舍：频率范围、种群大小、迭代次数与适应度设计

4.1 划分频带边界：滤波器的初值决定搜索空间

GWO-FMD的第一个坑往往不是算法本身，而是搜索边界没设好。FMD内部会先基于信号频谱确定初始化频带，但如果GWO的搜索边界比FMD的实际频带范围大很多，狼群会浪费大量迭代在无效区域，甚至收敛到频带外的假最优。常见做法是根据信号频谱特征手动限定优化变量的上下界。

%% GWO-FMD参数边界设置 VarMin = [20, 0.01, 10]; % [最小滤波器长度L, 最小步长lambda, 最小惩罚系数tau] VarMax = [200, 0.5, 100]; % 上界 %% 频带边界限定（以采样率为基准） % 滤波器长度L的上界由信号频率分辨率决定： % Lmax = fs / df，其中df为预期的频带最小带宽 % 比如fs=12000Hz，目标频带最窄200Hz，则Lmax=60 L_upper = fs / 200; VarMax(1) = min(200, L_upper); %% 步长的物理含义 % lambda = 单次迭代滤波器中心频率移动比例 % lambda过小<0.01，迭代缓慢; lambda过大>0.5，跳过窄带特征分量

参数说明：滤波器长度L理论上限不应超过信号长度的一半，实际比这严格得多——L决定了FIR滤波器阶数，L越大频带边缘越陡但数值稳定性越差。频带下边界设置可以参考目标故障特征频率或转频的倍频。轴承内圈故障特征频率通常在转频的整数倍附近，如果先验知道特征频率在500～2000Hz，就可以把搜索空间压到这一段。

4.2 种群大小和迭代次数：算力与精度的权衡

种群数NPop从10到50之间影响最大的是收敛稳定性和计算时间。实测经验：低于15时，收敛曲线波动大且不稳定，每次运行结果差异明显；在20～30之间，收敛结果基本稳定；50以上提升有限但耗时翻倍。迭代次数同理，30次的收敛曲线可能还在下降，50次基本平了，100次以上通常只是边际收益。

这里有个容易被忽略的细节：GWO每一轮迭代要评估NPop个个体的适应度，而每个适应度都包含一次完整FMD分解。也就是说，总FMD调用次数=NPop×MaxIter。如果FMD在1万点信号上需要5秒，那么20×50=1000次调用就是83分钟。所以第一次跑通一定要用小种群和少迭代，确认正确后再把参数拉满。

4.3 适应度函数不一定非要包络熵：多指标加权设计

包络熵能反映冲击特征的稀疏性，但它对周期性冲击和随机噪声的区分度不够。如果分解出来一个模态全是低幅值噪声，它的包络熵可能也很低，因为包络幅值均匀。这时GWO会误判成最优解。一种常见改进是让适应度函数同时考虑峭度和包络熵，或者说用“包络熵与峭度的比值”作为目标。

%% 复合适应度函数：包络熵 + 峭度调节 function fitness = FMD_Fitness_Improved(params, x) modes = FMD(x, params); kurt = zeros(size(modes,1), 1); env_entropy = zeros(size(modes,1), 1); for i = 1 : size(modes,1) kurt(i) = kurtosis(modes(i,:)); env = abs(hilbert(modes(i,:))); env = env / sum(env); env_entropy(i) = -sum(env .* log(env + eps)); end % 加权：峭度越大越好，包络熵越小越好 fitness = mean(env_entropy ./ (kurt + eps)); end

这里kurtosis函数返回信号的峭度值，健康振动信号的峭度约等于3，滚动轴承早期故障冲击信号峭度通常在5以上，但过高的峭度也可能对应单个突发冲击而非周期性故障。env_entropy归一化后取值范围在4～8之间，与峭度量级差异不大，可以直接做算术运算。加权比值的思路是：让GWO同时抑制“分解出单个尖刺”的极端情况。

4.4 信号长度与分帧处理：优化前的必要预处理

FMD对信号长度很敏感。信号太短（少于1000点），滤波器的频带分辨率不够，优化结果波动大；信号太长，每次FMD计算时间过长导致整体无法忍受。建议在主程序里先截取一段稳态信号再进入GWO循环。

%% 截取定长信号，平衡计算量与信息量 seg_len = 8192; % 2的整数次幂，便于FFT计算 if length(x) > seg_len x_seg = x(end - seg_len + 1 : end); % 取尾部稳定段 else x_seg = x; end % 去均值与去趋势，避免直流分量干扰 x_seg = x_seg - mean(x_seg); x_seg = detrend(x_seg);

这段预处理的逻辑是：截取尾段是因为旋转机械启动或停机阶段信号不平稳，尾段往往是稳态工况。detrend去线性趋势项能避免数据缓慢漂移被FMD当成低频模态。需要注意，截取长度不一定是越短越好——要保证截取长度大于目标特征频率周期的10倍以上才能看到完整的故障冲击序列。

5. 常见的坑与排查路径：五条踩过的记录

5.1 现象：GWO收敛曲线不下降，适应度值始终在同一个水平

原因：搜索范围设置不合理，所有狼群的初始位置都落在平缓区域内；或者适应度函数本身返回值有问题。对比两个信号的特征：如果适应度值在优化前后几乎完全一样，那基本可以断定是适应度函数短路了——比如env_entropy计算时忘记做归一化，导致所有模态的熵值都被压缩到同一个量级。

解决：先用一组随机参数人工调用五次FMD，打印五组不同参数下的适应度，确认适应度有波动。如果没波动，检查FMD函数是否真的接收了params中的每个字段，或者是否在FMD内部硬编码了参数导致外部传参失效。

5.2 现象：最优解跑出来滤波器长度L是上限值，再大会出内存不足

原因：L越大，FMD构建的滤波器越复杂，计算量指数上升。当GWO发现L越大适应度越好时，会一直把L推向边界，最终得到一个物理上不可用的参数。

解决：遇到L顶到边界的情况，说明目标信号里存在非常窄的频带特征，或者适应度函数过度奖励模态稀疏度。先用频谱图确认目标频带宽度，把VarMax(1)改成基于实际带宽的合理上限。我一般会强制约束L <= round(fs / (2 * target_bw))，其中target_bw是以Hz为单位的预期最窄频带宽度。

5.3 现象：分解出的第一阶模态全是噪声，有效成分跑到第五阶之后

原因：FMD的模态输出顺序是按中心频率从高到低排列的，如果你的目标特征频率较低，对应的模态往往排在后面。新手经常只看前两阶模态，误以为分解失败。

解决：先查看各模态的中心频率，再定位目标模态。具体实现是对每个模态做FFT，求频谱峰值位置，然后对比目标特征频率。

%% 计算各模态中心频率 for i = 1 : size(modes,1) spec = abs(fft(modes(i,:))); [~, idx] = max(spec(1:length(spec)/2)); cf = (idx-1) * fs / length(spec); fprintf('IMF%d 中心频率 = %.2f Hz\n', i, cf); end

5.4 现象：同一组数据两次运行结果不一致，甚至模态数都变了

原因：GWO的种群初始化是随机生成随机数，没有固定随机种子，导致两次运行的最优解不同。如果第二次运行只是略有差异但模态结构一致，属于正常随机波动；如果模态结构和数量完全变了，说明优化收敛不稳定。

解决：在主程序开头固定随机种子：

rng(2024); % 固定随机种子，保证可复现

做实验对比或写论文时，固定随机种子是必须的。否则审稿人或同门复现你的结果时会得到完全不同的图，这会被质疑实验可信度。

5.5 现象：Matlab报错“数组索引超出范围”，定位在FMD_FLFFilter相关函数里

原因：滤波器长度L设置大于信号长度，或者smooth类函数的平滑跨度过大。FMD内部调用的零相位滤波函数对数据长度有硬性要求。

解决：在调用FMD之前显式检查参数合法性：

assert(params.L < length(x_seg), '滤波器长度L必须小于信号长度'); assert(params.lambda > 0, '步长lambda必须为正数');

另外，Matlab版本差异也会在这里暴露：2021a之后某些内置函数行为有调整，旧的源码在2023b上可能报错。遇到这种情况先看报错的具体函数名，再用edit <函数名>打开检查，通常把滤波实现从filtfilt换成filter加上手动边缘延拓就能解决。

6. 进阶技巧：定阶、模态合并与结果验证的三个习惯

FMD玩到后期，真正拉开差距的往往不是优化算法本身，而是分解后的处理。第一个进阶点是定阶。前面提到params.K不要设太大，但实际问题中并不知道真实分量数。我的习惯是做一次K=8的分解，然后看各模态间的相关系数矩阵，把相关系数大于0.9的相邻模态合并。因为FMD在K设置过大时，会把一个物理分量按频率切片拆成多段，这些相邻片段的波形高度相关。

合并模态在Matlab里就是简单的相加：

%% 合并高度相关的相邻模态 modes_merged = modes(1,:) + modes(2,:); % 若相关系数>0.9

第二个进阶是结果验证。GWO得到的适应度曲线收敛到0轴不下降，不代表你的故障特征一定能被提取出来。我会在分解完成后强制验证两步：第一步，对最优模态做包络谱分析，看故障特征频率及其倍频是否清晰谱线；第二步，把原始信号的能量与各模态能量之和做对比，正常情况下能量损失控制在5%以内才算分解完整。这两个验证代码量不大，但能让你在报告里理直气壮地说“故障特征提取有效”。

第三个习惯是批量回测。单条信号的优化结果没有统计意义。我会把同工况下的多段信号依次跑一遍GWO-FMD，统计最优参数分布区间的标准差。标准差小说明优化稳定、参数可复用；标准差大则说明信号非平稳性强，需要按工况分段重新寻优。这一趟流程走下来，比单纯贴一张分解图有说服力得多。

从第一次跑通到批量回测，我现在的习惯是每改一个参数就在代码里留一行注释记录改动原因，哪怕只是把步长从0.1改成0.08这种微调也会记一笔。因为GWO-FMD的参数耦合性很强，过两周回头看，你根本记不清当时那组好结果用的到底是哪个L值。把适应度曲线最终收敛值、最优参数向量、随机种子编号都记录下来，这个习惯帮我在复现实验时少走了很多弯路，希望也能帮到你。

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