Matlab实现序贯变分模态分解（SVMD）

        大家好，我是带我去滑雪！

       序贯变分模态分解(SVMD) 是一种信号处理和数据分析方法。它可以将复杂信号分解为一系列模态函数，每个模态函数代表信号中的特定频率分量。 SVMD 的主要目标是提取信号中的不同频率分量并将其重构为原始信号。SVMD的基本原理是通过变分模态分解的方式将信号分解为多个模态函数。在每个迭代步骤中，SVMD 通过最小化信号和模态函数之间的差异来更新模态函数。重复这个过程直到收敛。得到的模态函数可用于重建原始信号。

        SVMD 的另一个关键特征是连续分解。在每个迭代步骤中，SVMD 从信号中提取主频率分量并将其从信号中删除。这样，每次迭代步骤都会提取信号中的一个频率分量，直到提取完所有频率分量。这种逐次分解方法可以更好地捕获信号中的不同频率分量。SVMD 在信号处理和数据分析方面有着广泛的应用。它可用于去噪、特征提取、频谱分析等多个领域。通过将信号分解为模态函数，SVMD可以更好地理解和描述信号的频率特性。这对于信号处理和数据分析非常重要。SVMD的数据重构是将分解后的模态函数重新组合成原始信号的过程。通过各模态函数的加权相加即可得到重构信号。这个过程可以用来恢复原始信号的频率特性，并且可以根据需要进一步分析和处理。

         综上所述，逐次变分模态分解是一种有效的信号处理和数据分析方法。它可以将复杂信号分解为多个模态函数，并可以通过数据重构将它们重新组合成原始信号。 SVMD有着广泛的应用范围，对于理解和描述信号的频率特性非常有帮助。通过深入研究和应用SVMD，我们可以更好地处理和分析各类信号和数据。下面开始代码实战。

（1）SVMD实现

function [u,u_hat,omega]=svmd(signal,maxAlpha,tau,tol,stopc,init_omega)%% ------------ Part 1: Start initializingy = sgolayfilt(signal,8,25); %--filtering the input to estimate the noise
signoise=signal-y; %-estimating the noisesave_T = length(signal);
fs = 1/save_T;%______________________________________________________________________
%
%           Mirroring the signal and noise part to extend
%______________________________________________________________________
T = save_T;
f_mir=zeros(1,T/2);
f_mir_noise=zeros(1,T/2);
f_mir(1:T/2) = signal(T/2:-1:1);
f_mir_noise(1:T/2) = signoise(T/2:-1:1);
f_mir(T/2+1:3*T/2) = signal;
f_mir_noise(T/2+1:3*T/2) = signoise;
f_mir(3*T/2+1:2*T) = signal(T:-1:T/2+1);
f_mir_noise(3*T/2+1:2*T) = signoise(T:-1:T/2+1);f = f_mir;
fnoise=f_mir_noise;
%______________________________________________________________________
%______________________________________________________________________T = length(f);%------------- time domain (t -->> 0 to T)
t = (1:T)/T;udiff = tol+eps; %------ update stepomega_freqs = t-0.5-1/T;%------------- discretization of spectral domain%______________________________________________________________________
%
%     FFT of signal(and Hilbert transform concept=making it one-sided)
%______________________________________________________________________f_hat = fftshift((fft(f)));
f_hat_onesided = f_hat;
f_hat_onesided(1:T/2) =0;
f_hat_n = fftshift((fft(fnoise)));
f_hat_n_onesided = f_hat_n;
f_hat_n_onesided(1:T/2) =0;
%______________________________________________________________________
%______________________________________________________________________noisepe=norm(f_hat_n_onesided,2).^2;%------------- noise power estimationN = 300;%------------ Max. number of iterations to obtain each modeomega_L = zeros(N, 1);%----------- Initializing omega_dswitch nargincase 6if init_omega == 0omega_L(1) = 0;elseomega_L(1) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,1)));endotherwiseinit_omega = 0;omega_L(1) = 0;
endminAlpha=10; %------ the initial value of alpha
Alpha=minAlpha; %------ the initial value of alpha
alpha=zeros(1,1);
%----------- dual variables vector
lambda = zeros(N, length(omega_freqs));%---------- keeping changes of mode spectrum
u_hat_L = zeros(N, length(omega_freqs));n = 1; %------------------ main loop counterm=0;      %------ iteration counter for increasing alpha
SC2=0; % ------ main stopping criteria index
l=1;  %------ the initial number of modes
bf=0;  % ----- bit flag to increase alpha
BIC=zeros(1,1);  % ------- the initial value of Bayesian indexh_hat_Temp=zeros(2, length(omega_freqs));%-initialization of filter matrixu_hat_Temp=zeros(1,length(omega_freqs),1);%- matrix1 of modes 
u_hat_i=zeros(1, length(omega_freqs));%- matrix2 of modesn2=0; % ---- counter for initializing omega_Lpolm=zeros(2,1);  % ---- initializing Power of Last Mode indexomega_d_Temp=zeros(1,1);%-initialization of center frequencies vector1
sigerror=zeros(1,1);%initializing signal error index for stopping criteria
gamma=zeros(1,1);%----initializing gamma
normind=zeros(1,1);%% ---------------------- Part 2: Main loop for iterative updates
while (SC2~=1)while (Alpha(1,1)<(maxAlpha+1)) while ( udiff > tol &&  n < N ) %------------------ update uLu_hat_L(n+1,:)= (f_hat_onesided+...((Alpha(1,1).^2)*(omega_freqs - omega_L(n,1)).^4).*u_hat_L(n,:)+...lambda(n,:)/2)./(1+(Alpha(1,1).^2)*(omega_freqs - omega_L(n,1)).^4 ....*((1+(2*Alpha(1,1))*(omega_freqs - omega_L(n,1)).^2))+sum(h_hat_Temp));%------------------ update omega_Lomega_L(n+1,1) = (omega_freqs(T/2+1:T)*(abs(u_hat_L(n+1, T/2+1:T)).^2)')/sum(abs(u_hat_L(n+1,T/2+1:T,1)).^2);%------------------ update lambda (dual ascent)lambda(n+1,:) = lambda(n,:) + tau*(f_hat_onesided...-(u_hat_L(n+1,:) + (((Alpha(1,1).^2)*(omega_freqs - omega_L(n,1)).^4.....*(f_hat_onesided - u_hat_L(n+1,:)-sum(u_hat_i)+lambda(n,:)/2)-sum(u_hat_i))..../(1+(Alpha(1,1).^2)*(omega_freqs - omega_L(n,1)).^4 ))+...sum(u_hat_i)));udiff = eps;%------------------ 1st loop criterionudiff = udiff + (1/T*(u_hat_L(n+1,:)-u_hat_L(n,:))*conj((u_hat_L(n+1,:)-u_hat_L(n,:)))').../ (1/T*(u_hat_L(n,:))*conj((u_hat_L(n,:)))');udiff = abs(udiff);n = n+1;end%% ---- Part 3: Increasing Alpha to achieve a pure modeif abs(m-log(maxAlpha))> 1m=m+1;elsem=m+.05;bf=bf+1;endif  bf>=2Alpha=Alpha+1;endif   Alpha(1,1)<=(maxAlpha-1)  %exp(SC1)<=(maxAlpha)if (bf ==1)Alpha(1,1)=maxAlpha-1;elseAlpha(1,1)=exp(m);endomega_L=omega_L(n,1);% ------- Initializingudiff = tol+eps; % update steptemp_ud = u_hat_L(n,:);%keeping the last update of obtained moden = 1; % loop counterlambda = zeros(N, length(omega_freqs));u_hat_L = zeros(N, length(omega_freqs));u_hat_L(n,:)=temp_ud;endend%% Part 4: Saving the Modes and Center Frequenciesomega_L=omega_L(omega_L>0);u_hat_Temp(1,:,l)=u_hat_L(n,:);omega_d_Temp(l)=omega_L(n-1,1);alpha(1,l)=Alpha(1,1);Alpha(1,1)=minAlpha;bf=0;%------------------------------initializing omega_Lif init_omega >0ii=0;while (ii<1 && n2 < 300)omega_L = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,1)));checkp=abs(omega_d_Temp-omega_L);if (size(find(checkp<0.02),2)<=0) % it will continue if difference between previous vector of omega_d and the current random omega_plus is about 2Hzii=1;endn2=n2+1;endelseomega_L=0;endudiff = tol+eps; % update steplambda = zeros(N, length(omega_freqs));gamma(l)=1;h_hat_Temp(l,:)=gamma(l) ./((alpha(1,l)^2)*...(omega_freqs - omega_d_Temp(l)).^4);%---------keeping the last desired mode as one of the extracted modesu_hat_i(l,:)=u_hat_Temp(1,:,l);%%  Part 5: Stopping Criteria:if nargin >=5 % checking input of the functionswitch stopccase 1%-----------------In the Presence of Noiseif size(u_hat_i,1) == 1sigerror(l)=  norm((f_hat_onesided-(u_hat_i)),2)^2;elsesigerror(l)=  norm((f_hat_onesided-sum(u_hat_i)),2)^2;endif ( n2 >= 300 || sigerror(l) <= round(noisepe))SC2=1;endcase 2%-----------------Exact Reconstructionsum_u=sum(u_hat_Temp(1,:,:),3); % -- sum of current obtained modesnormind(l)=(1/T) *(norm(sum_u-f_hat_onesided).^2)..../((1/T) * norm(f_hat_onesided).^2);if( n2 >= 300 || normind(l) <.005 )SC2=1;endcase 3%------------------Bayesian Methodif size(u_hat_i,1) == 1sigerror(l)= norm((f_hat_onesided-(u_hat_i)),2)^2;elsesigerror(l)= norm((f_hat_onesided-sum(u_hat_i)),2)^2;endBIC(l)=2*T*log(sigerror(l))+(3*l)*log(2*T);if(l>1)if(BIC(l)>BIC(l-1))SC2=1;endendotherwise%------------------Power of the Last Modeif (l<2)polm(l)=norm((4*Alpha(1,1)*u_hat_i(l,:)./(1+2*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:)',2);polm_temp=polm(l);polm(l)=polm(l)./max(polm(l));elsepolm(l)=norm((4*Alpha(1,1)*u_hat_i(l,:)./(1+2*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:)',2);polm(l)=polm(l)./polm_temp;endif (l>1 && (abs(polm(l)-polm(l-1))<0.001) )SC2=1;endendelse%------------------Power of the Last Modeif (l<2)polm(l)=norm((4*Alpha(1,1)*u_hat_i(l,:)./(1+2*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:)',2);polm_temp=polm(l);polm(l)=polm(l)./max(polm(l));elsepolm(l)=norm((4*Alpha(1,1)*u_hat_i(l,:)./(1+2*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:)',2);polm(l)=polm(l)./polm_temp;endif (l>1 && (abs(polm(l)-polm(l-1))<tol) )SC2=1;endend%% Part 6: Resetting the counters and initializations  u_hat_L = zeros(N, length(omega_freqs));n = 1; % ----- reset the loop counterl=l+1; %---(number of obtained modes)+1m=0;n2=0;
end%% ------------------ Part 7: Signal Reconstructionomega = omega_d_Temp;
L=length(omega); %------number of modesu_hat = zeros(T, L);
u_hat((T/2+1):T,:) = squeeze(u_hat_Temp(1,(T/2+1):T,:));
u_hat((T/2+1):-1:2,:) = squeeze(conj(u_hat_Temp(1,(T/2+1):T,:)));
u_hat(1,:) = conj(u_hat(end,:));u = zeros(L,length(t));for l = 1:Lu(l,:)=real(ifft(ifftshift(u_hat(:,l))));
end[omega,indic]=sort(omega);
u=u(indic,:);
%---------- remove mirror part
u = u(:,T/4+1:3*T/4);%--------------- recompute spectrum
clear u_hat;
for l = 1:Lu_hat(:,l)=fftshift(fft(u(l,:)))';
endend

（2）SVMD绘图

function Huatu_svmd(emd_imf,signal,t,Fs)
if nargin <4
figure('Name','SVMD分解与各IMF分量时域图','Color','white');set(gcf, 'Position', [400 100 600 700]);
subplot(size(emd_imf,1)+1,1,1);
plot(t,signal,'k');grid on;
ylabel('\fontname{宋体}原始数据');title('\fontname{Times new roman}SVMD\fontname{宋体}分解');
set(gca,'XTick',[]);
for i = 2:size(emd_imf,1)+1subplot(size(emd_imf,1)+1,1,i);plot(t,emd_imf(i-1,:),'k');ylabel(['\fontname{Times new roman}IMF',num2str(i-1)]);if (i ~= size(emd_imf,1)+1)set(gca,'XTick',[]);endif (i == size(emd_imf,1)+1)ylabel('\fontname{Times new roman}RSE');xlabel('\fontname{Times new roman}Time/\it{s}');endgrid on;
end
else
figure('Name','SVMD分解与各IMF分量频谱对照图','Color','white');set(gcf, 'Position', [400 100 600 700]);
subplot(size(emd_imf,1)+1,2,1);
plot(t,signal,'k');grid on;
ylabel('\fontname{宋体}原始数据');
title('\fontname{Times new roman}SVMD\fontname{宋体}分解');
set(gca,'XTick',[]);
subplot(size(emd_imf,1)+1,2,2);
pFFT(signal,Fs);grid on;
title('\fontname{宋体}对应频谱');
set(gca,'XTick',[]);
for i = 2:size(emd_imf,1)+1subplot(size(emd_imf,1)+1,2,i*2-1);plot(t,emd_imf(i-1,:),'k');ylabel(['\fontname{Times new roman}IMF',num2str(i-1)]);if (i ~= size(emd_imf,1)+1)set(gca,'XTick',[]);endif (i == size(emd_imf,1)+1)ylabel('\fontname{Times new roman}RSE');xlabel('\fontname{Times new roman}Time/\it{s}');endgrid on;subplot(size(emd_imf,1)+1,2,i*2);pFFT(emd_imf(i-1,:),Fs);if (i ~= size(emd_imf,1)+1)set(gca,'XTick',[]);endif (i == size(emd_imf,1)+1)xlabel('\fontname{Times new roman}Frequency/\it{Hz}');endgrid on;
end
end

（3）测试

clc
clear
close all
SIG=importdata('NASA电容量.csv');
sig=SIG(2:161,2);%%想要分解哪一列就填几
%(1)导入时间数据来设置时间
t=SIG(2:161,2);
Fs=1/(t(2)-t(1));
%(2)设置采样率来设置时间
% N=length(sig);
% Fs=1000;%%采样频率自己设置
% t1=((0:N-1)*1/Fs)';
% SNR = 10;    
% sig = awgn(sig,SNR,'measured');   
figure('Name','原始信号');
% subplot(211);
plot(t,sig,'k');
title('\fontname{宋体}原始信号');
ylabel('\fontname{宋体}幅值');
xlabel('\fontname{Times new roman}Time/\it{s}');
% subplot(212);pFFT(sig,Fs)
% title('\fontname{宋体}频谱图');
% ylabel('\fontname{宋体}幅值');
% xlabel('\fontname{Times new roman}Frequency/\it{Hz}');
maxAlpha=1000; %compactness of mode
tau=0;%time-step of the dual ascent
tol=1e-6; %tolerance of convergence criterion;
stopc=4;%the type of stopping criteria
[svmd_imf,uhat,omega]=svmd(sig,maxAlpha,tau,tol,stopc);
Huatu_svmd(svmd_imf,sig,t);
svmd_imf=svmd_imf';

需要数据集的家人们可以去百度网盘（永久有效）获取：

链接：https://pan.baidu.com/s/173deLlgLYUz789M3KHYw-Q?pwd=0ly6
提取码：2138 

---

更多优质内容持续发布中，请移步主页查看。

博主的WeChat:TCB1736732074

   点赞+关注,下次不迷路！