信号处理——EMD、VMD的⼀点⼩思考芯撑
作者：桂。
时间：2017-03-06 20:57:22
链接：
前⾔
本⽂为的内容补充，主要内容为：
1）EMD原理介绍
2）代码分析
approval
3）⼀种权衡的⼩trick
4）问题补充
内容主要为⾃⼰的学习总结，并多有借鉴他⼈，最后⼀并给出链接。
⼀、EMD原理介绍
A-EMD的意义
很多⼈都知道EMD(Empirical Mode Decomposition)可以将信号分解不同频率特性，并且结合Hilbert求解包络以及瞬时频率。EMD、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基本原理
此处给⼀个原理图：
fuck you 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;
htuclo 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 purpos
wsub{1} = 2*pi*f_1;
wsub{2} = 2*pi*f_2;
wsub{3} = 2*pi*f_3;
% composite signal, including noi
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; % noi-tolerance (no strict fidelity enforcement) K = 4; % 4 modes
DC = 0; % no DC part impod
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 --- www.math.ucla.edu/~zosso
% Initial relea 2013-12-12 (c) 2013
%
% Input and Parameters:
% ---------------------
% signal - the time domain signal (1D) to be decompod
% alpha - the balancing parameter of the data-fidelity constraint
% tau - time-step of the dual ascent ( pick 0 for noi-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 decompod modes
% u_hat - spectra of the modes
% omega - estimated mode center-frequencies
%
% When using this code, plea do cite our paper:
% -----------------------------------------------
% K. Dragomiretskiy, D. Zosso, Variational Mode Decomposition, IEEE Trans. % on Signal Processing (in press)
% plea check here for update reference:
% dx.doi/10.1109/TSP.2013.2288675
%---------- Preparations
% Period and sampling frequency of input signal
save_T = length(signal);
continue的用法fs = 1/save_T;
% extend the signal by mirroring
T = save_T;
f_mirror(1:T/2) = signal(T/2:-1:1);
f_mirror(T/2+1:3*T/2) = signal;
f_mirror(3*T/2+1:2*T) = signal(T:-1:T/2+1);
f = f_mirror;
% Time Domain 0 to T (of mirrored signal)
T = length(f);
t = (1:T)/T;
% Spectral Domain discretization
freqs = t-0.5-1/T;
% Maximum number of iterations (if not converged yet, then it won't anyway) N = 500;
% For future generalizations: individual alpha for each mode
Alpha = alpha*ones(1,K);
% Construct and center f_hat
f_hat = fftshift((fft(f)));
f_hat_plus = f_hat;
f_hat_plus(1:T/2) = 0;
% matrix keeping track of every iterant // could be discarded for mem
u_hat_plus = zeros(N, length(freqs), K);
% Initialization of omega_k
没问题英文
omega_plus = zeros(N, K);
switch init
ca 1
for i = 1:K
omega_plus(1,i) = (0.5/K)*(i-1);
end
ca 2
omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));
otherwi
omega_plus(1,:) = 0;
end
% if DC mode impod, t its omega to 0
if DC
omega_plus(1,1) = 0;
end
% start with empty dual variables
lambda_hat = zeros(N, length(freqs));
% other inits
uDiff = tol+eps; % update step
n = 1; % loop counter
every dog has his day
sum_uk = 0; % accumulator
% ----------- Main loop for iterative updates
while ( uDiff > tol && n < N ) % not converged and below iterations limit
% update first mode accumulator
k = 1;
sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);
% 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(1,k)*(freqs - omega_plus(n,k)).^2);
% update first omega if not held at 0
if ~DC
omega_plus(n+1,k) = (freqs(T/2+1:T)*(abs(u_hat_plus(n+1, T/2+1:T, k)).^2)')/sum(abs(u_hat_plus(n+1,T/2+1:T,k)).^2); end
% update of any other mode
for k=2: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(1,k)*(freqs - omega_plus(n,k)).^2);
% center frequencies
omega_plus(n+1,k) = (freqs(T/2+1:T)*(abs(u_hat_plus(n+1, T/2+1:T, k)).^2)')/sum(abs(u_hat_plus(n+1,T/2+1:T,k)).^2); end
% Dual ascent
lambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus);
% loop counter
n = n+1;
% converged yet?
uDiff = eps;
for i=1:K
uDiff = uDiff + 1/T*(u_hat_plus(n,:,i)-u_hat_plus(n-1,:,i))*conj((u_hat_plus(n,:,i)-u_hat_plus(n-1,:,i)))';
end
uDiff = abs(uDiff);
end戴帽子的英文
%------ Postprocessing and cleanup
% discard empty space if converged early
N = min(N,n);
omega = omega_plus(1:N,:);
% Signal reconstruction
u_hat = zeros(T, K);
u_hat((T/2+1):T,:) = squeeze(u_hat_plus(N,(T/2+1):T,:));
u_hat((T/2+1):-1:2,:) = squeeze(conj(u_hat_plus(N,(T/2+1):T,:)));
u_hat(1,:) = conj(u_hat(end,:));
u = zeros(K,length(t));
