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分解展示。代码如下：

clc
clear
fs=12000;%采样频率
Ts=1/fs;%采样周期
L=2000;%采样点数
t=(0:L-1)*Ts;%时间序列
%----------------导入内圈故障的数据-----------------------------------------
load 105.mat
X = 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 on
xlabel('Time Values')
ylabel('Mode Number')
zlabel('Mode Amplitude')

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

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

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

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

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

clc
clear
fs=12000;%采样频率
Ts=1/fs;%采样周期
L=2000;%采样点数
t=(0:L-1)*Ts;%时间序列
%----------------导入内圈故障的数据-----------------------------------------
load 105.mat
X = 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)")
end
title(t1,"Variational Mode Decomposition")

结果如下：

VMD去噪

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

cleanX = sum(imf(:,2:8),2);
figure
plot(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为例：

clear
clc
load 105.mat
L=2000; %采样点数
y = X105_DE_time(1:L); %这里可以选取DE(驱动端加速度)、FE(风扇端加速度)、BA(基座加速度)，直接更改变量名，挑选一种即可。
n=40;
%=============================
%==== 奇异值分解 ====
for i=1:L/2+1
t=i:i+L/2-1;
for j=1:L/2
A(j,i)=y(t(j)); %把一维信号重构为矩阵做奇异值分解
end
end
[U,S,V] = svd(A);
%=============================
%==== 重构信号 =====
X=zeros(size(A));
for i=1:n  %选取前n个大奇异值
X=X+U(:,i)*S(i,i)*V(:,i)';
end
JG=zeros(1,L);
for i=1:L
a=0;m=0;
for j1=1:L/2
for j2=1:L/2+1
if i+1==j1+j2
a=a+X(j1,j2); %把矩阵重构回一维信号
m=m+1;
end
end
end
JG(i)=a/m;
end
% JG(L/2+1:end)=X(:,end);
figure
fs=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