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

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/770913.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

异地两台电脑如何共享文件?

在当前数字化时代&#xff0c;人们对于数据的使用和管理变得越来越便捷。由于工作和生活的需要&#xff0c;我们常常需要在异地的电脑间共享文件。这给我们的工作和生活带来了一定程度的不便。有没有一种便捷的方法可以让异地的电脑实现文件的共享呢&#xff1f;答案是肯定的。…

flutter 父组件调用子组件方法

标题在Flutter中&#xff0c;父组件可以通过GlobalKey来引用子组件&#xff0c;并调用子组件的方法。以下是一个简单的例子&#xff1a; 在这个例子中&#xff0c;ParentComponent 有一个GlobalKey&#xff0c;它被传递给了ChildComponent。当按钮被点击时&#xff0c;通过chi…

06 mybatis </sql>

文章目录 products.sqlpom.xmlmybatis-config.xmlProductsMapper.xmlProductsMapperImpl.javaProducts.javaDButil.javaProductsMapperImplTest.javaMapperTest.java products.sql create table products (product_id int auto_increment comment 产品IDprimary key,prod…

知识图谱-图数据库-neo4j (1)踩坑记录

1、neo4j 安装 材料 &#xff1a; openjdk11 (neo4j 最低jdk版本要求) neo4j-community-4.4.30 CentOS 7.8 Release Date: 25 January 2024 Neo4j 4.4.30 is a maintenance release with many important improvements and fixes. Neo4j Deployment Center - Graph Database…

vuex状态管理的使用

一、创建store,单个store的使用 1、 /*** 该文件用于创建vuex中最核心的store*///引入Vuex import Vuex from vuex; import Vue from "vue";//使用vuex来集中管理状态,必要 //new store的前提是必须要使用Vuex插件 Vue.use(Vuex);//创建actions(本质就是对象) 用于…

【前端面试3+1】01闭包、跨域

一、对闭包的理解 定义&#xff1a; 闭包是指在一个函数内部定义的函数&#xff0c;并且该内部函数可以访问外部函数的变量。闭包使得函数内部的变量在函数执行完后仍然可以被访问和操作。 特点&#xff1a; 闭包可以访问外部函数的变量&#xff0c;即使外部函数已经执行完毕。…

vue项目中使用vue-pdf或pdf.Js,实现在页面上预览pdf内容

一。vue-pdf 1. 安装vue-pdf npm install --save vue-pdf2.页面引入 js部分 import pdf from "vue-pdf";data(){return {pdfUrl: "",pageTotal: 0,} }mounted(){this.pdfUrl pdf.createLoadingTask(pdf文件路径url);// 获取页码this.pdfUrl.promise…

Linux 搭建jenkins docker

jekin docker gitee docker 安装 jenkins docker run -d --restartalways \ --name jenkins -uroot -p 10340:8080 \ -p 10341:50000 \ -v /home/docker/jenkins:/var/jenkins_home \ -v /var/run/docker.sock:/var/run/docker.sock \ -v /usr/bin/docker:/usr/bin/docker je…

QT数据类型和容器用法

Qt库提供了基于通用模板的容器类, 这些类可用于存储指定类型的数据项&#xff0c;Qt中这些容器类的设计比STL容器更轻&#xff0c;更安全且更易于使用。容器类也都是隐式共的&#xff0c;它们是可重入的&#xff0c;并且已针对速度/低内存消耗和最小的内联代码扩展进行了优化&a…

【解析几何】 【多源路径】 【贪心】1520 最多的不重叠子字符串

作者推荐 视频算法专题 本身涉及知识点 解析几何 图论 多源路径 贪心 LeetCode1520. 最多的不重叠子字符串 给你一个只包含小写字母的字符串 s &#xff0c;你需要找到 s 中最多数目的非空子字符串&#xff0c;满足如下条件&#xff1a; 这些字符串之间互不重叠&#xff0…

简单讲讲spring事务的传播机制

事务传播机制就像是一个指挥家&#xff0c;控制着程序中的各种操作&#xff08;比如修改数据库&#xff09;何时开始、何时结束&#xff0c;以及如何处理错误。 保证数据一致性&#xff1a;想象一下你在网上购物&#xff0c;你需要先从银行账户扣款&#xff0c;然后再把商品加入…

Wireshark使用实训---分析IP包

1.Wireshark简介和作用 Wireshark是一个开源的网络分析工具&#xff0c;用于捕捉和分析网络数据包。它可以帮助网络管理员和安全专家监控和解决网络问题&#xff0c;同时也可以用于学习和教学网络通信原理。 Wireshark可以在网络中捕获和分析传输的数据包&#xff0c;包括协议…

【Java初阶(五)】类和对象

❣博主主页: 33的博客❣ ▶文章专栏分类: Java从入门到精通◀ &#x1f69a;我的代码仓库: 33的代码仓库&#x1f69a; 目录 1. 前言2.面向对象的认识3.类的认识4. 类的实例化4.1什么是实例化4.2类和对象的说明 5.this引用6.对象初始化6.1 构造方法 7.static关键字8.代码块8.1 …

探索Python中的集成方法:Stacking

在机器学习领域&#xff0c;Stacking是一种高级的集成学习方法&#xff0c;它通过将多个基本模型的预测结果作为新的特征输入到一个元模型中&#xff0c;从而提高整体模型的性能和鲁棒性。本文将深入介绍Stacking的原理、实现方式以及如何在Python中应用。 什么是Stacking&…

算法-数据结构

算法-数据结构 金无足赤人无完人&#xff0c;在处理实际问题的时候我们可以使用到很多合适的数据结构&#xff0c;但目前还没有一个数据结构可以称的上完美。查询速度快的&#xff0c;插入的速度就会慢&#xff1b;插入速度和查询速度都快得&#xff0c;占用的空间就会多&…

PTAxt的考研路

xt是我院19级专业第一&#xff0c;但他认为保研并不能展示他全部的实力&#xff0c;所以他在22年初试一结束就加入了23考研的队伍中&#xff0c;并且他为了填补我院近些年来无北大研究生的空白&#xff0c;毅然决然决定扛起19级的大旗&#xff0c;在学校百年华诞之际献上他最诚…

Springboot项目之mybatis-plus多容器分布式部署id重复问题之源码解析

mybatis-plus 3.3.2 部署多个pod id冲突问题 配置&#xff1a; # 设置随机 mybatis-plus.global-config.worker-id: ${random.int(1,31)} mybatis-plus.global-config.datacenter-id: ${random.int(1,31)}源码解析&#xff1a;MybatisSqlSessionFactoryBean 重点&#xff1a…

mysql数据库查询

MYSQL数据库的搭建 今日目标&#xff1a; 1.搭建数据库 2.实现数据库的增删改查 00-回顾 #dos的常用指令 1. 切换盘符&#xff1a; 盘符名: 2. 切换上一级&#xff1a; cd ../ 3. 切换下一级&#xff1a; cd 目录名 4. 查看当前目录下的所有子目录和子文件&#xff1a; di…

光明源@智慧公厕赋能“厕所革命”主要体现在哪些方面?

当我们提及厕所&#xff0c;不再仅是简单的卫生设施&#xff0c;而是一种对生活品质的关怀与呵护。智慧公厕&#xff0c;作为厕所革命的引领者&#xff0c;以其独特的拟人魅力&#xff0c;彰显着人性化关怀的新风尚。今日&#xff0c;让我们一同探索&#xff0c;智慧公厕是如何…

数据库备份工具(实现数据定时覆盖)

数据库备份工具&#xff08;实现数据定时覆盖&#xff09; 永远热爱&#xff0c;永远执着&#xff01; 工具介绍 自动化测试数据库更新调度程序 这段 Python 脚本自动化了每天定时从生产数据库更新测试数据库的过程。它利用了 schedule 库来安排并执行每天指定时间的更新任务…