1.VMD分解原理
很多博客都说的比较好了
2.VMD源码
vrcarva/vmdpy: Variational mode decomposition (VMD) in Python (github.com)
https://github.com/vrcarva/vmdpyGitHub已经将matlab里的实现方法用Python写出来了
并且也已经有了大佬做解读
(10条消息) 变分模态分解（VMD）运算步骤及源码解读_comli_cn的博客-CSDN博客_变分模态分解https://blog.csdn.net/comli_cn/article/details/110647039
3.源码的改进
来自
zhangzhaoweii/Variational-Modal-Decomposition-for-1D-and-nD-signal: Variational mode decomposition (VMD) for 1D and nD signal in Python (github.com)https://github.com/zhangzhaoweii/Variational-Modal-Decomposition-for-1D-and-nD-signal
从源码里得知在初始化和计算迭代
的时候，预分配的矩阵都将迭代次数Niter加入到了计算矩阵的维度里，导致空间复杂度提高
omega_plus = np.zeros([Niter, K])
lambda_hat = np.zeros([Niter, len(freqs)], dtype = complex)
u_hat_plus = np.zeros([Niter, len(freqs), K],dtype=complex) 因此考虑使用暂存矩阵的方法不断迭代交替，这样在预分配矩阵的时候就不需要考虑迭代次数这一维度
u_hat_plus_cur=np.zeros([f.shape[0], T, K],dtype=complex)
lambda_hat_cur=np.zeros([f.shape[0], T], dtype = complex)
omega_plus_cur=np.zeros([f.shape[0], K])
u_hat_plus_pre=copy.deepcopy(u_hat_plus)
lambda_hat_pre=copy.deepcopy(lambda_hat)
omega_plus_pre=copy.deepcopy(omega_plus)
u_hat_plus=copy.deepcopy(u_hat_plus_cur)
lambda_hat=copy.deepcopy(lambda_hat_cur)
omega_plus=copy.deepcopy(omega_plus_cur)import numpy as np
import copy
import scipy.fftpack as fftpack
#f中的行向量为一维信号
def VMD_for_1D_signal(f,alpha,tau,K,tol,Niter):
ltemp = len(f)//2
fs=1./len(f)
fMirr = np.append(np.flip(f[:ltemp],axis = 0),f)
fMirr = np.append(fMirr,np.flip(f[-ltemp:],axis = 0))
T = len(fMirr)
t = np.arange(1,T+1)/T
freqs = t-0.5-(1/T)
f_hat=np.zeros([1, K])
f_hat = fftpack.fftshift((fftpack.fft(fMirr)))
#对fMirr进行快速傅立叶变换，然后将变换后的直流分量移到频谱中央
f_hat_plus = np.copy(f_hat)
f_hat_plus[:T//2] = 0
omega_plus = np.zeros([1, K])
lambda_hat = np.zeros([1, T], dtype = complex)
uDiff = tol+np.spacing(1)
n = 0
sum_uk = 0
u_hat_plus = np.zeros([1, T, K],dtype=complex)
u_hat_plus_cur=np.zeros([1, T, K],dtype=complex)
lambda_hat_cur=np.zeros([1, T], dtype = complex)
omega_plus_cur=np.zeros([1, K])
#预分配矩阵空间，抛弃迭代维数
while ( uDiff > tol and n < Niter-1 ): # not converged and below iterations limit
k = 0
sum_uk = u_hat_plus[:,:,K-1] + sum_uk - u_hat_plus[:,:,0]
a=(f_hat_plus - sum_uk - lambda_hat[:]/2)
b=(1.+alpha*(freqs - omega_plus[:,k])**2)
u_hat_plus_cur[:,:,k]= a/b
# 通过残差的维纳滤波器更新第一模态的频谱
c1=abs(u_hat_plus_cur[:,T//2:T,k])**2
d1=freqs[T//2:T]
e1=np.dot(c1,d1.T)
f1=np.sum(abs(u_hat_plus_cur[:,T//2:T,k])**2)
omega_plus_cur[:,k] = e1/f1
#更新第一个omega
for k in np.arange(1,K):
sum_uk = u_hat_plus_cur[:,:,k-1] + sum_uk - u_hat_plus[:,:,k]
u_hat_plus_cur[:,:,k] = (f_hat_plus - sum_uk - lambda_hat[:]/2)/(1+alpha*
(freqs - omega_plus[:,k])**2)
c2=abs(u_hat_plus_cur[:,T//2:T,k])**2
d2=freqs[T//2:T]
e2=np.dot(c2,d2.T)
f2=np.sum(abs(u_hat_plus_cur[:,T//2:T,k])**2)
omega_plus_cur[:,k] = e2/f2
#更新k个模态的频谱和频率
lambda_hat_cur[:,:] = lambda_hat[:,:] + tau*(np.sum(u_hat_plus_cur[:,:,:],axis =
2) - f_hat_plus)
uDiff = np.spacing(1)
for i in range(K):
uDiff = uDiff + (1/T)*np.dot((u_hat_plus_cur[:,:,i]-
u_hat_plus[:,:,i]),np.transpose(np.conj((u_hat_plus_cur[:,:,i]-
u_hat_plus[:,:,i]))))
uDiff = np.abs(uDiff)
#n=n+1
u_hat_plus_pre=copy.deepcopy(u_hat_plus)
lambda_hat_pre=copy.deepcopy(lambda_hat)
omega_plus_pre=copy.deepcopy(omega_plus)
u_hat_plus=copy.deepcopy(u_hat_plus_cur)
lambda_hat=copy.deepcopy(lambda_hat_cur)
omega_plus=copy.deepcopy(omega_plus_cur)
n=n+1
#矩阵暂存迭代，节省空间，注意区分需要的矩阵
omega = omega_plus_pre
idxs = np.flip(np.arange(1,T//2+1),axis = 0)
u_hat = np.zeros([1, T, K],dtype=complex)
u_hat[:,T//2:T,:] = u_hat_plus_pre[:,T//2:T,:]
u_hat[:,idxs,:] = np.conj(u_hat_plus_pre[:,T//2:T,:])
u_hat[:,0,:] = np.conj(u_hat_plus_pre[:,-1,:])
u = np.zeros([1,K,T])
for k in range(K):
u[:,k,:] = np.real(np.fft.ifft(np.fft.ifftshift(u_hat[:,:,k])))
u = u[:,:,T//4:3*T//4]
u_hat = np.zeros([f.shape[0],u.shape[2], K],dtype=complex)
for k in range(K):
u_hat[:,:,k] = fftpack.fftshift(fftpack.fft(u[:,k,:]), axes=(1,))
#后处理与信号还原
return u, u_hat, omega同样地，利用这个思想，设计计算多维信号的方法，每个一维信号作为矩阵的一组
import numpy as np
import copy
import scipy.fftpack as fftpack
#f中的每一行为一维信号
def VMD_for_nD_signal(f,alpha,tau,K,tol,Niter):
ltemp = f.shape[1]//2
fs=1./f.shape[1]
fMirr1=np.flip(f[:,ltemp:],1)
fMirr2=np.flip(f[:,:ltemp],1)
fMirr =np.concatenate((fMirr2,f),axis=1)
fMirr =np.concatenate((fMirr,fMirr1),axis=1)
T=fMirr.shape[1]
t1 = np.arange(1,T+1)/T
t=np.tile(t1,(f.shape[0],1))
freqs = t-0.5-(1/T)
f_hat=np.zeros([f.shape[0], K])
f_hat= fftpack.fftshift((fftpack.fft(fMirr)), axes=(1,))
f_hat_plus = np.copy(f_hat) #copy f_hat
f_hat_plus[:,:T//2] = 0
omega_plus = np.zeros([f.shape[0], K])
lambda_hat = np.zeros([f.shape[0], T], dtype = complex)
lambda_hat_temp=copy.deepcopy(lambda_hat)
uDiff = tol+np.spacing(1)
uDiff_max=uDiff
n = 0
sum_uk = 0
u_hat_plus = np.zeros([f.shape[0], T, K],dtype=complex)
u_hat_plus_cur=np.zeros([f.shape[0], T, K],dtype=complex)
lambda_hat_cur=np.zeros([f.shape[0], T], dtype = complex)
omega_plus_cur=np.zeros([f.shape[0], K])
#预分配矩阵空间，抛弃迭代维数,需要 f.shape[0]组
while ( uDiff_max > tol and n < Niter-1 ):
k = 0
sum_uk = u_hat_plus[:,:,K-1] + sum_uk - u_hat_plus[:,:,0]
a=(f_hat_plus - sum_uk - lambda_hat[:]/2)
b=(1.+alpha*(freqs - np.tile(omega_plus[:,k].reshape(f.shape[0],1),(1,T)))**2)
u_hat_plus_cur[:,:,k]= a/b
c1=abs(u_hat_plus_cur[:,T//2:T,k])**2
d1=freqs[:,T//2:T]
e1=np.diagonal(np.dot(c1,d1.T))
f1=np.sum(abs(u_hat_plus_cur[:,T//2:T,k])**2,1)
omega_plus_cur[:,k] = e1/f1
for k in np.arange(1,K):
sum_uk = u_hat_plus_cur[:,:,k-1] + sum_uk - u_hat_plus[:,:,k]
a=(f_hat_plus - sum_uk - lambda_hat[:]/2)
b=(1.+alpha*(freqs - np.tile(omega_plus[:,k].reshape(f.shape[0],1),(1,T)))**2)
u_hat_plus_cur[:,:,k]= a/b
c2=abs(u_hat_plus_cur[:,T//2:T,k])**2
d2=freqs[:,T//2:T]
e2=np.diagonal(np.dot(c2,d2.T))
f2=np.sum(abs(u_hat_plus_cur[:,T//2:T,k])**2,1)
omega_plus_cur[:,k] = e2/f2
lambda_hat_cur[:,:] = lambda_hat[:,:] + tau*(np.sum(u_hat_plus_cur[:,:,:],axis = 2) - f_hat_plus)
uDiff = np.spacing(1)
for i in range(K):
uDiff = uDiff + (1/T)*np.dot((u_hat_plus_cur[:,:,i]-u_hat_plus[:,:,i]),np.transpose(np.conj((u_hat_plus_cur[:,:,i]-u_hat_plus[:,:,i]))))
uDiff = np.abs(uDiff)
uDiff_max=max(np.diagonal(uDiff))
u_hat_plus_pre=copy.deepcopy(u_hat_plus)
lambda_hat_pre=copy.deepcopy(lambda_hat)
omega_plus_pre=copy.deepcopy(omega_plus)
u_hat_plus=copy.deepcopy(u_hat_plus_cur)
lambda_hat=copy.deepcopy(lambda_hat_cur)
omega_plus=copy.deepcopy(omega_plus_cur)
n=n+1
omega = omega_plus_pre
idxs = np.flip(np.arange(1,T//2+1),axis = 0)
u_hat = np.zeros([f.shape[0], T, K],dtype=complex)
u_hat[:,T//2:T,:] = u_hat_plus_pre[:,T//2:T,:]
u_hat[:,idxs,:] = np.conj(u_hat_plus_pre[:,T//2:T,:])
u_hat[:,0,:] = np.conj(u_hat_plus_pre[:,-1,:])
u = np.zeros([f.shape[0],K,T])
for k in range(K):
u[:,k,:] = np.real(fftpack.ifft(fftpack.ifftshift(u_hat[:,:,k], axes=(1,))))
u = u[:,:,T//4:3*T//4]
u_hat = np.zeros([f.shape[0],u.shape[2], K],dtype=complex)
for k in range(K):
u_hat[:,:,k] = fftpack.fftshift(fftpack.fft(u[:,k,:]), axes=(1,))
return u, u_hat, omega
对于需要一次性计算的多组信号数据，将每一组作为矩阵f的行向量
