上一篇谈了RCADA的eemd函数存在的问题。这一篇谈谈我的eemd方法。我还是打算用G.Rilling编写的emd程序作eemd的基程序。但分析信号叠加各次白噪声再经emd分解后，所得imf分量个数不一样，怎么样把它们合成在一起呢？我的方法是分组：分量个数相同的imf放在同一组，总体计算只在组内进行。至于在相同的初始条件下为什么会分解出不同个数的imf分量来，这本身就是一个问题，应该另行研究，不应该回避。因为“真理只有一个”。

        先交代一个小问题。我在《29、黄锷院士的一个疏忽》中提到，用Rilling的emd函数计算所得到的jmf，偶尔会出现大于0的极小值或小于0的极大值情况，我称其为非严格的imf，否则，若所有的极大值都大于0，所有的极小值都小于0，则称其为严格的imf。由于我以后使用的我自编的瞬时频率估计函数，是针对严格的imf，因此我对Rilling的emd函数作了一点点修改，使其计算所得imf都是严格的imf。修改方法是将其反映“过零点与极值点数目相差不超过1”的相关代码改为“所有的极大值都大于0，所有的极小值都小于0”代码。为了以示区别，改后函数以emdb表示。

         将《29、黄锷院士的一个疏忽》中相关程序再运行一遍：

imf_MB2917_rd = emdb(MB2917_rd);

x1=imf_MB2917_rd(1,:);

[env,envmoy] = cenvelope(x1,1);%求上下包络线及其均值

env=env';

z=zeros(1,length(x1));

plot(x1,'r-')

hold on

plot(env)

plot(envmoy,'k--')

plot(z,'m-')

xlim([33900 34100])

        得：

        图51-1.

        与图29-1比较，可知已解决了“非严格imf”问题。 

        取一段数据MB1387，为本人2003、6、17～2007、4、3之间的脉搏信号。先用emdb计算其imf， 根据所得imf分量个数，设定加噪信号分解得到的的imf分组条件。

 程序如下：

x=MB1387;

nstd=0.1;

n=100;

[imf,o,v]=emdb(x);

STR0.imf0=imf;
STR0.ort0=o;

dv=length(v);%分解所得imf个数

lx=length(x);
sd=std(x);

i=1;

STR=[];
STR3=[];
STR2=[];
STR1=[];
STR_1=[];
STR_2=[];
STR_3=[];

m3=1;
m2=1;
m1=1;
m_1=1;
m_2=1;
m_3=1;

while i<=n
    r=randn(1,lx);
    ran=r*sd*nstd;
    xr=x+ran;

    [imf o v]=emdb(xr);    

    %分解所得imf个数等于原分析信号分解所得imf个数时
    if length(v)==dv

        STR(i).imf=imf;
        STR(i).ort=o;
        STR(i).ran=ran;
        i=i+1;

    %分解所得imf个数比原分析信号分解所得imf个数多出3个时
    elseif length(v)==dv+3
        STR3(m3).imf=imf;
        STR3(m3).ort=o;
        STR3(m3).ran=ran;
        m3=m3+1;

    %分解所得imf个数比原分析信号分解所得imf个数多出2个时   
    elseif length(v)==dv+2
        STR2(m2).imf=imf;
        STR2(m2).ort=o;
        STR2(m2).ran=ran;
        m2=m2+1;

    %分解所得imf个数比原分析信号分解所得imf个数多出1个时        
    elseif length(v)==dv+1
        STR1(m1).imf=imf;
        STR1(m1).ort=o;
        STR1(m1).ran=ran;
        m1=m1+1;
                
    %分解所得imf个数比原分析信号分解所得imf个数少出1个时
    elseif length(v)==dv-1
        STR_1(m_1).imf=imf;
        STR_1(m_1).ort=o;
        STR_1(m_1).ran=ran;
        m_1=m_1+1;
       
    %分解所得imf个数比原分析信号分解所得imf个数少出2个时
    elseif length(v)==dv-2
        STR_2(m_2).imf=imf;
        STR_2(m_2).ort=o;
        STR_2(m_2).ran=ran;
        m_2=m_2+1;
       
    %分解所得imf个数比原分析信号分解所得imf个数少出3个时
    elseif length(v)==dv-3
        STR_3(m_3).imf=imf;
        STR_3(m_3).ort=o;
        STR_3(m_3).ran=ran;
        m_3=m_3+1;
       
    end
end
end

        看看运行结果：

         图51-2

         加噪信号与原分析信号分解的imf分量个数不一样的，只有STR1与STR_1有数据，其余都是空集。经过大量运行，该数量差异，极少有2个以上的。说明程序中设置最多差3个，足够适应实际情况。

         将前面程序改成函数形式，调用格式为

[STR0,STR,STR3,STR2,STR1,STR_1,STR_2,STR_3]=imfstr(x,nstd,n);

        根据变量STR，计算一下总体imf：

allmode=zeros(size(STR(1).imf));
for j=1:n
    allmode=allmode+(STR(j).imf);
end
allmode=allmode/n;

%画图：

zs=zeros(1,lx);
W=1080;
H=720;
w=0.80;
h=0.70;
m=5;n=3;k=13;
fg=figure(3);
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);
            set(s(p),'position',[0.04+(j-1)/n,(1-0.95*i/m),w/n,h/m]);
            x=allimf(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
title('allmode_MB1387_01_100')

      图51-3 allmode_MB1387_01_100

        用肉眼直接可以看到，第9、第10分量仍旧存在“非imf极值点（大于0的极小值与小于0的极大值）”， 说明出现非严格的imf总体，主要是因为总体合成的方法，而不是因为emd函数。

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

首发时间：2015-02-03 14:12:37)）

（注：原创文章，转载请注明出处。欢迎对本博主题有兴趣的朋友加本人QQ554240962，注明“科网博友”）