EMD——学习笔记

简介

  语音信号是非平稳信号，对非平稳信号分析比较直观的方法是使用具有局限性的基本量和基本函数，例如瞬时频率。1998年，Norden E. Huang 等人对瞬时频率做了深入的研究，提出了本征模式函数(IMF)和将任意信号分解为IMF组成的新方法——经验模式分解(EMD)。

基本概念

  在了解EMD算法之前，需要先了解两个基本概念：瞬时频率；基本模式分量。当然如果了解这两个概念可以直接跳过。

瞬时频率

对任意的时间序列$x(t)$,可以得到Hilbert变换
$$y(t)=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\frac{x(\tau )}{t-\tau }d\tau \text{\hspace{1em}(1)}$$
构造解析函数（满足C-R条件的可微的函数）
$$z(t)=x(t)+iy(t)=a(t)e^{i\Phi (t)}\text{\hspace{1em}(2)}$$
其中幅值函数
$$a(t)=\sqrt{x^2(t)+y^2(t)}\text{\hspace{1em}(3)}$$
相位函数
$$\Phi (t)=arctan\frac{y(t)}{x(t)}tag4$$
显然相位函数对时间的导数就是瞬时频率$$w(t)=\frac{d\Phi (t)}{dt}\text{\hspace{1em}(5)}$$
即
$$f(t)=\frac{1}{2\pi }\frac{d\Phi (t)}{dt}\text{\hspace{1em}(6)}$$
  显然瞬时频率会出现负值的情况，可以得出，对任意信号做Hilbert变换会出现无法解释，缺乏实际物理意义的频率成分。后来经过Norden E. Huang 等人的深入研究，发现只有满足一定条件的信号才能求求得具有物理意义的频率，并且将此类信号成为基本模式分量（IMF）。

基本模式分量

上文已经说了，对IMF进行Hilbert变换出来的时频谱才有具体的物理意义。基本模式分量$f(t)$需要满足下面两个条件：

• 在整个数据序列中，极值点个数($N_e$)与过零点的个数($N_s$)必须相等或者相差不多于一个
  $$|N_e-N_s|\le 1\text{\hspace{1em}(7)}$$
• 在任意时间点$t_i$上，信号局部极大值上的包络线$f_{max}(t)$和局部极小值确定的包络线$f_{min}(t)$的均值为零，即
  $$[f_{max}(t_i)+f_{min}(t_i)]/2=0,t_i\in [t_a,t_b]\text{\hspace{1em}(8)}$$

基本原理

  对于满足基本模式分量的两个限定条件的信号，可以通过Hilbert变换求出其瞬时频率。但是大多数信号都不是基本模式分量。所以我们就要把信号变成基本模式分量，Norden E. Huang 等人创造性的提出了假设：任何信号都是又一些不同的基本模式分量组成。并且进一步提出可以通过EMD将信号的基本模式提取出来。
  信号模式分解的目的：得到使得瞬时频率有意义的时间序列——基本模式分量。
$f_{max}(t)$ 和 $f_{min}(t)$分别是$x(t)$的上下包络

• $$m(t)=[f_{max}(t)+f_{min}(t)]/2\text{\hspace{1em}(9)}$$
• $$h(t)=x(t)-m(t)\text{\hspace{1em}(10)}$$
  当$h(t)$不是基本分量时，$x(t)=h(t)$,重复以上步骤，直到$h(t)$为基本分量时继续以下步骤
• $$r(t)=x(t)-h(t)\text{\hspace{1em}(11)}$$
  当$r(t)$不是基本模式分量时，$x(t)=r(t)$ 重复以上步骤，直到满足停止准则为止。
  停止准则：
• 最后一个基本模式分量$h(t)$或者剩余分量$r(t)$比预期值小时停止
• 剩余分量$r(t)$变成单调函数，不能再从中筛选出基本模式分量
  当然要达到上面这两种情况之一，可能要很多次的重复筛选，可能会导致基本模式分量失去实际的物理意义，因此定义一个准则来实现
  $$S_d=\sum _{t=0}^T\frac{|h_{k-1}(t)-h_k(t){|}^2}{h_k^2(t)}\text{\hspace{1em}(12)}$$
  $T$表示信号的时间长度；$h_{k-1}(t)$ 和 $h_k(t)$表示在筛选过程中两个连续的处理结果的时间序列。$S_d$的取值通常为0.2-0.3。

MATLAB 代码

clear all; clc; close all;
fs=5000;                       % 采样频率
N=500;                         % 样点数
n=1:N;
t1=(n-1)/fs;                   % 设置时间
x1=sin(2*pi*50*t1);            % 产生笫1个正弦信号
x2=(1/3)*sin(2*pi*150*t1);     % 产生笫2个正弦信号
z=x1+x2;                       % 把两个信号叠加
imp=emd(z);                    % 对叠加信号进行EMD分解
[n,m]=size(imp);               % 求取EMD分解成几个分量
% 作图
subplot(m+1,1,1);              % 画叠加信号
plot(t1,z,'k');title('原始信号'); ylabel('幅值')
subplot (m+1,1,2);                   % 画第1个正弦信号
line(t1,x2,'color',[.6 .6 .6],'linewidth',5); hold on
subplot (m+1,1,3);                   % 画第2个正弦信号
line(t1,x1,'color',[.6 .6 .6],'linewidth',5); hold on
for i=1:m
   subplot(m+1,1,i+1);        % 画EMD分解后的信号
   plot(t1,imp(:,i),'k','linewidth',1.5); ylabel('幅值')
   title(['imf' num2str(i)]);
end
xlabel('时间/s');

如图所示，粗线是原始信号，黑线是分出来的信号，很好的重合了，两端没有重合是因为端点效应。
参考书籍： MATLAB语音信号分析与合成