# 写在前面

介绍变分模态分解(Variational Mode Decomposition, VMD)[1]。

前几年大家都在做经验模态分解(EMD)的时候，我找到一个极度符合我胃口的变分模型，当时看得一知半解（功力不够），反复看了几遍都有新的认识，这次做个记录，后面会持续跟进。

# 变分模态分解

基于以下假设构造变分模型

  • 通过Hilbert变换得到每个模态的解析信号，从而得到单边频谱
  • 通过与调谐到各自估计的中心频率的指数混合，将模态的频谱移至基带
  • 通过解调信号的H1H^1高斯光滑性，即梯度的L2L^2范数的平方，来估计带宽

min⁡uk,ωk∑k∥∂t[(δ(t)+jπt)∗uk(t)]e−jωkt∥22 s.t. ∑kuk=f\min _{u_k,\omega_k}\sum_{k}\left\|\partial_t\left[\left(\delta(t)+\frac{j}{\pi t}\right) * u_{k}(t)\right] e^{-j \omega_{k} t}\right\|_{2}^{2} \text { s.t. } \sum_{k} u_{k}=f

下面将剖解这一模型

# 固有模态

调频调幅信号(AM-FM)

uk(t)=Ak(t)cos⁡(ϕk(t))u_k (t) = A_k (t) \cos (\phi_k (t))

# Hilbert变换

其本质是具有平移不变性的线性算子H\mathcal H，将余弦函数映射至正弦函数。从信号角度看，Hilbert变换是一种全通滤波器(all-pass filter)，可由转递函数刻画：

h^(w)=−jsgn(w)=−jww\hat h (w) = -j \text{sgn} (w) = -j\frac{w}{w}

对应的脉冲相应为

h(t)=1πth(t) = \frac{1}{\pi t}

与h(w)h(w)做卷积是不可积的，信号f(t)f(t)的Hilbert变换Hf(t)\mathcal H f(t) 通过卷积积分的柯西主值(Cauchy principal value )得到

Hf(t)=1πp.v.∫Rf(v)t−vdv{\mathcal H} f(t) = \frac{1}{\pi} \text{p.v.} \int_{\mathbb R} \frac{f(v)}{t-v} dv

Hilbert变换的主要用途是将实信号转化为解析信号。设f(t)f(t)是一个纯实信号，复值解析信号定义如下：

fA(t)=f(t)+jHf(t)=A(t)ejϕ(t)f_A(t) = f(t) +j {\mathcal H} f(t) = A(t) e^{j \phi(t)}

其中ejϕ(t)e^{j \phi(t)}描述复信号在时间上旋转的相量，ϕ(t)\phi(t)表示相位，振幅A(t)A(t)描述了实包络。

  • 解析信号直接继承相同的振幅信号

uk,A=Ak(t)(cos⁡(ϕ(t))+(jsin⁡(ϕ(t))))=Ak(t)ejϕ(t)u_{k,A} = A_k(t)(\cos (\phi(t)) + (j \sin(\phi(t)))) = A_k(t) e^{j \phi(t)}

  • 解析信号的单边频谱由非负频率构成

f(t)=ℜ{fA(t)}f(t) = \Re \{f_A(t)\}

  • 解析信号得Fourier变换

u^A(w)=∫−∞∞uA(t)e−jwtdt=(1+sgn(w))u^(w)\hat u_A (w) = \int_{-\infty}^{\infty} u_A(t) e^{-jwt}dt = (1+\text{sgn}(w))\hat u(w)

# 混频和外差解调

混合两个实信号会产生两组混合频率

2cos⁡(w1t)cos⁡(w2t)=cos⁡((w1+w2)t)+cos⁡((w1−w2)t)2\cos(w_1t)\cos(w_2t) = \cos((w_1+w_2)t)+\cos((w_1-w_2)t)

混合两个解析信号后频率仅由单个频率组成

ejw1tejw2t=ej(w1+w2)te^{j w_1 t} e^{j w_2 t} = e^{j (w_1 + w_2) t}

运用Fourier变换

fA(t)e−jω0t⟷Ff^A(ω)∗δ(ω+ω0)=f^A(ω+ω0)f_{A}(t) e^{-j \omega_{0} t} \stackrel{\mathcal{F}}{\longleftrightarrow} \hat{f}_{A}(\omega) * \delta\left(\omega+\omega_{0}\right)=\hat{f}_{A}\left(\omega+\omega_{0}\right)

因此，将解析信号与纯指数相乘会导致简单的频移。

# Wiener滤波

原始信号f(t)f(t)受加性零均值高斯噪声的观察信号f0(t)f_0(t)，即

f0=f+ηf_0 = f + \eta

信号的去噪是一个典型病态的可逆问题，通常用Tikhonov正则化来处理

min⁡f∥f−f0∥22+α∥∂tf∥22\min_f \|f - f_0\|_2^2 + \alpha \|\partial_t f\|_2^2

求解其Euler-Lagrange方程可得到最优解，在Fourier域上

f^(w)=f0^1+αw2\hat f(w) = \frac{\hat {f_0}}{1+\alpha w^2}

其中f^\hat f是信号ff做Fourier变换所得。通常恢复的信号是观察信号在频率w=0w=0附近的低通窄带选择，即Wiener滤波可视为低通滤波器，α\alpha表示白噪声的方差，信号具有低通w−2w^{-2}频谱先验。

# 模型求解

VMD模型包含线性约束条件，通常做法是转化为无约束优化问题。

  • 构造增广拉格朗日函数

L(uk,ωk,λ)=α∑k∥∂t((δ(t)+jπt)∗uk(t))e−jωkt∥22+∥f(t)−∑kuk(t)∥22+⟨λ(t),f(t)−∑kuk(t)⟩\begin{aligned} \mathcal{L} ( u_k, \omega_{k}, \lambda ) &= \alpha \sum_{k}\left\|\partial_t((\delta(t)+\frac{j}{\pi t}) * u_{k}(t)) e^{-j \omega_{k} t}\right\|_{2}^{2}\\ &+\Bigg\Vert f(t)-\sum_{k}u_k(t)\Bigg\Vert_2^2+ \Bigg\langle\lambda(t),f(t)-\sum_{k}u_k(t)\Bigg\rangle \end{aligned}

下面使用ADMM算法来寻找增广拉格朗日函数的鞍点。

  • 更新uku_k

固定$w_k,\lambda $后，对应子问题为

min⁡ukα∑k∥∂t[(δ(t)+jπt)∗uk(t)]e−jωkt∥22+∥f(t)−∑iui(t)+λ(t)2∥22\min_ {u_k}\alpha \sum_{k}\left\|\partial_{t}\left[\left(\delta(t)+\frac{j}{\pi t}\right) * u_{k}(t)\right] e^{-j \omega_{k} t}\right\|_{2}^{2}+\Bigg\Vert f(t)-\sum_{i}u_i(t) + \frac{\lambda(t)}{2}\Bigg\Vert_2^2

利用L2L^2范数意义下的Fourier等距性，该子问题可在谱域求解

min⁡u^kukα∥jw[(1+sgn(w+wk))u^k(w+wk)]∥22+∥f^(t)−∑iu^i(t)+λ^(t)2∥22\min_ {\hat u_k u_k}\alpha \left\|jw[(1+\text{sgn}(w+w_k))\hat u_k (w+w_k)]\right\|_{2}^{2}+\Bigg\Vert \hat f(t)-\sum_{i}\hat u_i(t) + \frac{\hat\lambda(t)}{2}\Bigg\Vert_2^2

令w←w−wkw \leftarrow w-w_k

min⁡u^kukα∥j(w−wk))[(1+sgn(w))u^k(w)]∥22+∥f^(t)−∑iu^i(t)+λ^(t)2∥22\min_ {\hat u_k u_k}\alpha \left\|j(w-w_k))[(1+\text{sgn}(w))\hat u_k (w)]\right\|_{2}^{2}+\Bigg\Vert \hat f(t)-\sum_{i}\hat u_i(t) + \frac{\hat\lambda(t)}{2}\Bigg\Vert_2^2

重构保真项利用实信号的Hermitian对称性，可将两项写为非负频率上的半空间积分：

min⁡u^kuk∫0∞4α(w−wk)2∣u^k(w)∣2+2∣f^(t)−∑iu^i(t)+λ^(t)2∣2dw\min_ {\hat u_k u_k} \int_0^{\infty} 4\alpha (w-w_k)^2|\hat u_k (w)|^{2}+ 2\left|\hat f(t)-\sum_{i}\hat u_i(t) + \frac{\hat\lambda(t)}{2}\right|^2 dw

对于该二次优化问题，在正频率积分下，令一阶导数为零，得到

u^kn+1(ω)=x^(ω)−∑i≠ku^i(ω)+λ^(ω)21+2α(ω−ωk)2.\begin{aligned} \hat{u}_k^{n+1}(\omega) = \frac{\hat{x}(\omega)-\sum_{i\neq k}\hat{u}_i(\omega)+\frac{\hat{\lambda}(\omega)}{2}}{1+2\alpha(\omega-\omega_k)^2}. \end{aligned}

该式是当前残差对应先验(w−wk)−2(w-w_k)^{-2}的Wiener滤波形式。对u^k\hat u_k进行Fourier逆变换可得到解析信号的实部分uku_k。

  • 更新wkw_k

中心频率仅出现在带宽先验项，因此固定uk,λu_k,\lambda后子问题如下

min⁡wk∥∂t[(δ(t)+jπt)∗uk(t)]e−jωkt∥22\min_ {w_k} \left\|\partial_{t}\left[\left(\delta(t)+\frac{j}{\pi t}\right) * u_{k}(t)\right] e^{-j \omega_{k} t}\right\|_{2}^{2}

同样转换到Fourier域

min⁡w^k∫0∞(w−wk)2∣u^k(w)∣2dw\min_ {\hat w_k} \int_0^{\infty} (w-w_k)^2|\hat u_k (w)|^{2} dw

对于该二次优化问题，同样令一阶导数为零可得最优解

ωkn+1=∫0∞ω∣u^k(ω)∣2dω∫0∞∣u^k(ω)∣2dω,\omega_k^{n+1}=\frac{\int^{\infty}_0 \omega|\hat{u}_k(\omega)|^2 d\omega}{\int^{\infty}_0 |\hat{u}_k(\omega)|^2 d\omega},

该平均载频是模态所观察到的瞬时相位的最小二乘线性回归的频率。

  • 更新λ\lambda

拉格朗日乘子的更新通过对偶上升实现

λ^n+1(w)=λ^n(w)+τ(f^(w)−∑ku^kn+1(w))\hat\lambda ^{n+1} (w) = \hat\lambda ^{n} (w) + \tau (\hat f(w)-\sum_{k}\hat u_k^{n+1}(w))

其中参数τ\tau权衡了正则项和保真项，若需要精确恢复信号，即不存在高斯噪声情况下，设置τ=0\tau = 0。

# 算法流程

# 多元变分模态分析

下面介绍多元模型[2]。对于多元振荡，考虑CC个实值AM-FM信号{ui(t)}i=1C\{u_i(t)\}_{i=1}^C，用向量形式表示

u(t)=[u1(t)u2(t)⋮uC(t)]=[a1(t)cos⁡(ϕ1(t))a2(t)cos⁡(ϕ2(t))⋮aC(t)cos⁡(ϕC(t))]\mathbf{u}(t)=\left[\begin{array}{c} u_{1}(t) \\ u_{2}(t) \\ \vdots \\ u_{C}(t) \end{array}\right]=\left[\begin{array}{c} a_{1}(t) \cos \left(\phi_{1}(t)\right) \\ a_{2}(t) \cos \left(\phi_{2}(t)\right) \\ \vdots \\ a_{C}(t) \cos \left(\phi_{C}(t)\right) \end{array}\right]

对该向量信号u(t)\mathbf{u}(t)进行Hilbert变换，即对各个分量进行Hilbert变换，产生解析信号向量

uA(t)=u(t)+jHu(t)=[u1,A(t)u2,A(t)⋮uC,A(t)]=[a1(t)eϕ1(t)a2(t)eϕ2(t)⋮aC(t)eϕC(t)]\mathbf{u}_A(t)= \mathbf{u}(t) + j{\mathcal H} \mathbf{u}(t) =\left[\begin{array}{c} u_{1,A}(t) \\ u_{2,A}(t) \\ \vdots \\ u_{C,A}(t) \end{array}\right]=\left[\begin{array}{c} a_{1}(t) e^{\phi_{1}(t)} \\ a_{2}(t) e^{\phi_{2}(t)} \\ \vdots \\ a_{C}(t) e^{\phi_{C}(t)} \end{array}\right]

反之，取解析信号向量的实部即可恢复原信号

u(t)=ℜ{uA(t)}\mathbf{u}(t) = \Re \{\mathbf{u}_A(t)\}

# 模型构建

  • 目的：从输入信号x(t)=[x1(t),x2(t),…,xC(t)]\mathbf{x}(t) = [x_1(t),x_2(t),\dots,x_C(t)]提取KK个多元模态uk(t)=[u1(t),u2(t),…,uC(t)]\mathbf {u}_k (t) = [u_1(t),u_2(t),\dots,u_C(t)]
  • 先验假设：
    • 提取模态的带宽之和最小
    • 提取模态之和能精确恢复原信号
  • 通过采用谐波位移的梯度函数的L2范数来估算模态带宽，对应的目标函数为

∑k∥∂t[uk,A(t)e−jωkt]∥22\sum_{k}\left\|\partial_{t}\left[ \mathbf{u}_{k,A}(t) e^{-j \omega_{k} t}\right] \right\|_{2}^{2}

上式在整个信号向量uk,A(t)\mathbf{u}_{k,A}(t)的谐波混合中使用单个频率成分wkw_k，多元实信号uk(t)\mathbf{u}_{k}(t)的所有通道均具有频率wkw_k。因此，通过将uk,A(t)\mathbf{u}_{k,A}(t)所有通道的单边频谱偏移wkw_k并利用Frobenius范数来估算调制多元振荡的带宽。其目标函数按分量表示为

∑k∑c∥∂t[uk,c,A(t)e−jωkt]∥22\sum_{k}\sum_{c}\left\|\partial_t\left[ u_{k,c,A}(t) e^{-j \omega_{k} t}\right] \right\|_{2}^{2}

其中uk,c,A(t)u_{k,c,A}(t)表示通道cc模态kk的解析信号。

# 多元变分模态模型

min⁡uk,c,wk∑k∑c∥∂t[uk,c,A(t)e−jωkt]∥22 s.t. ∑kuk,c(t)=xc,∀c\min_{ u_{k,c}, w_k } \sum_{k}\sum_{c}\left\|\partial_{t}\left[ u_{k,c,A}(t) e^{-j \omega_{k} t}\right] \right\|_{2}^{2} \text { s.t. } \sum_k u_{k,c}(t) = x_c, \forall c

# 模型求解

与VMD模型一样，首先转化约束优化问题至无约束优化问题，再使用ADMM算法求解

  • 增广的拉格朗日函数

L(uk,c,ωk,λc)=α∑k∑c∥∂t[uk,c,A(t)e−jωkt]∥22+∑c∥xc(t)−∑kuk,c(t)∥22+∑c⟨λc(t),xc(t)−∑kuk,c(t)⟩\begin{aligned} \mathcal{L}({u_{k,c},\omega_{k},\lambda_c}) &= \alpha\ {\sum_k\sum_c}\Bigg\Vert\partial_t\Big[u_{k,c,A}(t) e^{-j\omega_kt}\Big]\Bigg\Vert^2_2\\ &+{\sum_c}\Bigg\Vert x_c(t)-\sum_{k}u_{k,c}(t)\Bigg\Vert_2^2+ {\sum_c}\Bigg\langle\lambda_c(t),x_c(t)-\sum_{k}u_{k,c}(t)\Bigg\rangle \end{aligned}

  • 更新模态

固定{ωk},λc\{\omega_{k}\},\lambda_c，可得如下子问题

min⁡uk,cα∥∂t[uk,c,A(t)e−jωkt]∥22+∥xc(t)−∑iui,c(t)+λc(t)2∥22\min_{u_{k,c}} \alpha\Bigg\Vert\partial_t\Big[u_{k,c,A}(t) e^{-j\omega_kt}\Big]\Bigg\Vert^2_2+ \Bigg\Vert x_c(t)-{\sum_i} u_{i,c}(t)+\frac{\lambda_c(t)}{2}\Bigg\Vert^2_2

其显示解为

u^k,cn+1(ω)=x^c(ω)−∑i≠ku^i,c(ω)+λ^c(ω)21+2α(ω−ωk)2\hat{u}_{k,c}^{n+1}(\omega) = \frac{\hat{x}_c(\omega)-\sum_{i\neq k}\hat{u}_{i,c}(\omega)+\frac{\hat{\lambda}_c(\omega)}{2}}{1+2\alpha(\omega-\omega_k)^2}

  • 更新中心频率

固定{uk,c},λc\{u_{k,c}\},\lambda_c，可得如下子问题

min⁡ωk∑c∥∂t[uk,c,A(t)e−jωkt]∥22\min_{\omega_k}{\sum_c}\Bigg\Vert\partial_t\Big[u_{k,c,A}(t) e^{-j\omega_kt}\Big]\Bigg\Vert^2_2

利用Fourier变换转为频域上的优化问题为

min⁡ωk∑c∫0∞(ω−ωk)2∣u^k,c(ω)∣2dω\min_{\omega_k}{\sum_c}\int_{0}^{\infty}(\omega-\omega_k)^2\Big| \hat{u}_{k,c}(\omega)\Big|^2d\omega

其显示解为

ωkn+1=∑c∫0∞ω∣u^k,c(ω)∣2dω∑c∫0∞∣u^k,c(ω)∣2dω\omega_k^{n+1}=\frac{\sum_c\int^\infty_0 \omega|\hat{u}_{k,c}(\omega)|^2 d\omega}{\sum_c\int^\infty_0 |\hat{u}_{k,c}(\omega)|^2 d\omega}

# 算法流程

# 参考文献

  1. Dragomiretskiy K , Zosso D . Variational Mode Decomposition[J]. IEEE Transactions on Signal Processing, 2014, 62(3):531-544. ↩︎

  2. Rehman N U, Aftab H. Multivariate Variational Mode Decomposition[J]. IEEE Transactions on Signal Processing, 2019, 67(23): 6039-6052. ↩︎

变分非线性线性调频模态分解