VMD的目标是将实值输入信号 f f f分解为离散数量的子信号(模态) u k u_k uk 。我们先假设每个模态在一个中心频率 ω k \omega_k ωk周围是紧密的，也就是说这个模态包含的频率分量都在 ω k \omega_k ωk附近，而 ω k \omega_k ωk是随着分解来确定。
为了评估一个模态的带宽，提出以下方案:1)对于每个模态，通过希尔伯特变换计算相关的分析信号，以便获得单向频谱。2)对于每种模态，通过与调谐到相应估计中心频率的指数混频，将模态的频谱移至“基带”。3)现在通过解调信号的高斯平滑度，即梯度的平方范数来估计带宽。
得到的约束变分问题如下:
求解上述方程，得到模态
u
k
u_k
uk的求解公式为：
中心频率
ω
k
\omega_k
ωk的求解公式为：
经过VMD分解之后的k个模态直接相加可以得到一个原信号的近似信号，两信号相减会有一个残差，这是因为对于一个实际信号不管分解多少次都不可能完全用分解出来的模态完全代表原信号。
所以在对分解出来的模态操作时不能忘记残差。
# -*- coding: utf-8 -*-
"""
Created on Wed Feb 20 19:24:58 2019
@author: Vinícius Rezende Carvalho
"""
import numpy as np
def VMD(f, alpha, tau, K, DC, init, tol):
"""
u,u_hat,omega = VMD(f, alpha, tau, K, DC, init, tol)
Variational mode decomposition
Python implementation by Vinícius Rezende Carvalho - vrcarva@gmail.com
code based on Dominique Zosso's MATLAB code, available at:
https://www.mathworks.com/matlabcentral/fileexchange/44765-variational-mode-decomposition
Original 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 - 即将被分解的一维时域信号
alpha - 数据保真度约束的平衡参数
tau - 双重上升的时间步长（为零噪声选择0）
K - 要被恢复的模态数量
DC - 如果将第一种模态置于并保持为直流（0频率），则为true
init - 0 = all omegas start at 0
1 = all omegas start uniformly distributed
2 = all omegas initialized randomly
tol - tolerance of convergence criterion; typically around 1e-6
收敛准则的容忍度； 通常在1e-6左右
Output:
-------
u - the collection of decomposed modes
分解模态的集合
u_hat - spectra of the modes
模态的频谱
omega - estimated mode center-frequencies
被估计的模态中心频率
omega 是一个矩阵，他储存了Niter-1组中心频率值，
形状为(Niter-1, K),Niter在vmdpy.py中定义，K为分解数量
omega矩阵储存了中心频率收敛过程中的数据，
所以一般取最后一行来用，这是最终的中心频率
"""
#
if len(f)%2:
f = f[:-1]
# 输入信号的周期和采样频率
fs = 1./len(f)
ltemp = len(f)//2
# 将f的前半截沿着axis = 0反转，再将f append到反转后的数组中
# 例如f = np.array([1,2,3,4,5,6,7,8]),处理后变成了
# array([4, 3, 2, 1, 1, 2, 3, 4, 5, 6, 7, 8])
fMirr = np.append(np.flip(f[:ltemp],axis = 0),f)
# 接着上面的例子，处理之后变成了array([4, 3, 2, 1, 1, 2, 3, 4, 5, 6, 7, 8,8,7,6,5])
fMirr = np.append(fMirr,np.flip(f[-ltemp:],axis = 0))
"""
# Time Domain 0 to T (of mirrored signal)
# 时域0到T（镜像信号的）
# 再接上面的例子，T=16
t =
array([0.0625, 0.125 , 0.1875, 0.25 , 0.3125, 0.375 , 0.4375, 0.5 ,
0.5625, 0.625 , 0.6875, 0.75 , 0.8125, 0.875 , 0.9375, 1. ])
"""
T = len(fMirr)
t = np.arange(1,T+1)/T
# Spectral Domain discretization 谱域离散化
"""
freqs =
array([-0.5 , -0.4375, -0.375 , -0.3125, -0.25 , -0.1875, -0.125 ,
-0.0625, 0. , 0.0625, 0.125 , 0.1875, 0.25 , 0.3125,
0.375 , 0.4375])
"""
freqs = 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 mode
# 为了将来的概括：每种模态都有单独的Alpha
Alpha = alpha*np.ones(K)
"""
# Construct and center f_hat 构建并居中f_hat
# np.fft.fft(fMirr)表示对fMirr进行快速傅立叶变换
# numpy.fft模块中的fftshift函数可以将FFT输出中的直流分量移动到频谱的中央。
# ifftshift函数则是其逆操作。
# 所以这一句的意思是先对fMirr进行快速傅立叶变换，然后将变换后的直流分量移到频谱中央
"""
f_hat = np.fft.fftshift((np.fft.fft(fMirr)))
f_hat_plus = np.copy(f_hat) #copy f_hat
# 把数组前半截元素都置为0
f_hat_plus[:T//2] = 0
# Initialization of omega_k
omega_plus = np.zeros([Niter, K])
if init == 1: # 1 = all omegas start uniformly distributed
for i in range(K):
omega_plus[0,i] = (0.5/K)*(i)
elif init == 2: # 2 = all omegas initialized randomly
omega_plus[0,:] = np.sort(np.exp(np.log(fs) + (np.log(0.5)-np.log(fs))*np.random.rand(1,K)))
else: # all omegas start at 0
omega_plus[0,:] = 0
# if DC mode imposed, set its omega to 0
if DC:
omega_plus[0,0] = 0
# start with empty dual variables 从空对偶变量开始
lambda_hat = np.zeros([Niter, len(freqs)], dtype = complex)
# other inits
uDiff = tol+np.spacing(1) # update step 更新步长？
n = 0 # loop counter 循环计数器
sum_uk = 0 # accumulator 累加器
# matrix keeping track of every iterant // could be discarded for mem
# 跟踪每个巡回剂的矩阵//可被mem丢弃
u_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 accumulator 更新第一模态累加器
k = 0
sum_uk = u_hat_plus[n,:,K-1] + sum_uk - u_hat_plus[n,:,0]
# update spectrum of first mode through Wiener filter of residuals
# 通过残差的维纳滤波器更新第一模态的频谱
u_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 0 如果没有保持0，则更新第一个omega
if 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 mode
for k in np.arange(1,K):
#accumulator
sum_uk = u_hat_plus[n+1,:,k-1] + sum_uk - u_hat_plus[n,:,k]
# mode spectrum
u_hat_plus[n+1,:,k] = (f_hat_plus - sum_uk - lambda_hat[n,:]/2)/(1+Alpha[k]*(freqs - omega_plus[n,k])**2)
# center frequencies
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)
# Dual ascent 双重上升
lambda_hat[n+1,:] = lambda_hat[n,:] + tau*(np.sum(u_hat_plus[n+1,:,:],axis = 1) - f_hat_plus)
# loop counter
n = 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 early
Niter = np.min([Niter,n])
omega = omega_plus[:Niter,:]
idxs = np.flip(np.arange(1,T//2+1),axis = 0)
# Signal reconstruction
u_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 part
u = u[:,T//4:3*T//4]
# recompute spectrum
u_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
（1）np.flip：沿着指定的轴反转元素
>>> A = np.arange(8).reshape((2,2,2))
>>> A
array([[[0, 1],
[2, 3]],
[[4, 5],
[6, 7]]])
>>> flip(A, 0)
array([[[4, 5],
[6, 7]],
[[0, 1],
[2, 3]]])
>>> flip(A, 1)
array([[[2, 3],
[0, 1]],
[[6, 7],
[4, 5]]])
>>>np.flip(A, 2)
array([[[1, 0],
[3, 2]],
[[5, 4],
[7, 6]]])
参考论文链接：https://pan.baidu.com/s/1BCKjc5paLkAOtM1lVlXSjQ
提取码：105f