中医现代化研究分享 http://blog.sciencenet.cn/u/baishp 当你们都还在想象时,我已经在路上了。这注定是一场一个人的战争吗?

博文

51、根据imf分量个数分组_我的eemd方法(一)

已有 5271 次阅读 2018-3-27 15:05 |个人分类:斤斤计较|系统分类:科研笔记| 总体经验模式分解, eemd, 本征模函数, imf

 

       上一篇谈了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据imf分量ä¸a数分组_我的eemdæ–1æ3•(一)

        图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、æ 1据imf分量ä¸a数分组_我的eemdæ–1æ3•(一)

         图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、æ 1据imf分量ä¸a数分组_我的eemdæ–1æ3•(一)

      图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,注明“科网博友”)




https://blog.sciencenet.cn/blog-825323-1105941.html

上一篇:50、RCADA程序包中eemd函数存在的问题
下一篇:52、宇宙膨胀之我见
收藏 IP: 219.133.173.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-11-9 07:26

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部