简介VMD变分模态分解算法是信号处理与故障诊断中的常用分解工具这份资源面向使用MATLAB开展模态分解、中心频率分析的研究者与工程师重点解决VMD分解阶数难以确定的问题——通过观察各模态中心频率的相近程度可据此选择合理的分解层数从而获得具有明确物理含义的模态分量。压缩包共包含14个文件整体大小约1.33MB其中4个m文件为主要实现脚本含主程序VMD.m与多个测试用例3个fig与3个jpg为结果图与对比图可直接查看分解效果另含3个asv自动备份文件及1个mat数据样例便于运行验证无需额外准备数据。资源附带的图片与fig文件直观展示了VMD分解前后的波形与频谱适合用于复现算法流程、调整参数或撰写相关实验报告。自发布以来已有622人学习参考对入门VMD算法和掌握按中心频率确定最佳阶数的思路具有实用价值。1. 中心频率相近原则VMD最佳阶数判据的切入点标题里的“最佳阶数”原题写作“阶说”说的是 VMD 算法中的模态分量个数 K。VMD 最折磨人的参数不是惩罚系数 alpha而是这个 KK 取小了欠分解两个真实分量挤进同一个模态带里K 取大了过分解一个真实分量被劈成两个形状怪异的伪模态。中心频率相近原则的切入点是VMD 每个模态分量在收敛后都会得到一个可确定的中心频率当某两个中心频率明显贴近说明频域通带已经重叠当前 K 已经超出了信号实际承载的模态数反之相邻中心频率间隔清楚这个 K 才配得上“最佳”。这篇文章我会拿 MATLAB 把这条路径完整走通先讲 VMD 给出可确定模态分量的原理再给中心频率扫描与相近量化方法然后写自动定阶函数最后用合成信号验证边界。适合做振动故障诊断、地震信号处理和金融时序分解的工程师。2. VMD模型与MATLAB迭代模态分量的可确定性来自哪2.1 变分约束每个模态分量都是一条窄带信号VMD 先把原始信号 f 看成 K 个模态分量 u_k 的叠加再构造一个带约束的变分问题目标函数是所有模态分量的估计带宽之和最小等式约束是 K 个模态重构后等于原信号。写成优化问题就是min ∑_k ‖∂_t[(δ(t)j/πt) * u_k(t)]‖₂²s.t. ∑_k u_k f。这个形式说明 VMD 不是直接做滤波而是把“模态分量带宽要窄”和“重构误差要小”两个目标同时放进一个优化问题。求解时用 ADMM 交替方向乘子法把问题拆成 u_k、中心频率 ω_k、拉格朗日乘子 λ 三个子问题循环更新。模态分量 u_k 的更新在频域等价于一个中心频率可变的带通滤波所以每次更新完 u_k中心频率也顺带更新一次最终两者同时收敛。2.1.1 中心频率为什么是“可确定”的量中心频率 ω_k 在 VMD 里的定义是模态频谱的“重心”。因为 u_k 是带限信号频谱能量集中在一个频带内重心位置就决定这个带通滤波器落在哪。内侧问题里它跟着频谱能量走每次迭代用的是更新后的功率谱所以只要 K、alpha、init 固定同一段信号的 ω_k 收敛值几乎不变这就构成了“可确定模态分量”的直接依据。2.2 ADMM 迭代与 vmd 调用的 MATLAB 实现先给一个能跑出数值的核心骨架帮助理解内建函数到底在迭代什么% 说明用内核骨架只展开 ADMM 频域更新部分 % f_hat 是输入信号 FFTomega_grid 是归一化角频率列向量K 是模态数 N length(f_hat); omega_grid (0:N-1) / N * 2 * pi; % 列向量单位是 rad/sample omega_k linspace(0, pi, K); % 中心频率初始值 for n 1:200 % 一般 tol 会先触发退出 for k 1:K others sum(u_hat, 2) - u_hat(:, k); % 其余模态的频域和 numerator f_hat - others lambda_hat / 2; u_hat(:, k) numerator ./ ... (1 2*alpha * (omega_grid - omega_k(k)).^2); % 频域维纳滤波 omega_k(k) sum(abs(u_hat(:,k)).^2 .* omega_grid) ... / sum(abs(u_hat(:,k)).^2); % 频谱重心 end lambda_hat lambda_hat tau * (f_hat - sum(u_hat, 2)); if norm(f_hat - sum(u_hat, 2)) / norm(f_hat) tol break; end end逻辑说明u_k 的更新是用当前残差除以一个与 alpha 和中心频率距离有关的惩罚项本质就是把残差信号往 ω_k 附近收紧中心频率 ω_k 再按新功率谱重心移动。两行代码交替执行模态分量和中心频率互相牵引着收敛。拉格朗日乘子项把重构误差压到 tol 以下保证分解不过度偏离原信号。参数说明alpha 越大滤波器越窄模态分量带宽越小中心频率区分度越高tau 取 0 是噪声鲁棒模式对重构精度略松tol 控制收敛精度别的参数不变时只影响中心频率的小数位和耗时。实际工程里你不需要自己写这套循环MATLAB 从 R2019b 开始内置vmd(x, NumIMF, K, Alpha, alpha)输出信息info.CentrFreq直接给中心频率File Exchange 上常见的旧版 vmd.m 则返回角频率向量 omega。2.3 alpha、tau、tol 如何影响中心频率分布下面这张表是我做中心频率相近判据时的一组稳定起点参数常用值对中心频率和模态分量的影响alpha2000带宽窄ω_k 稳定太大时低频分量容易被吞掉tau0噪声鲁棒模式ω_k 更接近真实谱峰tol1e-6决定 ω_k 末位抖动幅度一般不需要更小init1初始中心频率均匀铺开避免局部收敛DC0置 1 时让第一个模态固定在 0 频附近慢趋势分量不会被劈开注意alpha 和 K 是耦合的。alpha 设太小会让本应分开的两个真实频率互相粘连中心频率显示成“相近”但这不是过分解而是带宽太宽。我一般把 alpha 固定为 2000把 K 当作唯一扫描变量。3. 中心频率扫描按相近原则手动确定最佳K3.1 时域波形不能替代中心频率判据K 偏大时伪模态在时域里往往看不出异常IMF 波形依然平滑但频谱上两个模态分量的通带已经叠到一起。一个通带里出现两个中心频率就是过分解的最直接证据。所以中心频率扫描才比直接看 IMF 更可靠这也是“中心频率相近原则”在实践里的入口。3.2 一个 K3..8 的扫描脚本% 设 x 为 1×N 行向量dt 为采样间隔fs 1/dt dt 1 / 1000; % 假设采样率 1000 Hz K_list 3:8; centers cell(length(K_list), 1); alpha 2000; for i 1:length(K_list) K K_list(i); [~, ~, info] vmd(x, NumIMF, K, Alpha, alpha); centers{i} info.CentrFreq; % 内建函数直接给出中心频率单位 Hz end figure; hold on; for i 1:length(K_list) line_pos (length(K_list) - i 1) * ones(size(centers{i})); stem(centers{i}, line_pos, filled, MarkerSize, 4); end set(gca, YTick, 1:length(K_list), YTickLabel, string(K_list)); xlabel(中心频率 / Hz); ylabel(扫描档位 K);逻辑说明扫描 K 是观察 dMin 曲线的第一步。每个 K 调用一次 vmd把内建返回值info.CentrFreq拿出来画成点。当某一档 K 的相邻两个点几乎重合或者上下两档之间突然插入一对靠得很近的点就要马上检查是不是发生了伪分解。参数说明NumIMF是内建 vmd 的模态数参数旧版 File Exchange 的 vmd.m 通常用位置参数传 K返回的 omega 是角频率需要换算成 Hzomega / (2*pi*dt)。扫描宽度建议取你预估模态数的两倍最低从 3 开始因为 K1 和 K2 的判据信息太少。3.3 相近程度量化用相对差还是绝对差判断“相近”不能直接用绝对频率差。低频段 0.1 Hz 和 0.15 Hz 相差 0.05 Hz看上去很小但相对间隔已经有 40%不该被当成一对伪模态高频段 38 Hz 和 40.5 Hz 绝对值更大相对间隔只有 6.4%反而是高度可疑的频带头碰撞。我统一用相邻中心频率的相对差d_i (ω_{i1} − ω_i) / mean(ω_{i1}, ω_i)把 d_i 里的最小值记录下来叫 dMin。当 dMin 小于 0.1 时认为出现了中心频率相近的过分解信号。下面是一个含 0.3 Hz 慢趋势、10 Hz 正弦和 40 Hz 调幅信号的扫描实例扫描档位中心频率排列判读K30.33 / 10.02 / 39.87 Hz三个模态分量间隔清晰dMin 约 1.2最佳K40.33 / 9.98 / 38.15 / 40.72 Hz40 Hz 峰被劈成两个dMin≈0.065过分解手动判断的完整过程是先跑 K3 到 K8画出每一档的中心频率计算每档 dMin找到 dMin 第一次掉到 0.1 以下的档位把它的上一档取为最佳阶数。4. 自动定阶的 MATLAB 实现与阈值优化4.1 固定阈值失灵的两个原因直接把 dMin 0.1 当硬阈值会出问题。第一dMin 和信号自身频率尺度有关低频信号本来就存在小绝对差第二alpha 改变带宽同一对频率距离在 alpha500 时是正常间隔在 alpha5000 时就变成通带重叠。更稳定的做法是同时看两件事本档 dMin 是否足够小以及它相对上一档是否发生骤降。只有两个条件同时满足才判定刚发生了一次伪分裂。4.2 基于 dMin 曲线的自动定阶函数function bestK autoVmdK(x, alpha, Kmax, dThresh) if nargin 3, Kmax 8; end if nargin 4, dThresh 0.1; end dMin nan(1, Kmax); for K 2:Kmax [~, ~, info] vmd(x, NumIMF, K, Alpha, alpha); om sort(info.CentrFreq(:)); d diff(om) ./ ((om(1:end-1) om(2:end)) / 2); dMin(K) min(d); end bestK 1; for K 3:Kmax if dMin(K) dThresh * 0.8 dMin(K) dMin(K-1) * 0.7 bestK K - 1; % 本档伪分裂上一档就是最佳 break; end end if bestK 1 bestK Kmax; % 全部档位都未触发相近警告 end fprintf(dMin: ); fprintf(%.4f , dMin(2:end)); fprintf(\n); end逻辑说明函数先跑 K2 到 Kmax 的 VMD每档取排序后的中心频率计算相邻相对差记录最小值 dMin。判定循环里有硬阈值和趋势两个条件本档 dMin 低于 0.08dThresh 的 0.8 倍并且比上一档下降至少 30%。两个条件同时为真说明新增的模态把一个原本清晰的峰值劈开了此时返回上一档 K−1。参数说明dThresh默认 0.1对应中心频率间距低于 10% 时视为可疑强噪声场景建议收紧到 0.08而 alpha 很大的窄带场景可以放宽到 0.15。如果你用的是旧版 vmd.m把info.CentrFreq替换成omega / (2*pi*dt)即可其他逻辑不用动。4.3 中心频率不稳定的处理重复运行取中位数“可确定”不等于每一次都跳进同一个局部极小值。遇到强干扰或初始频率靠近真实峰时中心频率曲线会抖动。判断稳定性可以这样对同一个 K 重复调用 vmd 三次看最大偏差是否超过 0.5 Hz。稳妥的自动化做法是每个 K 跑 5 次取中心频率的中位数再算 dMinomAll zeros(5, K); for rep 1:5 [~, ~, info] vmd(x, NumIMF, K, Alpha, alpha); omAll(rep, :) sort(info.CentrFreq); end om median(omAll, 1); % 中位数比单次结果抗抖动 d diff(om) ./ ((om(1:end-1) om(2:end)) / 2);逻辑说明每次 vmd 的内部初始化是随机的重复运行得到一组中心频率。中位数对极端落点不敏感比均值更稳。代价是计算量变成 5 倍只在信噪比较低时才值得开。参数说明重复次数 5 是个经验值改成 3 会损失一点稳定性改成 10 则收益递减。对于采样点数超过 10 万的信号建议先用降采样或分段分解跑出 K 的大致范围再全段定阶。5. 验证技巧用合成信号校准你的最佳阶数判据5.1 先跑已知 K 的信号验证自动定阶结果fs 2560; dt 1/fs; t (0:fs*2-1) * dt; % 真实模态数 30.4 Hz 慢趋势 10 Hz 正弦 40 Hz 调幅 x 0.5*cos(2*pi*0.4*t) cos(2*pi*10*t) ... (1 0.3*cos(2*pi*3*t)).*cos(2*pi*40*t); bestK autoVmdK(x, 2000, 8, 0.1); % 期望输出 bestK 3再叠加幅度 0.2 的白噪声重复调用同一个函数。你会发现 dMin 从第三档开始不再单调下降而是出现上下抖动。这个抖动本身就是信号已经被噪声污染的信号VMD 会把噪声拆到高频小分量里产生大量低频假中心频率。5.2 两种边界条件的收尾手法第一类边界是信号里本来就有两个很近的真实分量比如 38 Hz 和 41 Hz相对差约为 7.6%会被 dThresh0.1 误判成过分解。对策是先把 alpha 调大到 5000让滤波器变窄到能区分这两个峰再看 dMin 是否仍有骤降。第二类边界是强噪声场景下 VMD 在高频端造出一串伪模态此时建议把扫描范围从高频段截断只保留占总能量 99% 的频带再重新跑 autoVmdK。最终我一般会留一个朴素的交叉验证取自动定阶给出的 K分解后把每个模态分量和原始信号的几个谱峰做相关核对最强相关的模态数等于 K才算彻底落定。本文还有配套的精品资源点击获取