基于G.Rilling的极值镜像延拓边界处理方法的改进

529899778

      之前看到《希尔伯特-黄变换方法的改进》中提出“平行延拓”的方法来对端点效应问题进行改善处理：
“利用端点处附近的两个相邻极值点（一极大值，以极小值）处斜率相等这一特性，认为在边端处定义出两个极值点，分别连接相邻的极大值与极小值，对包络线进行拟合”。

                  
       上图为第一个极值是极小值与最后一个极值是极大值的情况，其他三种情况类似。
       针对图2中的情况，G.Rilling方法（1）若第一个采样值小于第一个极大值，则以第一个极小值点所在时间点镜像延拓；      （2）若第一个采样值大于第一个极小值，则以第一个采样点所在时间点镜像延拓。本人欲结合两种方法，进行镜像延拓，例如情况（2）中，用平行延拓方法预测得到一个极小值点，镜像延拓后样条拟合再求均值。
        通过仿真信号发现并没有很好的改善端点效应问题，想请教论坛的朋友们，是本人理解有错还是程序有问题。

程序附上 emd_Rilling_paracon.txt (11.19 KB, 下载次数: 42)

2012-1-12 17:09 上传

点击文件名下载附件

（后缀改为.m即可运行）

529899778

对程序稍微进行了修改，还是直接将相对G.Rilling的程序修改的部分直接贴出来。

% boundary conditions for interpolations :
    %%%插值的边界条件   
     if indmax(1) < indmin(1)%第一个极值点是极大值
      if m(1) > m(indmin(1))%如果第一个采样点大于第一极小值点，以第一个极大值为对称中心
        lmax = fliplr(indmax(2:min(end,NBSYM+1)));
        lmin = fliplr(indmin(1:min(end,NBSYM)));
        lsym = indmax(1);
      else%如果第一个采样点小于第一个极小值点，则认为第一个采样点为极小值点，以第一个采样点的预测点为对称中心
        lmax = [fliplr(indmax(1:min(end,NBSYM-1))),1];
        lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];
        lsym = 1;
      end
     else%第一个极值点是极小值
      if m(1) < m(indmax(1))%以第一个极小值为对称中心
        lmax = fliplr(indmax(1:min(end,NBSYM)));
        lmin = fliplr(indmin(2:min(end,NBSYM+1)));
        lsym = indmin(1);
      else%以第一个采样点为对称中心
        lmax = [fliplr(indmax(1:min(end,NBSYM-1))),1];
        lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];
        lsym = 1;
      end
    end
    %%%末尾序列与开头序列类似
    if indmax(end) < indmin(end)%最后一个极值是极小值
      if m(end) < m(indmax(end))%最后一个采样点小于最后一个极大值点，以最后一个极小值点为对称中心
        rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
        rmin = fliplr(indmin(max(end-NBSYM,1):end-1));
        rsym = indmin(end);
      else%最后一个采样点大于最后一个极值点，认为该点为极大值点，以该点的预测点为对称中心
        rmax = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];
        rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];
        rsym = lx;
      end
    else
      if m(end) > m(indmin(end))
        rmax = fliplr(indmax(max(end-NBSYM,1):end-1));
        rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
        rsym = indmax(end);
      else
        rmax = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];
        rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];
        rsym = lx;
      end
    end
    %%%将序列根据中心对称镜像延拓到两边
    tlmin = 2*t(lsym)-t(lmin);
    tlmax = 2*t(lsym)-t(lmax);
    trmin = 2*t(rsym)-t(rmin);
    trmax = 2*t(rsym)-t(rmax);

    % in case symmetrized parts do not extend enough
    %%%如果对称的部分没有足够的极值点
    if tlmin(1) > t(1) || tlmax(1) > t(1)%对折后的序列没有超出原序列的范围
      if lsym == indmax(1)
        lmax = fliplr(indmax(1:min(end,NBSYM)));
      else
        lmin = fliplr(indmin(1:min(end,NBSYM)));
      end
      if lsym == 1%以第一个采样点为对称中心镜像，但对折后的序列不超出原序列的范围的这种情况不应该出现，若出现，则程序中止
        error('bug')
      end
      lsym = 1;%直接以第一采样点为对称中心取镜像
      tlmin = 2*t(lsym)-t(lmin);
      tlmax = 2*t(lsym)-t(lmax);
    end   
    %%%序列末尾情况与序列开头类似
    if trmin(end) < t(lx) || trmax(end) < t(lx)
      if rsym == indmax(end)
        rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
      else
        rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
      end
      if rsym == lx
        error('bug')
      end
      rsym = lx;
      trmin = 2*t(rsym)-t(rmin);
      trmax = 2*t(rsym)-t(rmax);
    end
    %%%延拓点上的取值      
    mlmax =m(lmax);
    mlmin =m(lmin);
    mrmax =m(rmax);
    mrmin =m(rmin);

    tmin = [tlmin t(indmin) trmin];
    tmax = [tlmax t(indmax) trmax];
    mmin = [mlmin m(indmin) mrmin];
    mmax = [mlmax m(indmax) mrmax];

    if indmax(1)<indmin(1)%第一个极值是极大值
        if m(1)<m(indmin(1))%第一个采样值小于第一个极小值，需要预测一个极小值
            slope=(m(indmin(1))-m(1))/(indmin(1)-1);
            b=m(indmax(1))-slope*indmax(1);
            y1_lpred=slope+b;
            mmax(length(mlmax))=y1_lpred;
        end
    else第一个极值是极小值
        if m(1)>m(indmax(1))%第一个采样值大于第一个极大值，需要预测一个极大值
            slope=(m(indmax(1))-m(1))/(indmax(1)-1);
            b=m(indmin(1))-slope*indmin(1);
            y2_lpred=slope+b;
            mmin(length(mlmin))=y2_lpred;
        end
    end
    %%%序列末尾情况与开头情况类似
    if indmax(end)<indmin(end)
        if m(end)>m(indmax(end))
            slope=(m(end)-m(indmax(end)))/(lx-indmax(end));
            b=x(indmin(end))-slope*indmin(end);
            y1_rpred=slope*lx+b;
            mmin(length(mlmin)+length(m(indmin))+1)=y1_rpred;
        end
    else
        if m(end)<m(indmin(end))
            slope=(m(end)-x(indmin(end)))/(lx-indmin(end));
            b=m(indmax(end))-slope*indmax(end);
            y2_rpred=slope*lx+b;
            mmax(length(mlmax)+length(m(indmax))+1)=y2_rpred;
        end
    end

完整程序 emd_Rilling_paracon.txt (10.84 KB, 下载次数: 20)

2012-1-13 09:50 上传

点击文件名下载附件

playfish

function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)
    lx=length(x);

    indmax=[1 indmax lx];
    indmin=[1 indmin lx];

    zlmax=x(1);
   
    if length(indmax)>=4
        slope1=(x(indmax(3))-x(indmax(2)))/(indmax(3)-indmax(2));
        tmp1=x(indmax(2))-slope1*(indmax(2)-indmax(1));
        if tmp1>x(1)
            zlmax=tmp1;
        end
    end


    zrmax=x(end);
   
    if length(indmax)>=4
        slope2=(x(indmax(end-1))-x(indmax(end-2)))/(indmax(end-1)-indmax(end-2));
        tmp2=x(indmax(end-1))+slope2*(indmax(end)-indmax(end-1));
        if tmp2>x(end)
            zrmax=tmp2;
        end
    end


    zlmin=x(1);
   
    if length(indmin)>=4
        slope3=(x(indmin(3))-x(indmin(2)))/(indmin(3)-indmin(2));
        tmp3=x(indmin(2))-slope3*(indmin(2)-indmin(1));
        if tmp3<x(1)
            zlmin=tmp3;
        end
    end


    zrmin=x(end);
   
    if length(indmin)>=4
        slope4=(x(indmin(end-1))-x(indmin(end-2)))/(indmin(end-1)-indmin(end-2));
        tmp4=x(indmin(end-1))+slope4*(indmin(end)-indmin(end-1));
        if tmp4<x(end)
            zrmin=tmp4;      
        end
    end

        tmin = t(indmin);
        tmax = t(indmax);
        zmin = [zlmin z(indmin(2:end-1)) zrmin];
        zmax = [zlmax z(indmax(2:end-1)) zrmax];
end

playfish

这是我根据Huang在《Ensemble Empirical Mode Decomposition: a Noise-Assisted Data Analysis Method》中给出的算法写的程序，感觉不错。发出来给大家用。

529899778

playfish 发表于 2012-2-4 12:45

                               
登录/注册后可看大图

这是我根据Huang在《Ensemble Empirical Mode Decomposition: a Noise-Assisted Data Analysis Method》中给 ...

非常感谢 先把您的程序尽量读懂 再请教您

529899778

playfish 发表于 2012-2-4 12:45

                               
登录/注册后可看大图

这是我根据Huang在《Ensemble Empirical Mode Decomposition: a Noise-Assisted Data Analysis Method》中给 ...

       您提到的文章之前下载过，Huang的EEMD的程序我也有，程序基本看懂了，就是用邻近端点的两个极值点求取一个值来代替端点，想必《希尔伯特-黄变换方法的改进》的作者也是看过这篇文章的，由于时间较紧张，文章没有细读，想请教您这样处理端点的理论依据是什么，您用该该方法所得记过的效果如何？

529899778

回复 4 # playfish 的帖子

我根据你的程序对G.Rilling的程序进行了修改，然后对下边的信号进行imf分解：
fs=1000;
f1=50;
f2=100;
t=1/fs:1/fs:2;
x1=2*sin(2*pi*f1*t+0.2*pi);
x2=5*cos(2*pi*f2*t+0.5*pi);
x=[x1 x2];
clear fs f1 f2 t x1 x2

分解得到5个imf分量，第一个imf分量基本与原信号一致，后4个imf分量都是低频的，但最后一个imf分量明显不是低频走势，因为不是单调的，越来越困惑了...........

529899778

529899778 发表于 2012-2-9 16:59

                               
登录/注册后可看大图

您提到的文章之前下载过，Huang的EEMD的程序我也有，程序基本看懂了，就是用邻近端点的两个极值点 ...

更正下，用的是与端点邻近的四个极值点。

playfish

回复 8 # 529899778 的帖子

g rilling的程序现在看已经不能用了。最好自己写一个。目前最新的进展指出，停止条件就是迭代10次就停止最好。然后NHT分离AM和FM。再对FM做DQ求出瞬时频率。

playfish

回复 6 # 529899778 的帖子

比rilling的镜像法更快收敛。边界抖动差不多。

529899778

回复 9 # playfish 的帖子

想请问playfish哪篇文献中有对“停止条件就是迭代10次就停止最好”这个问题的说明，确实，对仿真信号进行分解程序运行时间还是较短的，如果对实测信号分解时间就比较长了，这个需要改进下。

529899778

回复 10 # playfish 的帖子

对，我现在最头疼的问题就是端点效应，3个月的学习和编程，用6种方法（晚点我会发帖把我做的写出来，希望playfish能帮帮我）都无法得到一个好点的效果，以至于毕设无法往下做，真是悲了个剧啊......

playfish

回复 11 # 529899778 的帖子

注意看Huang的论文列表，2010年的两篇就是讨论这个的
http://rcada.ncu.edu.tw/research1_clip_reference.htm

529899778

回复 13 # playfish 的帖子

谢谢，确实，还是应该看看老Huang的文章，毕竟是他搞出来的算法，要权威些...

lxppopo

playfish 发表于 2012-2-4 12:44

                               
登录/注册后可看大图

function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)
    lx=length(x);
...

这个貌似出不来结果啊，亲