matlab实现VMD去噪、SVD去噪，源代码详解

为了更好的利用MATLAB自带的vmd、svd函数，本期作者将详细讲解一下MATLAB自带的这两个分解函数如何使用，以及如何画漂亮的模态分解图。

VMD函数用法详解

首先给出官方vmd函数的调用格式。

[imf,residual,info] = vmd(x)

函数的输入：

这里的x是待分解的信号，一行或者一列都可以！

函数的输出：

第一个参数位置--imf：vmd分解得到的IMF分量。

第二个参数位置--residual：残余分量，残差表示原始信号x中未被vmd分解的部分。

第三个参数位置--info，是一个结构体，包含以下信息：

ExitFlag——终止标志。值0表示算法在达到最大迭代次数时停止。值为1表示算法在满足绝对和相对公差时停止。

CentralFrequencies——IMF的中心频率。

NumIterations——迭代总数。

AbsoluteImprovement——在最后两次迭代之间，IMF收敛的均方绝对改进。

RelativeImprovement——在最后两次迭代之间，IMF收敛的平均相对改进。

LagrangeMultiplier ——上一次迭代时的频域拉格朗日乘数。

示例：

以西储大学轴承故障的105.mat为例，进行vmd分解展示。代码如下：

clcclearfs=12000;%采样频率Ts=1/fs;%采样周期L=2000;%采样点数t=(0:L-1)*Ts;%时间序列%----------------导入内圈故障的数据-----------------------------------------load 105.matX = X105_DE_time(1:L); %这里可以选取DE(驱动端加速度)、FE(风扇端加速度)、BA(基座加速度)，直接更改变量名，挑选一种即可。[imf,residual,info] = vmd(X);figure(1);[p,q] = ndgrid(t,1:size(imf,2));plot3(p,q,imf)grid onxlabel('Time Values')ylabel('Mode Number')zlabel('Mode Amplitude')

选取105.mat的2000个采样点，进行分解，结果图如下所示：

还可以指定vmd的模态分解数和惩罚因子：

[imf,residual,info] = vmd(X,'NumIMF',6,'PenaltyFactor',2500,')

可以看到，模态分解数变成了6个。

当分解模态数较多的时候，还可以这样画图：

clcclearfs=12000;%采样频率Ts=1/fs;%采样周期L=2000;%采样点数t=(0:L-1)*Ts;%时间序列%----------------导入内圈故障的数据-----------------------------------------load 105.matX = X105_DE_time(1:L); %这里可以选取DE(驱动端加速度)、FE(风扇端加速度)、BA(基座加速度)，直接更改变量名，挑选一种即可。[imf,residual] = vmd(X,'NumIMF',9);t1 = tiledlayout(3,3,TileSpacing="compact",Padding="compact");for n = 1:9    ax(n) = nexttile(t1);    plot(t,imf(:,n)')    xlim([t(1) t(end)])    txt = ["IMF",num2str(n)];    title(txt)    xlabel("Time (s)")endtitle(t1,"Variational Mode Decomposition")

结果如下：

VMD去噪

对分解的信号去除残余分量，然后进行重构，即为去噪，代码如下：

cleanX = sum(imf(:,2:8),2);figureplot(t,X,t,cleanX)legend("Original X","Clean X")xlabel("Time (s)")ylabel("Signal")

这里选用IMF2-IMF8进行相加，丢弃了IMF1和IMF9。为什么要这么做呢，以下MATLAB官方文档给出的解释：

第一种模式包含最多的噪声，第二种模式是以重要特征的频率振荡。通过对除第一个和最后一个VMD模式之外的所有模式求和来构建干净的X信号，从而丢弃低频基线振荡和大部分高频噪声。

结果如下：

105.mat下载地址：https://pan.baidu.com/s/1UX0dgLrW6SpjPFhML-OQeQ?pwd=puz2 

---

SVD去噪

SVD就是常见的奇异值分解，也算是一个古老的分解方法了。MATLAB官方关于svd的调用方法如下：

[U,S,V]=svd(A)

参数说明：A：要进行SVD分解的矩阵，且A不必为方阵。

U、S、V为返回的SVD分解的结果：

• U：A的奇异向量组成的矩阵，AA’的正交单位特征向量组成U；

• S：对角矩阵，对角元素是A的奇异值，非负且按降序排列；

• V：A的奇异向量组成的矩阵，A’A的正交单位特征向量组成V。

SVD去噪：同样以105.mat为例：

clearclcload 105.matL=2000; %采样点数y = X105_DE_time(1:L); %这里可以选取DE(驱动端加速度)、FE(风扇端加速度)、BA(基座加速度)，直接更改变量名，挑选一种即可。n=40;%=============================%==== 奇异值分解 ====for i=1:L/2+1t=i:i+L/2-1;for j=1:L/2A(j,i)=y(t(j)); %把一维信号重构为矩阵做奇异值分解endend[U,S,V] = svd(A);%=============================%==== 重构信号 =====X=zeros(size(A));for i=1:n  %选取前n个大奇异值X=X+U(:,i)*S(i,i)*V(:,i)';endJG=zeros(1,L);for i=1:La=0;m=0;for j1=1:L/2for j2=1:L/2+1if i+1==j1+j2a=a+X(j1,j2); %把矩阵重构回一维信号m=m+1;endendendJG(i)=a/m;end% JG(L/2+1:end)=X(:,end);figurefs=12000;%采样频率Ts=1/fs;%采样周期L=2000;%采样点数t=(0:L-1)*Ts;%时间序列plot(t,y,t,JG)legend("Original X","Clean X")xlabel("Time (s)")ylabel("Signal")title('奇异值降噪信号','FontSize',16)%set(gca,’linewidth’,1.5);set(gcf,'color','w')

105.mat下载地址：https://pan.baidu.com/s/1UX0dgLrW6SpjPFhML-OQeQ?pwd=puz2 

---

获取更多代码：

或者复制链接跳转：https://docs.qq.com/sheet/DU3NjYkF5TWdFUnpu