重头戏来了。
在以往的应用经验里，VMD方法在众多模态分解方法中可以说是非常好的。从催更力度上看，这个方法也是格外受关注。笔者决定加快进度快一些写完这个方法，十月份了有些同学要开始做毕设，希望这篇文能帮上忙。
1. VMD（变分模态分解）的概念
VMD（Variational Mode Decomposition）即变分模态分解，与2014年由Dragomiretskiy[1]等人提出，虽然它也叫模态分解，但是和之前介绍过的EMD、EEMD、CEEMD、CEEMDAN这些模态分解方法在原理上有本质区别。VMD的整体思路为：
该方法假设任何的信号都是由一系列具有特定中心频率、有限带宽的子信号组成（即IMF）。
以经典维纳滤波为基础，通过对变分问题进行求解，得到 中心频率与带宽限制，找到 各中心频率在频域中对应的有效成分，得到模态函数。其模型构建涉及到维纳滤波、希尔伯特变换和解析信号等知识点。 [2]
VMD 的分解过程就是变分问题的求解过程，其算法主要包括 变分问题的构造和 变分问题的求解。 [3]VMD的求解过程主要包含两点约束：（1）要求每个模态分量中心频率的带宽之和最小；（2）所有的模态分量之和等于原始信号。 [4]
1.1 关于内涵模态分量
不同于黄锷先生提出的内涵模态分量（IMF）概念，VMD算法重新定义了约束条件更为严格的有限带宽的本证模态函数，该内涵模态分量被定义为调幅调频的分量模态函数，数学表达式为：
为信号
的包络幅值，
为瞬时相位。需要注意该函数表征的分量也是同样满足EMD的约束条件的。
1.2 变分问题构造
由于变分问题涉及到泛函的相关知识，推导过程比较复杂，这里只介绍思路，对于详细过程有兴趣的同学可以直接看方法提出者的论文原文[1]。
所谓变分问题，就是求泛函的极值。在VMD中，泛函指的是VMD约束变分模型，而要求的极值，就是“每个模态分量中心频率的带宽之和最小”。
过去常遇到的是求函数极值，但有时我们需要对自变量也是函数的特殊函数求极值。
这种特殊函数即“函数的函数”，称为泛函，求泛函的极值问题称为变分问题。
构造出来的VMD约束变分模型是这样的[1]：
*上式表示的，就是求取“每个模态分量中心频率的带宽之和最小”时的模态函数
和中心频率
。
1.3 变分问题求解
上式的求解过程应用到了比较多的数学工具，包括二次惩罚项、拉格朗日算子以及增广拉格朗日函数等等，详细过程[1]比较复杂，同样不展开说了。
2.VMD的显著特点
2.1 VMD思路
VMD方法是一种非递归变分模式的信号分解方法，其整体框架为一个变分问题。
2.2 VMD的独特优势
（1）可以指定想要得到的模态数。
（2）通过VMD方法分解出来的IMF都具有独立的中心频率，并且在频域上表现出稀疏性的特征，具备稀疏研究的特质。
（3）在对IMF求解过程中，通过镜像延拓的方式避免了类似EMD分解中出现的端点效应。
（4）有效避免模态混叠（K值选取合适的情况下）。
2.3 VMD的重要参数——模态数K
在使用EMD方法时，有很多同学就经常问要怎样指定分解出的IMF数。VMD方法就可以很好地满足这个要求，因为VMD的一个明显特点，就是信号在做VMD分解前，需要先设定模态分量的个数K。但是成也萧何败也萧何，某些应用场合下可指定IMF数属于优点，但是对于某些不预知信号隐含模态数的场景，怎样设置这个K值反而会让人挠头。
若设定的K小于待分解信号中有用成分的个数（欠分解），会造成分解不充分，导致模态混叠；若设定的K值大于待分解信号中有用成分的个数（过分解），就导致产生一些没有用的虚假分量。因此，K值的确定对于VMD就非常重要。
至于如何确定K值，可以从以下几个方面考虑：
（1）某些场景下可以预知模态数目，比如对于某齿轮箱振动数据，通过分析可知信号数据中主要包含齿轮的啮合频率和主轴转动频率，那么K可以设置为2。
（2）可以通过K值从小到大逐个尝试，并结合分解结果分析确定K值：随着K值增大各主要频率段数据能分布到不同IMF分量中，并且没有产生虚假分量，此时K值就是比较合适的。这个方法虽然比较高效且靠谱，但是不便于写到论文中，这就需要用到方法（3）了。
（3）可以使用一些辅助算法，比如使用峭度最大原理[5]、能量差值原则[6]等，或者结合一些寻优算法对K值（也可以同时对α）寻优。
2.4 VMD的重要参数——惩罚系数α
VMD 分解的过程中，预设的 K 值决定着 IMF 分量的个数，惩罚系数 α 决定着 IMF 分量的带宽。惩罚系数越小， 各 IMF 分量的带宽越大，过大的带宽会使得某些分量包含其他分量信号；α值越大，各IMF分量的带宽越小，过小的带宽是使得被分解的信号中某些信号丢失。该系数常见取值范围为1000~3000。
2.5 VMD的另一个参数——收敛容差tol
是优化的停止准则之一，即在连续两次迭代中，当向 IMF 收敛的绝对平均平方改进小于 tol时，优化停止。通常可以取 1e-6~5e-6。
3. VMD的应用案例
vmd函数在MATLAB2020a版本中内置到官方库了，不过按照“类EMD”系列的代码的统一风格，笔者重新进行了封装，一来使该方法可以应用于MATLAB2020以前版本，二来增加了绘制IMF分量与频谱对照的绘图函数，同时做了一些其他的适应性调整。
下边使用一个例子来展示VMD方法的突出效果：
首先生成一组测试信号，由二次趋势项、啁啾信号和分段余弦共同组成，采样频率为1000Hz。
%% 1.生产仿真信号
fs = 1e3;
t = 0:1/fs:1-1/fs;
x = 6*t.^2 + cos(4*pi*t+10*pi*t.^2) + ...
[cos(60*pi*(t(t<=0.5))) cos(100*pi*(t(t>0.5)-10*pi))];
figure('color','w')
tiledlayout('flow')
nexttile
plot(t,[zeros(1,length(t)/2) cos(100*pi*(t(length(t)/2+1:end))-10*pi)])
xlabel('Time (s)')
ylabel('Cosine')
nexttile
plot(t,[cos(60*pi*(t(1:length(t)/2))) zeros(1,length(t)/2)])
xlabel('Time (s)')
ylabel('Cosine')
nexttile
plot(t,cos(4*pi*t+10*pi*t.^2))
xlabel('Time (s)')
ylabel('Chirp')
nexttile
plot(t,6*t.^2)
xlabel('Time (s)')
ylabel('Quadratic trend')
nexttile(5,[1 2])
plot(t,x)
xlabel('Time (s)')
ylabel('Signal')
四个信号分量及最终合成的Signal信号
这种程度的复合信号可以说是比较复杂的了，让我们看下对于该信号的VMD分解效果：
%% VMD分解参数设置
alpha=2000; % alpha - 惩罚因子
tol=1e-6; % tol - 收敛容差，是优化的停止准则之一，可以取 1e-6~5e-6
K=4; % K - 指定分解模态数
[imf,CenFs] = pVMD(x,fs, alpha, K, tol); %调用函数实现分解与画图
Signal信号的VMD分解结果
上图是对于Signal信号的VMD分解结果，可以说是近乎完美了。
作为对比，EMD分解结果是下图这样的，就十分差强人意了。
Signal信号的EMD分解结果
在EEMD、CEEMD等改进方法中得到的结果也同样不理想，也正因此VMD方法可以说是笔者在解决工程问题时最常用的分解方法。
4.关于封装函数
上边案例中用到的函数“pVMD”就是笔者封装好的vmd画图程序，只需要输入待分解信号、采样频率、α值、K值、tol值，就可以得到信号的分解结果IMF和中心频率CenFs。这个函数的介绍如下：
[imf,CenFs]=pVMD(x,fs, alpha, K, tol);
% function [imf,CenFs] = pVMD(y,FsOrT, alpha, K, tol)
% 画信号VMD分解图
% 输入：
% data：待分解的数据(一维)
% FsOrT：采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% alpha - 惩罚因子
% K - 指定分解模态数
% tol - 收敛容差，是优化的停止准则之一，可以取 1e-6~5e-6
% 输出：
% imf：内涵模态分量，统一为n*m格式，其中n为模态数，m为数据点数。例如 imf(1,:)即IMF1，imf(end,:)即为残差
% 注意，为了与其他“类EMD”方法分解出来的imf分量保持一致，本程序内将imf排序进行了翻转
% 即保证imf排列是从高频向低频排列
% CenFs：即CentralFrequencies，各imf分量的中心频率
% 注意：在使用该代码之前，请务必安装工具箱：http://www.khscience.cn/docs/index.php/2020/04/09/1/按照惯例，还有一个绘制VMD分解图及其频谱图的函数：
[imf,CenFs] = pVMDandFFT(x,fs, alpha, K, tol);
% function [imf,CenFs] = pVMDandFFT(y,FsOrT, alpha, K, tol)
% 画信号VMD分解与各IMF分量频谱对照图
% 输入：
% y：待分解的数据(一维)
% FsOrT： 采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% alpha - 惩罚因子
% K - 指定分解模态数
% tol - 收敛容差，是优化的停止准则之一，可以取 1e-6~5e-6
% 输出：
% imf：内涵模态分量，统一为n*m格式，其中n为模态数，m为数据点数。例如 imf(1,:)即IMF1，imf(end,:)即为残差
% 注意，为了与其他“类EMD”方法分解出来的imf分量保持一致，本程序内将imf排序进行了翻转
% 即保证imf排列是从高频向低频排列
% CenFs：即CentralFrequencies，各imf分量的中心频率
% 例1：（FsOrT为采样频率）
% fs = 100;
% t = 1/fs:1/fs:1;
% y = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pVMDandFFT(y,fs,2000,2,1e-7);
% 例2：（FsOrT为时间向量，需要注意此时FsOrT的长度要与y相同）
% t = 0:0.01:1;
% y = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pVMDandFFT(y,t,2000,2,1e-7);
% 注意：在使用该代码之前，请务必安装工具箱：http://www.khscience.cn/docs/index.php/2020/04/09/1/对于上边的案例，使用这个函数就可以得到各个imf分量的频谱图了：
Signal信号的VMD分解结果及其频谱图
