上一篇谈了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])

 

        得：

http://s4/mw690/001XhlqSgy6PBSTE00z53&690
        图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

 

        看看运行结果：

http://s12/mw690/001XhlqSgy6PFnXMEEX9b&690
        图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')

http://s12/mw690/001XhlqSgy6PGOftCVZ6b&690
         图51-3 allmode_MB1387_01_100

 

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