购物网站流量怎么做,医院网站域名备案,排版设计工作内容,wordpress自定义过滤大家好#xff0c;我是带我去滑雪#xff01; 序贯变分模态分解(SVMD) 是一种信号处理和数据分析方法。它可以将复杂信号分解为一系列模态函数#xff0c;每个模态函数代表信号中的特定频率分量。 SVMD 的主要目标是提取信号中的不同频率分量并将其重构为原始信号。SVMD的基… 大家好我是带我去滑雪 序贯变分模态分解(SVMD) 是一种信号处理和数据分析方法。它可以将复杂信号分解为一系列模态函数每个模态函数代表信号中的特定频率分量。 SVMD 的主要目标是提取信号中的不同频率分量并将其重构为原始信号。SVMD的基本原理是通过变分模态分解的方式将信号分解为多个模态函数。在每个迭代步骤中SVMD 通过最小化信号和模态函数之间的差异来更新模态函数。重复这个过程直到收敛。得到的模态函数可用于重建原始信号。 SVMD 的另一个关键特征是连续分解。在每个迭代步骤中SVMD 从信号中提取主频率分量并将其从信号中删除。这样每次迭代步骤都会提取信号中的一个频率分量直到提取完所有频率分量。这种逐次分解方法可以更好地捕获信号中的不同频率分量。SVMD 在信号处理和数据分析方面有着广泛的应用。它可用于去噪、特征提取、频谱分析等多个领域。通过将信号分解为模态函数SVMD可以更好地理解和描述信号的频率特性。这对于信号处理和数据分析非常重要。SVMD的数据重构是将分解后的模态函数重新组合成原始信号的过程。通过各模态函数的加权相加即可得到重构信号。这个过程可以用来恢复原始信号的频率特性并且可以根据需要进一步分析和处理。 综上所述逐次变分模态分解是一种有效的信号处理和数据分析方法。它可以将复杂信号分解为多个模态函数并可以通过数据重构将它们重新组合成原始信号。 SVMD有着广泛的应用范围对于理解和描述信号的频率特性非常有帮助。通过深入研究和应用SVMD我们可以更好地处理和分析各类信号和数据。下面开始代码实战。
1SVMD实现
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
signoisesignal-y; %-estimating the noisesave_T length(signal);
fs 1/save_T;%______________________________________________________________________
%
% Mirroring the signal and noise part to extend
%______________________________________________________________________
T save_T;
f_mirzeros(1,T/2);
f_mir_noisezeros(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/21:3*T/2) signal;
f_mir_noise(T/21:3*T/2) signoise;
f_mir(3*T/21:2*T) signal(T:-1:T/21);
f_mir_noise(3*T/21:2*T) signoise(T:-1:T/21);f f_mir;
fnoisef_mir_noise;
%______________________________________________________________________
%______________________________________________________________________T length(f);%------------- time domain (t -- 0 to T)
t (1:T)/T;udiff toleps; %------ update stepomega_freqs t-0.5-1/T;%------------- discretization of spectral domain%______________________________________________________________________
%
% FFT of signal(and Hilbert transform conceptmaking 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;
%______________________________________________________________________
%______________________________________________________________________noisepenorm(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;
endminAlpha10; %------ the initial value of alpha
AlphaminAlpha; %------ the initial value of alpha
alphazeros(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 counterm0; %------ iteration counter for increasing alpha
SC20; % ------ main stopping criteria index
l1; %------ the initial number of modes
bf0; % ----- bit flag to increase alpha
BICzeros(1,1); % ------- the initial value of Bayesian indexh_hat_Tempzeros(2, length(omega_freqs));%-initialization of filter matrixu_hat_Tempzeros(1,length(omega_freqs),1);%- matrix1 of modes
u_hat_izeros(1, length(omega_freqs));%- matrix2 of modesn20; % ---- counter for initializing omega_Lpolmzeros(2,1); % ---- initializing Power of Last Mode indexomega_d_Tempzeros(1,1);%-initialization of center frequencies vector1
sigerrorzeros(1,1);%initializing signal error index for stopping criteria
gammazeros(1,1);%----initializing gamma
normindzeros(1,1);%% ---------------------- Part 2: Main loop for iterative updates
while (SC2~1)while (Alpha(1,1)(maxAlpha1)) while ( udiff tol n N ) %------------------ update uLu_hat_L(n1,:) (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(n1,1) (omega_freqs(T/21:T)*(abs(u_hat_L(n1, T/21:T)).^2))/sum(abs(u_hat_L(n1,T/21:T,1)).^2);%------------------ update lambda (dual ascent)lambda(n1,:) lambda(n,:) tau*(f_hat_onesided...-(u_hat_L(n1,:) (((Alpha(1,1).^2)*(omega_freqs - omega_L(n,1)).^4.....*(f_hat_onesided - u_hat_L(n1,:)-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(n1,:)-u_hat_L(n,:))*conj((u_hat_L(n1,:)-u_hat_L(n,:)))).../ (1/T*(u_hat_L(n,:))*conj((u_hat_L(n,:))));udiff abs(udiff);n n1;end%% ---- Part 3: Increasing Alpha to achieve a pure modeif abs(m-log(maxAlpha)) 1mm1;elsemm.05;bfbf1;endif bf2AlphaAlpha1;endif Alpha(1,1)(maxAlpha-1) %exp(SC1)(maxAlpha)if (bf 1)Alpha(1,1)maxAlpha-1;elseAlpha(1,1)exp(m);endomega_Lomega_L(n,1);% ------- Initializingudiff toleps; % 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_Lomega_L(omega_L0);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;bf0;%------------------------------initializing omega_Lif init_omega 0ii0;while (ii1 n2 300)omega_L sort(exp(log(fs) (log(0.5)-log(fs))*rand(1,1)));checkpabs(omega_d_Temp-omega_L);if (size(find(checkp0.02),2)0) % it will continue if difference between previous vector of omega_d and the current random omega_plus is about 2Hzii1;endn2n21;endelseomega_L0;endudiff toleps; % 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))SC21;endcase 2%-----------------Exact Reconstructionsum_usum(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 )SC21;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(l1)if(BIC(l)BIC(l-1))SC21;endendotherwise%------------------Power of the Last Modeif (l2)polm(l)norm((4*Alpha(1,1)*u_hat_i(l,:)./(12*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:),2);polm_temppolm(l);polm(l)polm(l)./max(polm(l));elsepolm(l)norm((4*Alpha(1,1)*u_hat_i(l,:)./(12*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:),2);polm(l)polm(l)./polm_temp;endif (l1 (abs(polm(l)-polm(l-1))0.001) )SC21;endendelse%------------------Power of the Last Modeif (l2)polm(l)norm((4*Alpha(1,1)*u_hat_i(l,:)./(12*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:),2);polm_temppolm(l);polm(l)polm(l)./max(polm(l));elsepolm(l)norm((4*Alpha(1,1)*u_hat_i(l,:)./(12*Alpha(1,1)*...(omega_freqs-omega_d_Temp(l)).^2))*u_hat_i(l,:),2);polm(l)polm(l)./polm_temp;endif (l1 (abs(polm(l)-polm(l-1))tol) )SC21;endend%% Part 6: Resetting the counters and initializations u_hat_L zeros(N, length(omega_freqs));n 1; % ----- reset the loop counterll1; %---(number of obtained modes)1m0;n20;
end%% ------------------ Part 7: Signal Reconstructionomega omega_d_Temp;
Llength(omega); %------number of modesu_hat zeros(T, L);
u_hat((T/21):T,:) squeeze(u_hat_Temp(1,(T/21):T,:));
u_hat((T/21):-1:2,:) squeeze(conj(u_hat_Temp(1,(T/21):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);
uu(indic,:);
%---------- remove mirror part
u u(:,T/41:3*T/4);%--------------- recompute spectrum
clear u_hat;
for l 1:Lu_hat(:,l)fftshift(fft(u(l,:)));
endend2SVMD绘图
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
SIGimportdata(NASA电容量.csv);
sigSIG(2:161,2);%%想要分解哪一列就填几
%(1)导入时间数据来设置时间
tSIG(2:161,2);
Fs1/(t(2)-t(1));
%(2)设置采样率来设置时间
% Nlength(sig);
% Fs1000;%%采样频率自己设置
% 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});
maxAlpha1000; %compactness of mode
tau0;%time-step of the dual ascent
tol1e-6; %tolerance of convergence criterion;
stopc4;%the type of stopping criteria
[svmd_imf,uhat,omega]svmd(sig,maxAlpha,tau,tol,stopc);
Huatu_svmd(svmd_imf,sig,t);
svmd_imfsvmd_imf;
需要数据集的家人们可以去百度网盘永久有效获取
链接https://pan.baidu.com/s/173deLlgLYUz789M3KHYw-Q?pwd0ly6 提取码2138 更多优质内容持续发布中请移步主页查看。
博主的WeChat:TCB1736732074 点赞关注,下次不迷路
