大家好,我是带我去滑雪!
序贯变分模态分解(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
点赞+关注,下次不迷路!