本博有两年多没有更新了。现在继续来看看EMD（经验模态分解）。前年（2012）较早的时候，为了了解EEMD（总体经验模态分解），我本来进去过EMD的“官方网站”RCADA（国立中央大学数据分析方法研究中心）网站的，当时下载了6、7段程序代码，以为全都是关于EEMD的，就那么多东西了，所以也就没再多关注而继续我自己的“瞬时频率估计方法研究”。今年又去RCADA网站看了，下了一个程序包，有40多段程序代码，知道他们对“瞬时频率估计方法”作了详尽的研究，各种方法总共有6种之多。可以说，他们已经穷尽了他们所能想到的一切方法了。进一步的了解，知道这6种方法没有一种跟我的方法相同的，而且我的方法比他们的所有方法都要好。关于其运行结果的对比分析，以后我会详细地贴出。 

        先来看看eemd的情况。RCADA程序包中，函数allmode=eemd(Y,Nstd,NE)，Y是输入的一维分析数据，Nstd是所添加噪声的标准差与Y的标准差之比，NE是为获取总体imf（本征模函数）而进行的加噪imf的叠加次数。allmode是一个L*N矩阵，L是Y的长度，N＝fix(log2(L))+1，其第一列是原数据Y，其余列是总体imf（最后一列是残余项）。

        这个eemd函数是有很大问题的。最大的问题，就是由它计算得到的总体allmode并不是严格意义的imf。先来运行一段数据看看：

Y=GY1387;%（本人6/17/2003～4/3/2007之间的收缩压数据）
Nstd=0.1;
NE=100;
allmode=eemd(Y,Nstd,NE);

        运行得到一个16644×15的矩阵。作图：

lx=length(Y);
zs=zeros(1,lx);
W=1080;
H=720;
w=0.80;
h=0.70;
m=5;n=3;k=15;
fg=figure(1);
set(fg,'position',[50,50,W,H]);

for i=1:m
    for j=1:n
        p=j+(i-1)*n;
        if p<=k
            s(p)=subplot(m,n,p);
            x=allmode(:,p);
            plot(x)
            hold on
            plot(zs,'k--')
            mx=max(x);
            mn=min(x);
            xn=mx-mn;

            axis([fix(-0.1*lx) fix(1.1*lx) mn-0.1*xn mx+0.1*xn])
        end
    end
end

         图50-1

        很容易看出，第9个分量（（4，1）子图）、第11个分量（（4，3）子图）都含有非imf极值点（大于0的极小值）。

        再看看高频imf分量情况：

fg=figure(2);
set(fg,'position',[50,50,W,H]);
set(gca,'position',[0.035,0.06,0.95,0.88]);
x=allmode(:,3);
plot(x)
hold on
plot(zs,'k--')
xlim([5200,5400])
title('allmode(:,3) xlim([5200,5400])')

        图50-2

        这是第2个分量局部放大图，可以看出同样有很多的非imf极值点（大于0的极小值与小于0的极大值）。

        我们知道，成为imf的两个条件（1）上下包络线对称（2）极值点与过零点数相差不超过1，它的思维出发点，是将原信号分解成一个个绕原点作单向环周运动而不含任何骑波的信号。因此，经过噪声叠加、分解合成处理后得到的总体，也应该仍然是一组imf，即也应该符合上面两个条件才好，否则违背了初衷，整个EMD研究也就失去意义了。遗憾的是，eemd函数求出的总体，不是严格意义的imf。

        特别地，当Nstd=0,NE=1时，eemd就退化成一个纯粹的emd函数。比起我以前使用的G.Rilling编写的emd函数来，其计算过程要简单得多。概括说来，就是（1）每次分解的imf个数由fix(log2(L))-1确定，这里L是数据Y的长度，（2）每个imf的筛选迭代次数规定为10次，而不是象以前的emd程序那样由自适应条件决定。这样做的目的是为了确保每次分解得到的imf个数是一样的，因为只有各次imf的个数一样时，才可能最后叠加合成成总体imf（allmode）。

         当Nstd≠0,NE较小（如取“1”、“2”）时，可以看到输出结果，存在明显问题。

取：

Y=GY1387;

Nstd=0.1;
NE=1;
allmode=eemd(Y,Nstd,NE);

 运行并作图，得：

        图50-3

        可以看到，第12个分量（（5，1）子图）、第13个分量（（5，2）子图）全部在0值线以上。容易查知，这两个分量的首尾端点的函数值都等于0。它们都不是imf。

        另一种反常情况，就是在低频端出现幅值极小的高频分量，见图：

      图50-4

        图50-5

        这种非正常的分量，都出现在低频端，其具体位置与个数不定。

       当出现全非负或全非正分量时，使用程序包中的瞬时频率估计函数ifndq，无法对该分量作频率估计，代码不能运行。如图：

     图50-6

       图50-7

        说明eemd函数问题多多，不能当真的。

（本文首发于：http://blog.sina.com.cn/s/blog_6ad0d3de0102v6zi.html

首发时间：2014-12-21 14:26:38)