变分模态分解（VMD）运算步骤及源码解读

2021/9/5 17:06:49

编程Tag： 源码模态 plus np Omega hat 变分 Niter VMD

本文主要是介绍变分模态分解（VMD）运算步骤及源码解读，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

1. 简述

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的求解公式为：
在这里插入图片描述

2. 运算步骤

在这里插入图片描述

3.VMD分解之后的模态和原信号的关系

经过VMD分解之后的k个模态直接相加可以得到一个原信号的近似信号，两信号相减会有一个残差，这是因为对于一个实际信号不管分解多少次都不可能完全用分解出来的模态完全代表原信号。
所以在对分解出来的模态操作时不能忘记残差。

4. 源码解读

# -*- 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

5. 相关知识点

（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

这篇关于变分模态分解（VMD）运算步骤及源码解读的文章就介绍到这儿，希望我们推荐的文章对大家有所帮助，也希望大家多多支持为之网！