综述
作者：aresmiki
链接：https://www.zhihu.com/question/23701194/answer/167005497
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。
非平稳信号处理应该是现在信号处理技术最新的也是最热的研究方向了，信号处理方法从最早的时域统计到傅里叶变换的频域分析，是人们认识信号本质的一次巨大飞跃，给了分析人员换个角度看世界的方法，这个时期傅里叶变换在调和分析，谐波分析等领域得到的巨大发展，当然工程人员将傅里叶变换也迅速运用的工程运用中，得到了巨大发展，特别是在机械故障诊断和 功能性磁共振成像fMRI 中。但是现实却不是理论那么美好，人们发现，傅里叶变换作为一个全局变换，天然的少了另一个维度，如果将时间域信号比作一个平面中的物体的话，那么频域信号也同样是一个平面中的物体，只是给我们换了一个角度而已，而人们总是对三维世界的物体更具有直观了解，平面的东西始终是不具体不形象，信号一样，工程人员总是想知道信号有哪些频率，且这些频率在何时产生，而这个需求就给分析方法提出了一个要求，必须多一个维度，也就是频域变换需要保留时间信息，聪明的Gabor第一次提出了窗口傅里叶变换，这个信号分析带来了时间尺度，也第一次有了时频分析的概念，当然其理论就不介绍了，而其优点当然就是同时给了我们时间和频率的信息。科学的魅力就是，永远没有完美的理论，窗口福利叶变换的窗口如何选择成为了一个难题，选的太大，时间分辨率太大，选的太小，频率分辨率太大，如何能在该大的时候大，该小的时候小呢？伟大的Morlet告诉我们，傅里叶变换本身就是全局的，而我们为什么还非得坚守傅里叶变换这个全局框架下呢？这种变化快的信号本身就是短时局部的信号不断变换，为什么不是一个个局部的去分析信号呢，这也是小波变换的基本思想，小波将时频分析推向了研究高潮，这一时期小波理论不断发展，出现了许多小波，db系列小波，morlet小波，coif系列小波等等，可以看出来，小波研究都是基于小波基函数的研究，它不像傅里叶变换一样，基是固定的，而小波基函数有很多形式。这种灵活性给了小波广泛的运用优势，但是，科学就是这样，一种方法的诞生必然伴随着缺陷的诞生，小波基的这种灵活性却给分析人员带来困扰，如何针对不同信号选择这些基呢？是否有一个通用基来处理信号分解任务呢？还有一个关键问题是，小波变换受到测不准原理的限制，不能无限制的细分时间和频率。当然研究人员根据实际处理任务提出了相对通用的小波基，如图像处理中的曲波变换和脊波变换，这一时期小波研究含有另一个方向，不是一个小波不能完全符合处理任务吗？不是不好选择小波基吗？研究人员提出把多个小波基一起用，必然可以得到更好的结果，这典型的代表是多小波的发展。当然普通小波变换存在时移缺陷和频率混叠的缺陷，而双树复数小波DTCWT得提出基本解决了该问题。也有针对小波基难以构造的缺陷，提出了提升框架的小波变换，但是小波变换就是小波变换，本质还是没有改变，小波具有的问题，后面发展出来的方法依然存在这些问题，小波处理方法没有过时，但是其诸多问题却在工程运用中越来越明显，工程运用中信号千变万化，小波基选择始终是摆在分析人员面前的问题，还有就是楼主提到的消噪，小波消噪软阈值和硬阈值的选择会对信号消噪能力有巨大影响，选择不好要么消除不够，有么太多。更直白的说，你说你消除了噪声，那么万一别人问你，你凭什么说你消除的是噪声而不是有用信号了，这个问题是一个无法回答的问题。当然后小波时代崛起了一批数据驱动的非线性非平稳信号处理方法，这些方法的提出就是摆脱小波基选择困难，测不准原理限制的问题，这类方法典型代表就是EMD，EMD方法根据信号自身特点不断分离开来，噪声和有用信号安照不同频带划分分离开，当然由于EMD方法本身诸多问题，后来发展出来了EEMD，CEEMD，CEEMDAN。CWEEMDAN等方法，但是其数学基础问题依然存疑，可是，EMD的提出给了人们一个启发，信号分解就是要让信号根据自己特点分开而不是人工的参与，一切应该是自驱动自适应才是最好的，受该思想的启发，小波研究者也希望将小波推广到这类方法中，其典型代表SWT和EWT，虽然方法很妙，但是其本质问题却并非数据驱动，有兴趣可以自行查阅文献，当然继承该思想的还有VMD、NSP、ITD方法，不得不提一下，这些方法都是非常适合非平稳非线性信号处理，但是分解中的一些参数设置却十分困难。关于信号分析最前沿的方法和理论，推荐专栏 信号处理与机器学习 - 知乎专栏 里面讲了很多信号处理方法的核心思想和最新方法。。。。
VMD（变分模态分解）（https://zhuanlan.zhihu.com/p/66898788?utm_source=qq）
这一篇写一下变分模态分解（原始论文：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函数包括希尔伯特变换和解析函数转换两部分，直接得到实信号的解析信号，其中希尔伯特变换
正文
接下来我们看看如何一步一步得到变分模态分解的思路
原论文通过一个信号降噪问题进行说明，现需要对采样信号
下面解这个式子，首先来看一下直接最小化泛函
由上面的傅立叶变换的微分性质知道，用傅立叶变换在复数域很方便可以把微分约掉，而且还有一个叫什么Plancherel傅立叶等距映射的东西，时域的L2范数与傅立叶变换到的频率域的L2范数等距，故上面的最小化的项可以直接用傅立叶变换转换到频域
可以看到得到的
再往下进行，模态分解需要把原始信号分解成多个子信号的和，我们为了和原文对应修改一下符号表示，
到现在为止我们发现每个基函数都会趋向于每次的剩余信号分量的低频部分，这与我们原始的假设“每个基函数都有不同的频率分量”是相悖的，但根据上面的低通滤波的性质，每个基函数进行特定频率的滤波应该就能解决这个问题了，那么上面的式子就简单的变为
由上面的一步一步的演化发现，基函数对剩余信号的低通滤波是由导数的L2正则最小化带来的，要得到基函数的中心频率约束也要从这个地方入手，由上面的推导可知每个基函数都会被约束到
至此，约束变为
最后，为了保证每个点处的重构信号与原信号尽可能相似，增加了每个点处的重构约束，其实这一项并不是必需的，最终的约束对象为
最后整理一下更新公式：
代码
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from scipy.signal import hilbertT = 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 =Kself.alpha = alphaself.tau = tauself.tol = tolself.maxIters = maxItersself.eps = epsdef __call__(self, f):T = f.shape[0]t = np.linspace(1, T, T) / Tomega = 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 = 0while n < self.maxIters and u_D > self.tol:for k in range(self.K):# u_hatsum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]res = f_hat - sum_u_hatu_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)# omegau_hat_k_2 = np.abs(u_hat[k, :]) ** 2omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)# lambda_hatsum_u_hat = np.sum(u_hat, axis=0)res = f_hat - sum_u_hatlambda_hat -= self.tau * resn += 1u_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 * Tidx = np.argsort(omega_K)omega_K = omega_K[idx]u = u[idx, :]return u, omega_KK = 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 hilbertT = 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 =Kself.alpha = alphaself.tau = tauself.tol = tolself.maxIters = maxItersself.eps = epsdef __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) / Tomega = 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 = 0while n < self.maxIters and u_D > self.tol:for k in range(self.K):# u_hatsum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]res = f_hat - sum_u_hatu_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)# omegau_hat_k_2 = np.abs(u_hat[k, :]) ** 2omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)# lambda_hatsum_u_hat = np.sum(u_hat, axis=0)res = f_hat - sum_u_hatlambda_hat -= self.tau * resn += 1u_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 / 2idx = np.argsort(omega_K)omega_K = omega_K[idx]u = u[idx, :]return u, omega_KK = 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
