关于EMD分解的边界效应处理

zhangcheng1108

       请教论坛上的各位高手：我打算采用AR模型处理EMD的边界效应，现在用AR模型进行端点延拓的部分已经完成，现在遇见的问题主要是两点：

1.是否应该先对原来的信号进行端点预测延拓以后在利用论坛上的emd.m进行分解，还是要对emd.m进行修改，把端点延拓的部分加到emd.m中，再对原信号分解？

2.在对信号进行端点预测的时候，应该延拓多少点为好？个人发现不同的文献记载不同的结论。有的是说要延拓一个采样周期，有的是说要在两侧各延拓一个极大值点和极小值点，还有的是说，要在信号两端各延拓两个极大值点和极小值点。不知哪种说法正确。

谢谢论坛上的各位高手予以解答

xray

问题1：需要对emd.m进行修改，把端点延拓的部分加到emd.m中。
理由：在EMD中，每次迭代过程都需要对端点进行延拓（这个是由emd.m自动调用的），如果仅仅是修改原始信号，那么在迭代过程中的端点延拓就无法使用AR模型进行处理了。

问题2：这个没有定论，可以根据试验效果来选择延拓点的个数。

顺便说一句：请不要重复发帖！

[ 本帖最后由 xray 于 2008-6-25 16:04 编辑 ]

zhangcheng1108

谢谢xray及时回复，我不是有意重复发贴的，之所以发了两个是因为我觉得这是两个不同的问题。我下次一定注意。避免类似现象的发生。
还想请教一下你，若对emd.m进行修改，应该怎样进行？我一直在看emd,m的源程序，确实有的地方不是很明白

xray

主要修改 emd.m 中下面这个函数的代码就可以了
function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)

zhangcheng1108

谢谢xray师兄的及时回复，我回去一定努力研究

cleverblue

又没有研究波形匹配法解决边界拓延问题的?

cleverblue

if (BOUNDMETHOD == 0)  % 0:边界延拓镜像法; 1:边界延自适应波形匹配法
     [tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,m,m,NBSYM);        % 边界延拓
else                   % 1:边界延自适应波形匹配法
     [tmin,tmax,mmin,mmax] = boundary_matchwave(indmin,indmax,t,m);
end;
  if (ENVMETHOD == 0 )                       % 0:上下包络法求均值; 1:用中值定理法求均值  
      envmin = interp1(tmin,mmin,t,INTERP);  % 插值函数为interp1,其中x，y为插值点，yi为在被插值点xi处的插值结果；x,y为向量， 表示采用的插值方法
      envmax = interp1(tmax,mmax,t,INTERP);
      envmoy = (envmin+envmax)/2;            % 计算包络均值
      if nargout > 3
         amp = mean(abs(envmax-envmin),1)/2; % 计算包络幅度
      end
  else                                       % 1:用中值定理法求均值
      [mvalue,indt] = compose(tmin,tmax,mmin,mmax);
      [env,tindex] = mvtd(m,mvalue,indt,0);    % 用中值定理法求均值
      envmoy = interp1(tindex,env,t,INTERP);   % 插值函数x = 0:10; y = sin(x); xx = 0:.25:10;yy = spline(x,y,xx); plot(x,y,'o',xx,yy)  CSAPI(0
      if nargout > 3
         amp = mean(abs(envmoy),1)/2; % 计算包络幅度
      end
  end;

mingtiangenghao

楼上的，请把boundary_matchwave子程序发上来，学习学习。

cleverblue

function [tmin,tmax,zmin,zmax] = boundary_matchwave(indmin,indmax,t,x)
%
%  一种自适应波形匹配端点拓延法
% 输入参数说明:  indmin      极小值地址序列
%                indmax      极大值地址序列
%                t           时间序列(时序)1:length(x)
%                x           时间序列(值)
%
% 输出参数说明:  tmin    极小值地址序列
%                tmax    极大值地址序列
%                zmin    极小值
%                zmax    极大值
%

lx = length(x);

% 判断极值点个数
if (length(indmin) + length(indmax) < 3)
  error('not enough extrema')
end

% 左端处理
leftmax = indmax(1);
leftmin = indmin(1);
if leftmin > leftmax
  left = leftmin;      % 左边处理长度
  leftindex = indmin;
  leftminind = indmax;
  leftdisc = leftmax;  % 末端长度
else
  left  = leftmax;
  leftindex = indmax;
  leftminind = indmin;
  leftdisc = leftmin;
end;
leftvalue = x(1:left);
indexlen = length(leftindex);
mintol = 99999;
for i = 2:indexlen
    len = leftindex(i);
    leftvalue1 = x(len-left+1:len);
    [tol,err,extrx1] = SelfAdapMatchWave(leftvalue,leftvalue1);
    if tol < mintol
       mintol = tol;
       avg=(sum(abs(extrx1))+sum(abs(rightvalue)))/(length(extrx1)*2-2);
       error = rightvalue - extrx1;
       avgerr=sum(abs(error))/(length(error)-1);
       minextrx1 = extrx1;
       newerr = err;
       minindex = len;
       indmini = i;
    end
end;
%minextrx1
if (mintol < avg/(10*2) || mintol < avgerr*2)  % 没有找到匹配子波,左端处理,取原始信号最左端的两个相邻极大值点的均值作为左端点的极大值，取信号最左端的两个相邻极小值点的均值作为左端点的极小值
    lindmin = indmin(1)-indmin(2);
    lvmin = (x(indmin(1))+x(indmin(2)))/2;
    lindmax = indmax(1)-indmax(2);
    lvmax = (x(indmax(1))+x(indmax(2)))/2;
else             % 找到匹配子波,左端处理
  leftindex = leftminind(indmini);
  leftx = x(leftindex-left-1:leftindex)+newerr;
  leftx = [leftx,x(1:2)];
  [lhindmin,lhindmax,lhindzer] = extr(leftx);
  lindmin = lhindmin(end)-(length(leftx)-2);
  lindmax = lhindmax(end)-(length(leftx)-2);
  lvmin = leftx(lhindmin(end));
  lvmax = leftx(lhindmax(end));
end;
% 右端处理
rightmax = indmax(end);
rightmin = indmin(end);
if rightmin < rightmax
  right = lx - rightmin;   % 右边处理长度
  rightindex = indmin;
  rightmaxind = indmax;
  rightdisc = lx - rightmax;  % 末端长度
else
  right  = lx - rightmax;
  rightindex = indmax;
  rightmaxind = indmin;
  rightdisc = lx - rightmin;
end;
rightvalue = x(end-right+1:end);
indexlen = length(rightindex)-1;
mintol = 99999;
for i = 1:indexlen
    len = rightindex(i);
    rightvalue1 = x(len:len+right-1);   
    [tol,err,extrx1] = SelfAdapMatchWave(rightvalue,rightvalue1);
    if tol < mintol
       mintol = tol;
       avg=(sum(abs(extrx1))+sum(abs(rightvalue)))/(length(extrx1)*2-2);
       error = rightvalue - extrx1;
       avgerr=sum(abs(error))/(length(error)-1);
       minextrx1 = extrx1;
       newerr = err;
       minindex = len;
       indmini = i;
    end
end;
%minextrx1
if (mintol < avg/(10*2) || mintol < avgerr*2) % 没有找到匹配子波,右端处理,取原始信号最右端的两个相邻极大值点的均值作为右端点的极大值，取信号最右端的两个相邻极小值点的均值作为右端点的极小值
    rindmin = lx + (indmin(end)-indmin(end-1));
    rvmin = (x(indmin(end-1))+x(indmin(end)))/2;
    rindmax = lx + (indmax(end)-indmax(end-1));
    rvmax = (x(indmax(end-1))+x(indmax(end)))/2;
else             % 找到匹配子波,右端处理
  rightindex = minindex;  % rightmaxind(indmini)
  rightx = x(rightindex+right:rightindex+2*right+5)+newerr;
  rightx =[x(rightmaxind(end)+1:end),rightx];
  [rhindmin,rhindmax,rhindzer] = extr(rightx);
  rvmin = rightx(rhindmin(1));
  rvmax = rightx(rhindmax(1));
  rindmin = lx + (rhindmin(1)-length(x(rightmaxind(end):end)));
  rindmax = lx + (rhindmax(1)-length(x(rightmaxind(end):end)));
end;
% 取拓延的点,拓延到两边(极大值与极小值)
tlmin = lindmin;
tlmax = lindmax;
trmin = rindmin;
trmax = rindmax;


% 延拓点上的取值(极大值与极小值)      
zlmax = lvmax;
zlmin = lvmin;
zrmax = rvmax;
zrmin = rvmin;
% 完成延拓
tmin = [tlmin t(indmin) trmin];
tmax = [tlmax t(indmax) trmax];
zmin = [zlmin x(indmin) zrmin];
zmax = [zlmax x(indmax) zrmax];

mingtiangenghao

先学习:lol
:lol

zhangcheng1108

最近在研究EMD方法的边界效应处理的时候发现：造成HHT边界效应主要有两个原因：
（1）在应用三次样条曲线求包络谱的时候，这时的端点效应应用emd.m中的镜像延拓的方法已经解决了，效果还不错，有兴趣的同门可试试去掉emd.m中的边界延拓部分的程序，再看看。分解出的IMF分量在端点处发散的很厉害
（2）在对IMF分量作hilbert变换的时候也会出现端点效应，具体造成的原因我还不是很清楚，若有师兄了解，请赐教。这部分的端点效应的克服可以应用基于AR模型以及其他各种方法对各个IMF分量向双方向延拓，将边界效应转移到实际信号外，再截断。

sssbbbhhh1

hilbert变换怪现象！
http://forum.vibunion.com/forum/viewthread.php?tid=41401
这里有说明 这是由于hilbert变换本身是通过DFT计算的，对于非周期采样，有泄露

graceward

cleverblue朋友所给的波形匹配算法缺少一个函数： [tol,err,extrx1] = SelfAdapMatchWave(rightvalue,rightvalue1);
能把这个函数发一下吗？谢谢   另外你的程序是按照什么算法编的啊？有文章吗？学习一下

cleverblue

function [tol,err,extrx1] = SelfAdapMatchWave(x0,x1)
%
%  一种自适应波形匹配方法
% 输入参数说明:  x0      需要匹配的波形数据,长度为l
%                x1      被匹配的波形数据,长度为l
%
% 输出参数说明:  tol    x0,x1的波形匹配度
%                extrx1 x1的极值(extremum)序列,平移后的
%

len0 = length(x0);
len1 = length(x1);
if (len0 ~= len1)
    error('x0,x1数据的长度不一致');
end;
% 数据平移操作,使x0,x1的极值重合
Max0 = max(x0);
Max1 = max(x1);
err = Max0 - Max1;
extrx1 = x1 + err;
x2 = extrx1;
% 对极值做匹配度的计算
tol = 0;
for i = 1:len0
    error1(i) = x0(i)-x2(i);
    tol = tol+(x0(i)-x2(i))*(x0(i)-x2(i));
end;
avg = (sum(abs(x0))+sum(abs(x2)))/(len0*2-2);
erravg = sum(abs(error1))/(len0-1);

cleverblue

一种自适应的EMD端点延拓方法
邵晨曦 ， ， t-，王 剑 ，范金锋 ，杨 明 ，王子才
(1．1}1国科学技术大学多媒体计算与通信教育部．微软重点实验室，安徽合肥230026；
2．中【日科学技l术大学计算机科学技术系，安徽合肥230o26；3．哈尔滨T业大学控制 j仿
真t{，心，黑龙江哈尔滨150001；4 安徽省计算与通讯软件重点实验室，安徽合肥230027)
摘要： 由美国国家航空航天局(NASA)的Huang等发明的经验模态分解(EMD)是一种先进的信号处理方法，能
够有效地获得非平稳信号的时频特征，但是其利用样条曲线构造信号上下包络线的过程中存在严重的端点问题．在研
究了该问题已有方法的基础上，提出1r一种基于波形匹配的自适应端点延拓方法，采用信号内部和端点处变化趋势最
为相似的子波来对端点处的信号进行延拓．该方法充分考虑了信号的内在特性以及边缘处的变化趋势，使端点处的延
拓更加合理，从而使得三次样条曲线在端点处不会发生大的摆动 实验表明该方法能够有效地抑制端点效应