标签(空格分隔)： 数学
这一篇写一下变分模态分解(原始论文：Variational Mode Decomposition)，跟原始论文思路思路一致但有一点点不太一样，原始论文写的很好，但我不是通信专业没有学过信号相关课程一开始看起来有点费劲。模态分解认为信号是由不同“模态”的子信号叠加而成的，而变分模态分解则认为信号是由不同频率占优的子信号叠加而成的，其目的是要把信号分解成不同频率的子信号。变分模态分解的分解结果如图所示
基础
一开始论文看不懂的原因是缺少相关前置知识，但一旦顺下来就会感觉其实没有那么难，难的是作者的思路很巧妙，先写下我遇到的这些知识盲点
第一点是傅立叶变换的微分性质，
 的傅立叶变换为 
 ，其导数 
 的傅立叶变换为 
另一点是解析信号，现实世界只能采集实信号，但实信号有很多不好用的性质，如存在负频率，无法直接得到调制频率后的实信号等。 设原始信号是一个实信号
为了方便表示 
 为随 
 变化的函数，相当于瞬时频率，解析信号是一个复信号，可以通过希尔伯特变换得到  
  解析信号的实部是原本的实数信号，并且经过调频之后复数信号的实部仍然是调频之后的实信号，如对信号 
 增加频率 
 ，只需要乘以 
 即可 
  由此，其实部相当于在原本频率 
 的基础上增加了频率 
 ，如下面的matlab脚本
clear;close all;clc;
t = 1:0.01:10;
%%
f1 = sin(20*t).*(t-5).^2;
subplot(3,1,1);
plot(f1);
ylim([-25 25]);
%%
f2 = sin(50*t).*(t-5).^2;
subplot(3,1,2);
plot(f2);
ylim([-25 25]);
%%
H = hilbert(f1);
f_hat = H.*exp(1i*30.*t);
subplot(3,1,3);
plot(real(f_hat));
ylim([-25 25]);
matlab的hilbert函数包括希尔伯特变换和解析函数转换两部分，直接得到实信号的解析信号，其中希尔伯特变换
正文
接下来我们看看如何一步一步得到变分模态分解的思路
原论文通过一个信号降噪问题进行说明，现需要对采样信号
 进行降噪重构，假设观测信号是由原始信号叠加一个独立的高斯噪音 
，需要求 
 ，又说该等式是一个不适定问题(ill-posed problem)，不满足识定问题的三个条件，所以要用一个正则化的方法 
 第一部分是对原信号进行重构，第二部分是为了解决不适定问题的解不唯一，而且不同于机器学习的建模问题一样 
 是一个权值形式可以直接加权值 
 的L1-norm或L2-norm，这里的 
 是一个纯函数的形式，其导数的L2-norm最小化感觉上应该是保证了函数 
 不会产生太大的波动，这里可不可以跟防止过拟合联系起来，接下来说明这个约束会使 
 在频率域产生什么影响。
下面解这个式子，首先来看一下直接最小化泛函
 能不能解出来 
 有 
 ，然后根据E-L方程的引理 
 得 
 ，往下忘了怎么求了。。。。。想起来怎么求再继续写-_-!，这个直接求的方法与原文没什么关系的
由上面的傅立叶变换的微分性质知道，用傅立叶变换在复数域很方便可以把微分约掉，而且还有一个叫什么Plancherel傅立叶等距映射的东西，时域的L2范数与傅立叶变换到的频率域的L2范数等距，故上面的最小化的项可以直接用傅立叶变换转换到频域
 ，这里需要注意的是利用上面的那个什么等距定理有 
 和 
 的 
 是同一个 
 ，这一点很重要，第二个是L2范数大于0，则把其展开为泛函 
 求最极值，由于这里面只有 
 直接求偏导即可也不用解微分方程 
可以看到得到的
 相当于对观测信号在 
 在频率段进行滤波，过滤掉了高频部分，这说明加了该导数的L1正则约束与上面的直观感觉是一样的，过滤掉了高频部分，减弱 
 的波动
再往下进行，模态分解需要把原始信号分解成多个子信号的和，我们为了和原文对应修改一下符号表示，
 表示观测的采样信号， 
 表示分解得到的基函数，则上面的约束对象变为 
 同样先转化为频率域再求极值 
 泛函 
 每个基函数基于其他的基函数更新，相当于每个基函数是原信号剩余部分的低通滤波，每次迭代都是保留剩余信号的低频率部分。
到现在为止我们发现每个基函数都会趋向于每次的剩余信号分量的低频部分，这与我们原始的假设“每个基函数都有不同的频率分量”是相悖的，但根据上面的低通滤波的性质，每个基函数进行特定频率的滤波应该就能解决这个问题了，那么上面的式子就简单的变为
 其中， 
 为每个基函数 
 的中心频率，该式就是变分模态分解的基函数的更新公式，我们来看一下这个式子应该如何得到，以便于找到中心频率的更新公式
由上面的一步一步的演化发现，基函数对剩余信号的低通滤波是由导数的L2正则最小化带来的，要得到基函数的中心频率约束也要从这个地方入手，由上面的推导可知每个基函数都会被约束到
 频率附近，那么我们把基函数的频率增加各自的中心频率 
 得到 
 ，并保证 
 即可，则相当于对每个基函数乘以了一个 
 ，这里需要对频率进行变换，我们沿用文章开头的解析信号的性质，认为 
 是一个复数的解析信号，同样也需要把观测信号 
 预先转化为解析信号，下文默认都是解析信号，则每一个基函数转换频率后的导数变为
傅立叶变换得 
 先来看看傅立叶变换为什么会得到这个，而不是 
 反傅立叶变换 
 感觉应该不是这样证的，我数学不是很好，不怎么会这个变量变换
至此，约束变为
 傅立叶变换 
 根据论文原文把 
 然后求最小就可以得到上面 
 的更新公式 
最后，为了保证每个点处的重构信号与原信号尽可能相似，增加了每个点处的重构约束，其实这一项并不是必需的，最终的约束对象为
 ，然后拉格朗日乘子法带进去，但其需要满足下式才有意义 
 这个式子还是符合Parseval定理，故整理一下 
最后整理一下更新公式：
 其中 
 使用梯度下降更新
代码
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.signal import hilbert
T = 1000
fs = 1./T
t = np.linspace(0, 1, 1000,endpoint=True)
f_1 = 10
f_2 = 50
f_3 = 100
mode_1 = (2 * t) ** 2
mode_2 = np.sin(2 * np.pi * f_1 * t)
mode_3 = np.sin(2 * np.pi * f_2 * t)
mode_4 = np.sin(2 * np.pi * f_3 * t)
f = mode_1 + mode_2 + mode_3 + mode_4 + 0.5 * np.random.randn(1000)
plt.figure(figsize=(6,3), dpi=150)
plt.plot(f, linewidth=1)
class VMD:
def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):
""":param K: 模态数:param alpha: 每个模态初始中心约束强度:param tau: 对偶项的梯度下降学习率:param tol: 终止阈值:param maxIters: 最大迭代次数:param eps: eps"""
self.K =K
self.alpha = alpha
self.tau = tau
self.tol = tol
self.maxIters = maxIters
self.eps = eps
def __call__(self, f):
T = f.shape[0]
t = np.linspace(1, T, T) / T
omega = t - 1. / T
# 转换为解析信号
f = hilbert(f)
f_hat = np.fft.fft(f)
u_hat = np.zeros((self.K, T), dtype=np.complex)
omega_K = np.zeros((self.K,))
lambda_hat = np.zeros((T,), dtype=np.complex)
# 用以判断
u_hat_pre = np.zeros((self.K, T), dtype=np.complex)
u_D = self.tol + self.eps
# 迭代
n = 0
while n < self.maxIters and u_D > self.tol:
for k in range(self.K):
# u_hat
sum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]
res = f_hat - sum_u_hat
u_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)
# omega
u_hat_k_2 = np.abs(u_hat[k, :]) ** 2
omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)
# lambda_hat
sum_u_hat = np.sum(u_hat, axis=0)
res = f_hat - sum_u_hat
lambda_hat -= self.tau * res
n += 1
u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)
u_hat_pre[::] = u_hat[::]
# 重构，反傅立叶之后取实部
u = np.real(np.fft.ifft(u_hat, axis=-1))
omega_K = omega_K * T
idx = np.argsort(omega_K)
omega_K = omega_K[idx]
u = u[idx, :]
return u, omega_K
K = 4
alpha = 2000
tau = 1e-6
vmd = VMD(K, alpha, tau)
u, omega_K = vmd(f)
omega_K
# array([0.85049797, 10.08516203, 50.0835613, 100.13259275]))
plt.figure(figsize=(5,7), dpi=200)
plt.subplot(4,1,1)
plt.plot(mode_1, linewidth=0.5, linestyle='--')
plt.plot(u[0,:], linewidth=0.2, c='r')
plt.subplot(4,1,2)
plt.plot(mode_2, linewidth=0.5, linestyle='--')
plt.plot(u[1,:], linewidth=0.2, c='r')
plt.subplot(4,1,3)
plt.plot(mode_3, linewidth=0.5, linestyle='--')
plt.plot(u[2,:], linewidth=0.2, c='r')
plt.subplot(4,1,4)
plt.plot(mode_4, linewidth=0.5, linestyle='--')
plt.plot(u[3,:], linewidth=0.2, c='r')
# []
可以看到有比较强的端点效应，端点处会有重叠，文章原始论文中采用的方法是对称拼接的方法，把原信号复制一份然后拼成两半，一半对称放前面，一般对称放后面
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.signal import hilbert
T = 1000
fs = 1./T
t = np.linspace(0, 1, 1000,endpoint=True)
f_1 = 10
f_2 = 50
f_3 = 100
mode_1 = np.sin(2 * np.pi * f_1 * t)
mode_2 = np.sin(2 * np.pi * f_2 * t)
mode_3 = np.sin(2 * np.pi * f_3 * t)
f = np.concatenate((mode_1[:301], mode_2[301:701], mode_3[701:])) + 0.1 * np.random.randn(1000)
plt.figure(figsize=(6,3), dpi=150)
plt.plot(f, linewidth=0.5)
# []
class VMD:
def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):
""":param K: 模态数:param alpha: 每个模态初始中心约束强度:param tau: 对偶项的梯度下降学习率:param tol: 终止阈值:param maxIters: 最大迭代次数:param eps: eps"""
self.K =K
self.alpha = alpha
self.tau = tau
self.tol = tol
self.maxIters = maxIters
self.eps = eps
def __call__(self, f):
N = f.shape[0]
# 对称拼接
f = np.concatenate((f[:N//2][::-1], f, f[N//2:][::-1]))
T = f.shape[0]
t = np.linspace(1, T, T) / T
omega = t - 1. / T
# 转换为解析信号
f = hilbert(f)
f_hat = np.fft.fft(f)
u_hat = np.zeros((self.K, T), dtype=np.complex)
omega_K = np.zeros((self.K,))
lambda_hat = np.zeros((T,), dtype=np.complex)
# 用以判断
u_hat_pre = np.zeros((self.K, T), dtype=np.complex)
u_D = self.tol + self.eps
# 迭代
n = 0
while n < self.maxIters and u_D > self.tol:
for k in range(self.K):
# u_hat
sum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]
res = f_hat - sum_u_hat
u_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)
# omega
u_hat_k_2 = np.abs(u_hat[k, :]) ** 2
omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)
# lambda_hat
sum_u_hat = np.sum(u_hat, axis=0)
res = f_hat - sum_u_hat
lambda_hat -= self.tau * res
n += 1
u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)
u_hat_pre[::] = u_hat[::]
# 重构，反傅立叶之后取实部
u = np.real(np.fft.ifft(u_hat, axis=-1))
u = u[:, N//2 : N//2 + N]
omega_K = omega_K * T / 2
idx = np.argsort(omega_K)
omega_K = omega_K[idx]
u = u[idx, :]
return u, omega_K
K = 3
alpha = 2000
tau = 1e-6
vmd = VMD(K, alpha, tau)
u, omega_K = vmd(f)
omega_K
# array([  9.68579292,  50.05232833, 100.12321047])
plt.figure(figsize=(5,7), dpi=200)
plt.subplot(3,1,1)
plt.plot(mode_1, linewidth=0.5, linestyle='--')
plt.plot(u[0,:], linewidth=0.2, c='r')
plt.subplot(3,1,2)
plt.plot(mode_2, linewidth=0.5, linestyle='--')
plt.plot(u[1,:], linewidth=0.2, c='r')
plt.subplot(3,1,3)
plt.plot(mode_3, linewidth=0.5, linestyle='--')
plt.plot(u[2,:], linewidth=0.2, c='r')
# []
好像结果要好一点，VMD的一个缺点是K的值对结果有很大影响，但这个迭代过程，怎么开始怎么迭代怎么结束都可以自己控制，感觉可以按照自己的需求来定制怎么动态决定K

标签（空格分隔）： 数学
这一篇写一下变分模态分解（原始论文：Variational Mode Decomposition），跟原始论文思路思路一致但有一点点不太一样，原始论文写的很好，但我不是通信专业没有学过信号相关课程一开始看起来有点费劲。模态分解认为信号是由不同“模态”的子信号叠加而成的，而变分模态分解则认为信号是由不同频率占优的子信号叠加而成的，其目的是要把信号分解成不同频率的子信号。变分模态分解的分解结果如图所示
基础
一开始论文看不懂的原因是缺少相关前置知识，但一旦顺下来就会感觉其实没有那么难，难的是作者的思路很巧妙，先写下我遇到的这些知识盲点
第一点是傅立叶变换的微分性质， 
 的傅立叶变换为 
 ，其导数 
 的傅立叶变换为 
另一点是解析信号，现实世界只能采集实信号，但实信号有很多不好用的性质，如存在负频率，无法直接得到调制频率后的实信号等。 设原始信号是一个实信号 
为了方便表示 
 为随 
 变化的函数，相当于瞬时频率，解析信号是一个复信号，可以通过希尔伯特变换得到  
  解析信号的实部是原本的实数信号，并且经过调频之后复数信号的实部仍然是调频之后的实信号，如对信号 
 增加频率 
 ，只需要乘以 
 即可 
  由此，其实部相当于在原本频率 
 的基础上增加了频率 
 ，如下面的
matlab脚本
clear;close all;clc;
t = 1:0.01:10;
%%
f1 = sin(20*t).*(t-5).^2;
subplot(3,1,1);
plot(f1);
ylim([-25 25]);
%%
f2 = sin(50*t).*(t-5).^2;
subplot(3,1,2);
plot(f2);
ylim([-25 25]);
%%
H = hilbert(f1);
f_hat = H.*exp(1i*30.*t);
subplot(3,1,3);
plot(real(f_hat));
ylim([-25 25]);
matlab的hilbert函数包括希尔伯特变换和解析函数转换两部分，直接得到实信号的解析信号，其中希尔伯特变换   
正文
接下来我们看看如何一步一步得到变分模态分解的思路
原论文通过一个信号降噪问题进行说明，现需要对采样信号 
 进行降噪重构，假设观测信号是由原始信号叠加一个独立的高斯噪音 
，需要求 
 ，又说该等式是一个不适定问题（ill-posed problem），不满足识定问题的三个条件，所以要用一个正则化的方法 
 第一部分是对原信号进行重构，第二部分是为了解决不适定问题的解不唯一，而且不同于机器学习的建模问题一样 
 是一个权值形式可以直接加权值 
 的L1-norm或L2-norm，这里的 
 是一个纯函数的形式，其导数的L2-norm最小化感觉上应该是保证了函数 
 不会产生太大的波动，这里可不可以跟防止过拟合联系起来，接下来说明这个约束会使 
 在频率域产生什么影响。
下面解这个式子，首先来看一下直接最小化泛函 
 能不能解出来 
 有 
 ，然后根据E-L方程的引理 
 得 
 ，往下忘了怎么求了。。。。。想起来怎么求再继续写-_-!，这个直接求的方法与原文没什么关系的
由上面的傅立叶变换的微分性质知道，用傅立叶变换在复数域很方便可以把微分约掉，而且还有一个叫什么Plancherel傅立叶等距映射的东西，时域的L2范数与傅立叶变换到的频率域的L2范数等距，故上面的最小化的项可以直接用傅立叶变换转换到频域 
，这里需要注意的是利用上面的那个什么等距定理有 
和 
 的 
 是同一个 
 ，这一点很重要，第二个是L2范数大于0，则把其展开为泛函 
 求最极值，由于这里面只有 
 直接求偏导即可也不用解微分方程 
可以看到得到的 
 相当于对观测信号在 
 在频率段进行滤波，过滤掉了高频部分，这说明加了该导数的L1正则约束与上面的直观感觉是一样的，过滤掉了高频部分，减弱 
 的波动
再往下进行，模态分解需要把原始信号分解成多个子信号的和，我们为了和原文对应修改一下符号表示， 
 表示观测的采样信号， 
 表示分解得到的基函数，则上面的约束对象变为 
 同样先转化为频率域再求极值 
 泛函 
 每个基函数基于其他的基函数更新，相当于每个基函数是原信号剩余部分的低通滤波，每次迭代都是保留剩余信号的低频率部分。
到现在为止我们发现每个基函数都会趋向于每次的剩余信号分量的低频部分，这与我们原始的假设“每个基函数都有不同的频率分量”是相悖的，但根据上面的低通滤波的性质，每个基函数进行特定频率的滤波应该就能解决这个问题了，那么上面的式子就简单的变为 
 其中， 
 为每个基函数 
 的中心频率，该式就是变分模态分解的基函数的更新公式，我们来看一下这个式子应该如何得到，以便于找到中心频率的更新公式
由上面的一步一步的演化发现，基函数对剩余信号的低通滤波是由导数的L2正则最小化带来的，要得到基函数的中心频率约束也要从这个地方入手，由上面的推导可知每个基函数都会被约束到 
 频率附近，那么我们把基函数的频率增加各自的中心频率 
 得到 
 ，并保证 
 即可，则相当于对每个基函数乘以了一个 
 ，这里需要对频率进行变换，我们沿用文章开头的解析信号的性质，认为 
 是一个复数的解析信号，同样也需要把观测信号 
 预先转化为解析信号，下文默认都是解析信号，则每一个基函数转换频率后的导数变为 
傅立叶变换得 
 先来看看傅立叶变换为什么会得到这个，而不是 
 反傅立叶变换 
 感觉应该不是这样证的，我数学不是很好，不怎么会这个变量变换
至此，约束变为 
 傅立叶变换 
 根据论文原文把 
 然后求最小就可以得到上面 
 的更新公式 
最后，为了保证每个点处的重构信号与原信号尽可能相似，增加了每个点处的重构约束，其实这一项并不是必需的，最终的约束对象为 
 ，然后拉格朗日乘子法带进去，但其需要满足下式才有意义 
 这个式子还是符合Parseval定理，故整理一下 
最后整理一下更新公式： 
 其中 
 使用梯度下降更新
代码
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.signal import hilbert
T = 1000
fs = 1./T
t = np.linspace(0, 1, 1000,endpoint=True)
f_1 = 10
f_2 = 50
f_3 = 100
mode_1 = (2 * t) ** 2
mode_2 = np.sin(2 * np.pi * f_1 * t)
mode_3 = np.sin(2 * np.pi * f_2 * t)
mode_4 = np.sin(2 * np.pi * f_3 * t)
f = mode_1 + mode_2 + mode_3 + mode_4 + 0.5 * np.random.randn(1000)

plt.figure(figsize=(6,3), dpi=150)
plt.plot(f, linewidth=1)
class VMD:
    def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):
        """
        :param K: 模态数
        :param alpha: 每个模态初始中心约束强度
        :param tau: 对偶项的梯度下降学习率
        :param tol: 终止阈值
        :param maxIters: 最大迭代次数
        :param eps: eps
        """
        self.K =K
        self.alpha = alpha
        self.tau = tau
        self.tol = tol
        self.maxIters = maxIters
        self.eps = eps

    def __call__(self, f):
        T = f.shape[0]
        t = np.linspace(1, T, T) / T
        omega = t - 1. / T
        # 转换为解析信号
        f = hilbert(f)
        f_hat = np.fft.fft(f)
        u_hat = np.zeros((self.K, T), dtype=np.complex)
        omega_K = np.zeros((self.K,))
        lambda_hat = np.zeros((T,), dtype=np.complex)
        # 用以判断
        u_hat_pre = np.zeros((self.K, T), dtype=np.complex)
        u_D = self.tol + self.eps

        # 迭代
        n = 0
        while n < self.maxIters and u_D > self.tol:
            for k in range(self.K):
                # u_hat
                sum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]
                res = f_hat - sum_u_hat
                u_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)

                # omega
                u_hat_k_2 = np.abs(u_hat[k, :]) ** 2
                omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)

            # lambda_hat
            sum_u_hat = np.sum(u_hat, axis=0)
            res = f_hat - sum_u_hat
            lambda_hat -= self.tau * res

            n += 1
            u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)
            u_hat_pre[::] = u_hat[::]

        # 重构，反傅立叶之后取实部
        u = np.real(np.fft.ifft(u_hat, axis=-1))

        omega_K = omega_K * T
        idx = np.argsort(omega_K)
        omega_K = omega_K[idx]
        u = u[idx, :]
        return u, omega_K
K = 4
alpha = 2000
tau = 1e-6
vmd = VMD(K, alpha, tau)
u, omega_K = vmd(f)
omega_K
# array([0.85049797, 10.08516203, 50.0835613, 100.13259275]))
plt.figure(figsize=(5,7), dpi=200)
plt.subplot(4,1,1)
plt.plot(mode_1, linewidth=0.5, linestyle='--')
plt.plot(u[0,:], linewidth=0.2, c='r')

plt.subplot(4,1,2)
plt.plot(mode_2, linewidth=0.5, linestyle='--')
plt.plot(u[1,:], linewidth=0.2, c='r')

plt.subplot(4,1,3)
plt.plot(mode_3, linewidth=0.5, linestyle='--')
plt.plot(u[2,:], linewidth=0.2, c='r')

plt.subplot(4,1,4)
plt.plot(mode_4, linewidth=0.5, linestyle='--')
plt.plot(u[3,:], linewidth=0.2, c='r')
# [<matplotlib.lines.Line2D at 0x7fac532f4dd8>]
可以看到有比较强的端点效应，端点处会有重叠，文章原始论文中采用的方法是对称拼接的方法，把原信号复制一份然后拼成两半，一半对称放前面，一般对称放后面
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.signal import hilbert

T = 1000
fs = 1./T
t = np.linspace(0, 1, 1000,endpoint=True)

f_1 = 10
f_2 = 50
f_3 = 100
mode_1 = np.sin(2 * np.pi * f_1 * t)
mode_2 = np.sin(2 * np.pi * f_2 * t)
mode_3 = np.sin(2 * np.pi * f_3 * t)
f = np.concatenate((mode_1[:301], mode_2[301:701], mode_3[701:])) + 0.1 * np.random.randn(1000)

plt.figure(figsize=(6,3), dpi=150)
plt.plot(f, linewidth=0.5)
# [<matplotlib.lines.Line2D at 0x7fc2134b8630>]
class VMD:
    def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):
        """
        :param K: 模态数
        :param alpha: 每个模态初始中心约束强度
        :param tau: 对偶项的梯度下降学习率
        :param tol: 终止阈值
        :param maxIters: 最大迭代次数
        :param eps: eps
        """
        self.K =K
        self.alpha = alpha
        self.tau = tau
        self.tol = tol
        self.maxIters = maxIters
        self.eps = eps

    def __call__(self, f):
        N = f.shape[0]
        # 对称拼接
        f = np.concatenate((f[:N//2][::-1], f, f[N//2:][::-1]))
        T = f.shape[0]
        t = np.linspace(1, T, T) / T
        omega = t - 1. / T
        # 转换为解析信号
        f = hilbert(f)
        f_hat = np.fft.fft(f)
        u_hat = np.zeros((self.K, T), dtype=np.complex)
        omega_K = np.zeros((self.K,))
        lambda_hat = np.zeros((T,), dtype=np.complex)
        # 用以判断
        u_hat_pre = np.zeros((self.K, T), dtype=np.complex)
        u_D = self.tol + self.eps

        # 迭代
        n = 0
        while n < self.maxIters and u_D > self.tol:
            for k in range(self.K):
                # u_hat
                sum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]
                res = f_hat - sum_u_hat
                u_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)

                # omega
                u_hat_k_2 = np.abs(u_hat[k, :]) ** 2
                omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)

            # lambda_hat
            sum_u_hat = np.sum(u_hat, axis=0)
            res = f_hat - sum_u_hat
            lambda_hat -= self.tau * res

            n += 1
            u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)
            u_hat_pre[::] = u_hat[::]

        # 重构，反傅立叶之后取实部
        u = np.real(np.fft.ifft(u_hat, axis=-1))
        u = u[:, N//2 : N//2 + N]

        omega_K = omega_K * T / 2
        idx = np.argsort(omega_K)
        omega_K = omega_K[idx]
        u = u[idx, :]
        return u, omega_K
K = 3
alpha = 2000
tau = 1e-6
vmd = VMD(K, alpha, tau)

u, omega_K = vmd(f)
omega_K
# array([  9.68579292,  50.05232833, 100.12321047])

plt.figure(figsize=(5,7), dpi=200)
plt.subplot(3,1,1)
plt.plot(mode_1, linewidth=0.5, linestyle='--')
plt.plot(u[0,:], linewidth=0.2, c='r')

plt.subplot(3,1,2)
plt.plot(mode_2, linewidth=0.5, linestyle='--')
plt.plot(u[1,:], linewidth=0.2, c='r')

plt.subplot(3,1,3)
plt.plot(mode_3, linewidth=0.5, linestyle='--')
plt.plot(u[2,:], linewidth=0.2, c='r')
# [<matplotlib.lines.Line2D at 0x7fc2134075c0>]
好像结果要好一点，VMD的一个缺点是K的值对结果有很大影响，但这个迭代过程，怎么开始怎么迭代怎么结束都可以自己控制，感觉可以按照自己的需求来定制怎么动态决定K

浅谈VMD(变分模态分解)
统一解释一下，我也是才学习这个，很多地方不懂。至于译文，也是自己根据原论文进行的一些翻译，很多地方不是很准确，研一的可以自己试着去读英文原文，而且大家学校应该都能下载这篇论文。研二的如果急着要用可以参考借鉴一下。

新的自适应信号处理方法，对非平稳、非线性信号具有良好处理效果。
文章说明
因为最近论文需要用到VMD，所以看了几篇关于VMD的论文，VMD在2014年提出，所以其论文比较新，而且知网上的论文都是基于VMD的应用，里面的原理都是从VMD原文摘抄，不太好理解，我也看了VMD这篇论文，我也进行了一定的翻译，译文在后面，但翻译总有差距，希望英语好的专业好的能提出建议，谢谢。
阅读说明
我知道好多人看着VMD看博客最想知道的就是这东西的应用和大概步骤原理，而具体原理算法不太感兴趣，而且也不太容易看懂。本文既然是浅谈，就从VMD的步骤和应用。
VMD工作原理步骤
步骤

VMD是通过迭代搜寻变分模型（具体怎么搜寻，请亲们自己看，我主要讲他的大概）最优解， 来确定我们所知的模态uk(t)及其对应的中心频率ωk和带宽。

每个模态都是具有中心频率的有限带宽（就是在频域中有在一定的宽度）。所有模态之和为源信号。
而对求最优解采用二次惩罚和拉格朗日乘数将上诉约束问题转换为非约束问题，并用交替方向乘子法求解这个非约束问题， 通过迭代更新最终得到信号分解的所有模态。分解的所有模态中有包含主要信号的模态和包含噪声的模态。将包含主要信号的模态进行重构，从而达到去噪的效果。
代码步骤思路（uk和ωk更新算法）
1、初始化uk、ωk、λ和n=0，k=0
2、n=n+1（迭代次数）
3、k=k+1，根据VMD算法公式更新uk、ωk
4、又根据相关的算法更新拉格朗日乘数λ
5、知直到满足一定条件，停止迭代，不然转到2步骤
以上只是求每一个模态的单步骤


我的理解总步骤：
1、初始化uk、ωk、λ和n=0，
2、n=n+1（迭代次数）
3、根据VMD算法公式更新uk、ωk
4、又根据相关的算法更新拉格朗日乘数λ
5、知直到满足一定条件根据（相似系数来判断），停止迭代，不然转到2步骤
6、k=k+1，将源信号减去分解出来的模态，并作为下次一循环的源信号，转到步骤1

如何判断相关模态
判断

用信号与模态的相似程度来判断信号与噪声 。推荐一篇论文，他对VMD进行了一些优化。例如：在VMD中一般采用局部重构，即将与原信号相似的模态就认为是信号，与原信号相差大的模态认为噪声，然而噪声模态中其实还含有一些信号，用一定方法提取信号，可增加信噪比和可信度。同理（我自己的看法），采用定的滤波器处理信号模态来去除其中的噪声会不会提高信噪比。这是一个方向。
推荐论文：基于VMD的激光雷达回波信号去噪方法研究
应用
其实看了上面的原理，其实大家都该有一个基本的想法，VMD将信号分解与重构，用的最多就是来去噪。

我准备用来进行心率提取

其他应用：谐波提取、用的最多的滚动轴承检测、故障检测等等
VMD原文
下面是原文链接，可以自行下载，也可以联系我，还有我自己的译文，翻译一班不是太准确，但能帮助你更好的理解。
缺点及解决方法
1、最大的局限性是边界效应和突发的信号。这与基于L2平滑阶段的使用密切相关，该阶段过度惩罚了域边界和内部的跳跃。
2、长期模态的光谱带会随着时间的推移而急剧变化，并且会在全局范围内重叠。在这里，直接的解决方案是将信号分解成较短的块，在这些块上信号足够稳定。
3、要求预先定义模态数K。这是我们与许多成功的聚类和分段算法（例如k-means）共享的缺点。
补充
在写论文，通常拿VMD与EMD(经验模态分解)、EWT（经验小波变换）进行性能对比。

对于预定义模态数K的解决方法，推荐一篇论文：基于反馈变分模式分解的单通道盲源分离算法，他采用反馈模式对K自动识别定义，解除了预定义的先验模式。
EMD以及VMD的matlab代码（求大家给个赞别）
VMD链接: https://pan.baidu.com/s/1tlPoq07HTsdhr1zEN7bzdw

提取码: wac3

EMD链接：https://pan.baidu.com/s/1sAvwEV8bOCannCjVBJNcpg

提取码：xww7

	

本文内容来源于《测绘通报》2020年第8期，审图号：GS(2020)3611号
变分模态分解及能量熵在地心运动降噪中的应用  
王庆余, 杜宁, 王莉, 张小东, 吴磊, 熬逍                                                                                                                            
贵州大学矿业学院, 贵州 贵阳 550025
摘要：
针对地心运动时间序列噪声种类复杂，随机性强，信号与噪声难以有效分离等问题，本文采用网平移法对IGS站周解进行解算，得到2012-2018年的地心运动时间序列，并提出了一种基于变分模态分解(VMD)及能量熵的地心运动时间序列降噪方法。首先，对各方向时间序列进行VMD分解，获得各方向高频依次到低频的时间序列分量；然后，计算每个变分模态分量的能量熵，辨识出噪声与信号的分界，并将信号分量进行重构，得到降噪后的地心运动时间序列；最后，通过与基于EMD和EEMD的降噪方法对比，从相关系数、信噪比、剩余能量百分比、方差贡献率等参数评价指标上定量说明该方法对地心运动时间序列降噪表现出更好的降噪效果。
关键词网平移法  降噪  变分模态分解  能量熵  信噪比 引文格式：王庆余, 杜宁, 王莉, 张小东, 吴磊, 熬逍. 变分模态分解及能量熵在地心运动降噪中的应用[J]. 测绘通报, 2020, 0(8): 59-64.基金项目：贵州省科技计划(黔科合基础〔2017〕1026)参考文献[1] 程鹏飞,成英燕.我国毫米级框架实现与维持发展现状和趋势[J].测绘学报,2017,46(10):1327-1335.
[2] 乔灵娜,赵春梅,马天明.基于EMD方法的地心运动时间序列分析[J].测绘通报,2019(2):6-11.
[3] 魏冠军,陈晨,黄逸宇,等.基于小波变换的CORS高程坐标时序周期性分析与降噪方法研究[J].全球定位系统,2019,44(4):61-67.
[4] 张双成,何月帆,李振宇,等.EMD用于GPS时间序列降噪分析[J].大地测量与地球动力学,2017,37(12):1248-1252.
[5] 戴建彪,岳东杰,汤同旭,等.桥梁索塔GPS监测信号的小波分析[J].测绘科学,2019,44(9):1-11.
[6] 张恒璟,龙安森,文汉江.EEMDAN的CORS站高程时间序列分析方法[J].测绘科学,2020,45(2):33-38.
[7] 张成龙,刘超,赵兴旺,等.基于改进EEMD的GPS多路径误差建模与削弱方法研究[J].大地测量与地球动力学,2018,38(10):1021-1026.
[8] DRAGOMIRETSKIY K,ZOSSO D.Variational mode deco-mposition[J].IEEE Transactions on Signal Processing,2014,62(3):531-544.
[9] 龙丹,牛聪,周怀来,等.基于VMD算法在地震数据时频分析中的应用[J].地球物理学进展,2020,35(1):166-173.
[10] 任刚,贾继德,梅检民,等.基于变分模态分解和去趋势波动分析的柴油机振动信号去噪方法[J].内燃机工程,2019,40(2):76-81.
[11] 张杏莉,卢新明,贾瑞生,等.基于变分模态分解及能量熵的微震信号降噪方法[J].煤炭学报,2018,43(2):356-363.
[12] 郭海荣,杨元喜,焦文海.地心运动时间序列的抗差谱分析[J].测绘学报,2003,32(4):308-312.
[13] 郭金运,常晓涛,韩延本,等.由SLR观测的地心周期性运动(1993-2006年)[J].测绘学报,2009,38(4):311-317.
[14] 秦显平,杨元喜.用SLR数据导出的地心运动结果[J].测绘学报,2003,32(2):120-124.
[15] 贾瑞生,赵同彬,孙红梅,等.基于经验模态分解及独立成分分析的微震信号降噪方法[J].地球物理学报,2015,58(3):1013-1023.
[16] 李成武,解北京,杨威,等.基于HHT法的煤冲击破坏SHPB测试信号去噪[J].煤炭学报,2012,37(11):1796-1802.
作 者 简 介： 王庆余(1995-),男,硕士生,主要从事空间大地测量学理论的学习和研究。E-mail:qywang17@163.com
通 信 作 者：杜宁。E-mail:ndu1@gzu.edu.cn											E-mail:ndu1@gzu.edu.cn
                                   
《测绘通报》杂志社
编辑部电话：010-68531192(金老师)
E-mail：chtb@chinajournal.net.cn
往期精彩推荐
结合图像增强和恢复的增强现实跟踪注册
多源点云数据融合的建筑物精细化建模
GEE云平台支持下的西天山森林遥感监测与时空变化分析
如需了解更多文章精彩内容敬请订阅《测绘通报》！邮发代号：2-223