信号处理——EMD、VMD的一点小思考

作者：桂。

时间：2017-03-06  20:57:22

链接：http://www.cnblogs.com/xingshansi/p/6511916.html 

---

前言

本文为Hilbert变换一篇的内容补充，主要内容为：

　　1）EMD原理介绍

　　2）代码分析

　　3）一种权衡的小trick

　　4）问题补充

内容主要为自己的学习总结，并多有借鉴他人，最后一并给出链接。

一、EMD原理介绍

　　A-EMD的意义

很多人都知道EMD(Empirical Mode Decomposition)可以将信号分解不同频率特性，并且结合Hilbert求解包络以及瞬时频率。EMD、Hilbert、瞬时频率三者有无内在联系？答案是：有。

按照Hilbert变换一篇的介绍，

$f(t) = \frac{ {d\Phi (t)}}{ {d(t)}}$

然而，这样求解瞬时频率在某些情况下有问题，可能出现$f(t)$为负的情况：我1秒手指动5下，频率是5Hz；反过来，频率为8Hz时，手指1秒动8下，可如果频率为-5Hz呢？负频率没有意义。

考虑信号

$x(t) = {x_1}(t) + {x_2}(t) = {A_1}{e^{j{\omega _1}t}} + {A_2}{e^{j{\omega _2}t}} = A(t){e^{j\varphi (t)}}$

为了简单起见，假设$A_1$和$A_2$恒定，且$\omega_1$和$\omega_2$是正的。信号$x(t)$的频谱应由两个在$\omega_1$和$\omega_2$的$\delta$函数组成，即

$X(\omega ) = {A_1}\delta (\omega  - {\omega _1}) + {A_2}\delta (\omega  - {\omega _2})$

因为假设$\omega_1$和$\omega_2$是正的，所以该信号解析。求得相位

$\Phi (t) = \frac{ { {A_1}\sin {\omega _1}t + {A_{\rm{2}}}\sin {\omega _{\rm{2}}}t}}{ { {A_1}\cos {\omega _1}t + {A_{\rm{2}}}\cos {\omega _{\rm{2}}}t}}$

分别取两组参数，对$t$求导，得到对应参数下的瞬时频率：

参数：

$\omega_1 = 10Hz$和$\omega_2 = 20Hz$.

• 组1：{$A_1 = 0.2, A_2 = 1$};
• 组2：{$A_1 = 1.2, A_2 = 1$}

对于组2，瞬时频率出现了负值。

可见：

对任意信号进行Hilbert变换，可能出现无法解释、缺乏实际意义的频率分量。Norden E. Hung等人对瞬时频率进行研究后发现，只有满足特定条件的信号，其瞬时频率才具有物理意义，并将此类信号成为：IMF/基本模式分量。　

　　B-EMD基本原理

此处给一个原理图：

　　C-基本模式分量（IMF）

EMD分解的IMF其瞬时频率具有实际物理意义，原因有两点：

• 限定1：
  • 在整个数据序列中，极值点的数量$N_e$（包括极大值、极小值点）与过零点的数量必须相等，或最多相差1个，即$(N_e-1) \le N_e \ge (N_e+1)$.
• 限定2：
  • 在任意时间点$t_i$上，信号局部极大值确定的上包络线$f_{max}(t)$和局部极小值确定的下包络线$f_{min}(t)$的均值为0.

限定1即要求信号具有类似传统平稳高斯过程的分布；限定2要求局部均值为0，同时用局部最大、最小值的包络作为近似，从而信号局部对称，避免了不对称带来的瞬时频率波动。

　　D-VMD

关于VMD(Variational Mode Decomposition)，具体原理可以参考其论文，这里我们只要记住一点：其分解的各个基本分量——即各解析信号的瞬时频率具有实际的物理意义。

 

二、代码分析

首先给出信号分别用VMD、EMD的分解结果:

给出对应的代码：

%--------------- Preparation
clear all;
close all;
clc;
% Time Domain 0 to T
T = 1000;
fs = 1/T;
t = (1:T)/T;
freqs = 2*pi*(t-0.5-1/T)/(fs);
% center frequencies of components
f_1 = 2;
f_2 = 24;
f_3 = 288;
% modes
v_1 = (cos(2*pi*f_1*t));
v_2 = 1/4*(cos(2*pi*f_2*t));
v_3 = 1/16*(cos(2*pi*f_3*t));
% for visualization purposes
wsub{1} = 2*pi*f_1;
wsub{2} = 2*pi*f_2;
wsub{3} = 2*pi*f_3;
% composite signal, including noise
f = v_1 + v_2 + v_3 + 0.1*randn(size(v_1));
% some sample parameters for VMD
alpha = 2000;        % moderate bandwidth constraint
tau = 0;            % noise-tolerance (no strict fidelity enforcement)
K = 4;              % 4 modes
DC = 0;             % no DC part imposed
init = 1;           % initialize omegas uniformly
tol = 1e-7;

%--------------- Run actual VMD code
[u, u_hat, omega] = VMD(f, alpha, tau, K, DC, init, tol);
subplot(size(u,1)+1,2,1);
plot(t,f,'k');grid on;
title('VMD分解');
subplot(size(u,1)+1,2,2);
plot(freqs,abs(fft(f)),'k');grid on;
title('对应频谱');
for i = 2:size(u,1)+1
    subplot(size(u,1)+1,2,i*2-1);
    plot(t,u(i-1,:),'k');grid on;
    subplot(size(u,1)+1,2,i*2);
    plot(freqs,abs(fft(u(i-1,:))),'k');grid on;
end

%---------------run EMD code
imf = emd(f);
figure;
subplot(size(imf,1)+1,2,1);
plot(t,f,'k');grid on;
title('EMD分解');
subplot(size(imf,1)+1,2,2);
plot(freqs,abs(fft(f)),'k');grid on;
title('对应频谱');
for i = 2:size(imf,1)+1
    subplot(size(imf,1)+1,2,i*2-1);
    plot(t,imf(i-1,:),'k');grid on;
    subplot(size(imf,1)+1,2,i*2);
    plot(freqs,abs(fft(imf(i-1,:))),'k');grid on;
end

　　附上两个子程序的code.

VMD:

function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol)
% Variational Mode Decomposition
% Authors: Konstantin Dragomiretskiy and Dominique Zosso
% zosso@math.ucla.edu --- http://www.math.ucla.edu/~zosso
% Initial release 2013-12-12 (c) 2013
%
% Input and Parameters:
% ---------------------
% signal  - the time domain signal (1D) to be decomposed
% alpha   - the balancing parameter of the data-fidelity constraint
% tau     - time-step of the dual ascent ( pick 0 for noise-slack )
% K       - the number of modes to be recovered
% DC      - true if the first mode is put and kept at DC (0-freq)
% 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
%
% Output:
% -------
% u       - the collection of decomposed modes
% u_hat   - spectra of the modes
% omega   - estimated mode center-frequencies
%
% When using this code, please d