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

带我去滑雪

已于 2024-03-25 17:49:26 修改

阅读量1.7k

点赞数 11

分类专栏：模态分解文章标签： matlab 开发语言

于 2024-03-25 17:49:15 首次发布

本文链接：https://blog.csdn.net/qq_45856698/article/details/137019612

版权

模态分解专栏收录该内容

2 篇文章

订阅专栏

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

序贯变分模态分解(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 initializing

y = sgolayfilt(signal,8,25); %--filtering the input to estimate the noise
signoise=signal-y; %-estimating the noise

save_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 step

omega_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 estimation


N = 300;%------------ Max. number of iterations to obtain each mode


omega_L = zeros(N, 1);%----------- Initializing omega_d

switch nargin
    case 6
        if init_omega == 0
            omega_L(1) = 0;
        else
            omega_L(1) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,1)));
        end
    otherwise
        init_omega = 0;
        omega_L(1) = 0;
end

minAlpha=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 counter

m=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 index

h_hat_Temp=zeros(2, length(omega_freqs));%-initialization of filter matrix

u_hat_Temp=zeros(1,length(omega_freqs),1);%- matrix1 of modes 
u_hat_i=zeros(1, length(omega_freqs));%- matrix2 of modes

n2=0; % ---- counter for initializing omega_L


polm=zeros(2,1);  % ---- initializing Power of Last Mode index

omega_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 uL
            u_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_L
            omega_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 criterion
            udiff = 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 mode
        
        if abs(m-log(maxAlpha))> 1
            m=m+1;
        else
            m=m+.05;
            bf=bf+1;
        end
        if  bf>=2
            Alpha=Alpha+1;
        end
        if   Alpha(1,1)<=(maxAlpha-1)  %exp(SC1)<=(maxAlpha)
            
            if (bf ==1)
                Alpha(1,1)=maxAlpha-1;
            else
                Alpha(1,1)=exp(m);
            end
            
            omega_L=omega_L(n,1);
            
            % ------- Initializing
            udiff = tol+eps; % update step
            
            temp_ud = u_hat_L(n,:);%keeping the last update of obtained mode
            
            n = 1; % loop counter
            
            lambda = zeros(N, length(omega_freqs));
            u_hat_L = zeros(N, length(omega_freqs));
            u_hat_L(n,:)=temp_ud;
            
        end
    end
    
    
    
    
    %% Part 4: Saving the Modes and Center Frequencies
    
    omega_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_L
    if init_omega >0
        ii=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 2Hz
                ii=1;
            end
            n2=n2+1;
        end
        
    else
        omega_L=0;
    end
    udiff = tol+eps; % update step

    lambda = 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 modes
    u_hat_i(l,:)=u_hat_Temp(1,:,l);
    
    
    
    
    
    %%  Part 5: Stopping Criteria:
    
    if nargin >=5 % checking input of the function
        
        switch stopc
            
            case 1
                %-----------------In the Presence of Noise
                if size(u_hat_i,1) == 1
                    sigerror(l)=  norm((f_hat_onesided-(u_hat_i)),2)^2;
                else
                    sigerror(l)=  norm((f_hat_onesided-sum(u_hat_i)),2)^2;
                end
                
                if ( n2 >= 300 || sigerror(l) <= round(noisepe))
                    SC2=1;
                end
            case 2
                %-----------------Exact Reconstruction
                sum_u=sum(u_hat_Temp(1,:,:),3); % -- sum of current obtained modes
                
                normind(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;
                end
            case 3
                %------------------Bayesian Method
                if size(u_hat_i,1) == 1
                    sigerror(l)= norm((f_hat_onesided-(u_hat_i)),2)^2;
                    
                else
                    sigerror(l)= norm((f_hat_onesided-sum(u_hat_i)),2)^2;
                    
                end
                BIC(l)=2*T*log(sigerror(l))+(3*l)*log(2*T);
                
                if(l>1)
                    if(BIC(l)>BIC(l-1))
                        SC2=1;
                    end
                end
                
            otherwise
                %------------------Power of the Last Mode
                if (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));
                else
                    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(l)=polm(l)./polm_temp;
                end
                
                if (l>1 && (abs(polm(l)-polm(l-1))<0.001) )
                    SC2=1;
                end
        end
    else
        %------------------Power of the Last Mode
        if (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));
        else
            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(l)=polm(l)./polm_temp;
        end
        
        if (l>1 && (abs(polm(l)-polm(l-1))<tol) )
            SC2=1;
        end
    end
    
    
    
    
    
 %% Part 6: Resetting the counters and initializations  
    u_hat_L = zeros(N, length(omega_freqs));
    
    
    n = 1; % ----- reset the loop counter
    
    
    l=l+1; %---(number of obtained modes)+1
    m=0;
    n2=0;
end





%% ------------------ Part 7: Signal Reconstruction

omega = omega_d_Temp;
L=length(omega); %------number of modes


u_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:L
    u(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:L
    u_hat(:,l)=fftshift(fft(u(l,:)))';
end

end

（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)+1
    subplot(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',[]);
    end
    if (i == size(emd_imf,1)+1)
        ylabel('\fontname{Times new roman}RSE');
        xlabel('\fontname{Times new roman}Time/\it{s}');
    end
    grid 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)+1
    subplot(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',[]);
    end
    if (i == size(emd_imf,1)+1)
        ylabel('\fontname{Times new roman}RSE');
        xlabel('\fontname{Times new roman}Time/\it{s}');
    end
    grid 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',[]);
    end
    if (i == size(emd_imf,1)+1)
        xlabel('\fontname{Times new roman}Frequency/\it{Hz}');
    end
    grid 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

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