信号处理VMD 变分模态分解，示例+完整代码

首先， pip install vmdpy

下面两段代码，1为VMD函数的定义

2为一个模态分解的示例

1111111111111111111111111

# -*- coding: utf-8 -*-
"""
Created on Wed Feb 20 19:24:58 2019@author: Vinícius Rezende Carvalho
"""
import numpy as npdef  VMD(f, alpha, tau, K, DC, init, tol):"""u,u_hat,omega = VMD(f, alpha, tau, K, DC, init, tol)Variational mode decompositionPython implementation by Vinícius Rezende Carvalho - vrcarva@gmail.comcode based on Dominique Zosso's MATLAB code, available at:https://www.mathworks.com/matlabcentral/fileexchange/44765-variational-mode-decompositionOriginal paper:Dragomiretskiy, K. and Zosso, D. (2014) ‘Variational Mode Decomposition’, IEEE Transactions on Signal Processing, 62(3), pp. 531–544. doi: 10.1109/TSP.2013.2288675.Input and Parameters:---------------------f       - the time domain signal (1D) to be decomposedalpha   - the balancing parameter of the data-fidelity constrainttau     - time-step of the dual ascent ( pick 0 for noise-slack )K       - the number of modes to be recoveredDC      - true if the first mode is put and kept at DC (0-freq)init    - 0 = all omegas start at 01 = all omegas start uniformly distributed2 = all omegas initialized randomlytol     - tolerance of convergence criterion; typically around 1e-6Output:-------u       - the collection of decomposed modesu_hat   - spectra of the modesomega   - estimated mode center-frequencies"""if len(f)%2:f = f[:-1]# Period and sampling frequency of input signalfs = 1./len(f)ltemp = len(f)//2 fMirr =  np.append(np.flip(f[:ltemp],axis = 0),f)  fMirr = np.append(fMirr,np.flip(f[-ltemp:],axis = 0))# Time Domain 0 to T (of mirrored signal)T = len(fMirr)t = np.arange(1,T+1)/T  # Spectral Domain discretizationfreqs = t-0.5-(1/T)# Maximum number of iterations (if not converged yet, then it won't anyway)Niter = 500# For future generalizations: individual alpha for each modeAlpha = alpha*np.ones(K)# Construct and center f_hatf_hat = np.fft.fftshift((np.fft.fft(fMirr)))f_hat_plus = np.copy(f_hat) #copy f_hatf_hat_plus[:T//2] = 0# Initialization of omega_komega_plus = np.zeros([Niter, K])if init == 1:for i in range(K):omega_plus[0,i] = (0.5/K)*(i)elif init == 2:omega_plus[0,:] = np.sort(np.exp(np.log(fs) + (np.log(0.5)-np.log(fs))*np.random.rand(1,K)))else:omega_plus[0,:] = 0# if DC mode imposed, set its omega to 0if DC:omega_plus[0,0] = 0# start with empty dual variableslambda_hat = np.zeros([Niter, len(freqs)], dtype = complex)# other initsuDiff = tol+np.spacing(1) # update stepn = 0 # loop countersum_uk = 0 # accumulator# matrix keeping track of every iterant // could be discarded for memu_hat_plus = np.zeros([Niter, len(freqs), K],dtype=complex)    #*** Main loop for iterative updates***while ( uDiff > tol and  n < Niter-1 ): # not converged and below iterations limit# update first mode accumulatork = 0sum_uk = u_hat_plus[n,:,K-1] + sum_uk - u_hat_plus[n,:,0]# update spectrum of first mode through Wiener filter of residualsu_hat_plus[n+1,:,k] = (f_hat_plus - sum_uk - lambda_hat[n,:]/2)/(1.+Alpha[k]*(freqs - omega_plus[n,k])**2)# update first omega if not held at 0if not(DC):omega_plus[n+1,k] = np.dot(freqs[T//2:T],(abs(u_hat_plus[n+1, T//2:T, k])**2))/np.sum(abs(u_hat_plus[n+1,T//2:T,k])**2)# update of any other modefor k in np.arange(1,K):#accumulatorsum_uk = u_hat_plus[n+1,:,k-1] + sum_uk - u_hat_plus[n,:,k]# mode spectrumu_hat_plus[n+1,:,k] = (f_hat_plus - sum_uk - lambda_hat[n,:]/2)/(1+Alpha[k]*(freqs - omega_plus[n,k])**2)# center frequenciesomega_plus[n+1,k] = np.dot(freqs[T//2:T],(abs(u_hat_plus[n+1, T//2:T, k])**2))/np.sum(abs(u_hat_plus[n+1,T//2:T,k])**2)# Dual ascentlambda_hat[n+1,:] = lambda_hat[n,:] + tau*(np.sum(u_hat_plus[n+1,:,:],axis = 1) - f_hat_plus)# loop countern = n+1# converged yet?uDiff = np.spacing(1)for i in range(K):uDiff = uDiff + (1/T)*np.dot((u_hat_plus[n,:,i]-u_hat_plus[n-1,:,i]),np.conj((u_hat_plus[n,:,i]-u_hat_plus[n-1,:,i])))uDiff = np.abs(uDiff)        #Postprocessing and cleanup#discard empty space if converged earlyNiter = np.min([Niter,n])omega = omega_plus[:Niter,:]idxs = np.flip(np.arange(1,T//2+1),axis = 0)# Signal reconstructionu_hat = np.zeros([T, K],dtype = complex)u_hat[T//2:T,:] = u_hat_plus[Niter-1,T//2:T,:]u_hat[idxs,:] = np.conj(u_hat_plus[Niter-1,T//2:T,:])u_hat[0,:] = np.conj(u_hat[-1,:])    u = np.zeros([K,len(t)])for k in range(K):u[k,:] = np.real(np.fft.ifft(np.fft.ifftshift(u_hat[:,k])))# remove mirror partu = u[:,T//4:3*T//4]# recompute spectrumu_hat = np.zeros([u.shape[1],K],dtype = complex)for k in range(K):u_hat[:,k]=np.fft.fftshift(np.fft.fft(u[k,:]))return u, u_hat, omega

2222222222222222222222222222222222

import numpy as np
import matplotlib.pyplot as plt
from vmdpy import VMD
# Time Domain 0 to T
T = 1000
fs = 1/T
t = np.arange(1,T+1)/T
freqs = 2*np.pi*(t-0.5-fs)/(fs)
# center frequencies of components
f_1 = 2
f_2 = 24
f_3 = 288
# modes
v_1 = (np.cos(2*np.pi*f_1*t))
v_2 = 1/4*(np.cos(2*np.pi*f_2*t))
v_3 = 1/16*(np.cos(2*np.pi*f_3*t))
#
#% for visualization purposes
fsub = {1:v_1,2:v_2,3:v_3}
wsub = {1:2*np.pi*f_1,2:2*np.pi*f_2,3:2*np.pi*f_3}# composite signal, including noisef = v_1 + v_2 + v_3 + 0.1*np.random.randn(v_1.size)
f_hat = np.fft.fftshift((np.fft.fft(f)))# some sample parameters for VMD
alpha = 2000       # moderate bandwidth constraint
tau = 0.            # noise-tolerance (no strict fidelity enforcement)
K = 3              # 3 modes
DC = 0             # no DC part imposed
init = 1           # initialize omegas uniformly
tol = 1e-7
# Run actual VMD code
u, u_hat, omega = VMD(f, alpha, tau, K, DC, init, tol)
#%%
# Simple Visualization of decomposed modes
plt.figure()
plt.plot(u.T)
plt.title('Decomposed modes')
# For convenience here: Order omegas increasingly and reindex u/u_hat
sortIndex = np.argsort(omega[-1,:])
omega = omega[:,sortIndex]
u_hat = u_hat[:,sortIndex]
u = u[sortIndex,:]
linestyles = ['b', 'g', 'm', 'c', 'c', 'r', 'k']fig1 = plt.figure()
plt.subplot(411)
plt.plot(t,f)
plt.xlim((0,1))
for key, value in fsub.items():plt.subplot(4,1,key+1)plt.plot(t,value)
fig1.suptitle('Original input signal and its components')fig2 = plt.figure()
plt.loglog(freqs[T//2:], abs(f_hat[T//2:]))
plt.xlim(np.array([1,T/2])*np.pi*2)
ax = plt.gca()
ax.grid(which='major', axis='both', linestyle='--')
fig2.suptitle('Input signal spectrum')fig3 = plt.figure()
for k in range(K):plt.semilogx(2*np.pi/fs*omega[:,k], np.arange(1,omega.shape[0]+1), linestyles[k])
fig3.suptitle('Evolution of center frequencies omega')fig4 = plt.figure()
plt.loglog(freqs[T//2:], abs(f_hat[T//2:]), 'k:')
plt.xlim(np.array([1, T//2])*np.pi*2)
for k in range(K):plt.loglog(freqs[T//2:], abs(u_hat[T//2:,k]), linestyles[k])
fig4.suptitle('Spectral decomposition')
plt.legend(['Original','1st component','2nd component','3rd component'])fig4 = plt.figure()for k in range(K):plt.subplot(3,1,k+1)plt.plot(t,u[k,:], linestyles[k])plt.plot(t, fsub[k+1], 'k:')plt.xlim((0,1))plt.title('Reconstructed mode %d'%(k+1))
plt.show()
print('finish')

信号处理VMD 变分模态分解，示例+完整代码相关推荐

matlab中使用VMD(变分模态分解)
最近我们被客户要求撰写关于VMD(变分模态分解)的研究报告,包括一些图形和统计输出. 拨号音信号的变模分解创建一个以4 kHz采样的信号,类似于拨打数字电话的所有键.将信号另存为MATLAB®时间数 ...
变分模态分解 python_浅谈VMD(变分模态分解)
学号:19011210554 姓名:袁博 [嵌牛导读]:好多人看着VMD看博客最想知道的就是这东西的应用和大概步骤原理,而具体原理算法不太感兴趣,而且也不太容易看懂.本文既然是浅谈,就讲解一下VM ...
vmd变分模态分解程序matlab论坛_博士兼职辅导员论坛分享会第三期
新一期经验分享报告会又来了!本次报告满满干货,快来看看都有些什么内容吧? 报告题目齿轮箱关键零部件复合故障特征提取方法研究报告摘要: 基于振动信号的复合故障特征提取技术一直以来都是旋转机械故障诊断 ...
linux环境vmd下载,VMD Linux版下载|VMD(变分模态分解程序) V1.9.3 Linux版下载_当下软件园_软件下载...
VMD Linux版是款适用于Linux操作系统的分子运动绘图分析模拟软件.它可以真实的模拟分析的运动场景,提供可视化的图形界面,帮助用户更好的理解分析变化,提高研究效率,操作简单,方便快捷,非常好用 ...
鲸鱼算法优化变分模态分解(VMD)包络熵和参数的特征提取及MATLAB代码实现
目录 1 简介 2 变分模态分解VMD原理 3 鲸鱼优化算法优化VMD原理 3.1. 鲸鱼优化算法优化VMD原理及流程 3.2. 特征提取流程 4 优化效果 4.1. VMD各分量信号时域图 4.2. ...
变分模态分解_Android小部件示例中的模态对话框（弹出）
变分模态分解在此示例中,我们将看到如何在主屏幕中创建一个可以打开弹出对话框的Android小部件. 如您所知,Android Widgets是小型应用程序,基本上可以做两件事. 按下时启动新的活动, ...
利用智能算法优化参数的自适应变分模态分解，VMD实现混合储能系统的分频
关键词:混合储能,VMD,麻雀搜索算法,遗传算法,混合储能容量配置优化,混合储能功率分配,利用智能算法优化参数的自适应变分模态分解,VMD实现混合储能系统的分频,高频分配给超级电容器,低频分配给蓄电池 ...
基于VMD变分模态分解算法Python程序
基于VMD变分模态分解算法Python程序可用于时间序列和其他领域特色:1.基于Python 2.数据从excel文件中读取,更换简单全部完整的代码,保证可以运行的代码看这里. http://t ...
分解得到的时频域特征_【推荐文章】基于变分模态分解和广义Warblet变换的齿轮故障诊断...
<机械传动>2018年第42卷第7期文章编号:1004-2539(2018)07-0157-05 DOI:10.16578/j.issn.1004.2539.2018.07.0 ...
【VMD-SSA-LSSVM】基于变分模态分解与麻雀优化Lssvm的负荷预测【多变量】（Matlab代码实现）

信号处理VMD 变分模态分解，示例+完整代码

信号处理VMD 变分模态分解，示例+完整代码相关推荐

最新文章

热门文章