tryemdt_h4.m

EMD分解（利用在极值点序列首尾添加极值来处理问题）

%EMD的分解过程

% dt=0.001;
% T=1;
% n=T/dt;
% t=0:dt:T;
% xt=2*sin(2*pi*15*t)+4*sin(2*pi*10*t).*sin(2*pi*t/10)+sin(2*pi*5*t);  %原始信号

function imf=TryEmdt(h1,t)
  flag=0;
  ll=0;
  xtT=h1;
 while max(abs(h1))>0.5
     uu=0;
     while(1)
      uu=uu+1
      [bg,bt,cg,ct,N2]=MaxMin(h1,t);
       l1=length(bt);
       l2=length(ct);
   
       if (l1==1)||(l2==1)
         flag=1;                     %只有一个极值点的情况 
       break
       end
   
        N1=ZeroNum(h1);
        [By,Cy]= Enovelope(bt,bg,ct,cg,t);
        Avg=(By+Cy)/2;                   %极大极小平均
        h=h1-Avg;
       if ((abs(N2-N1)<=1)&&(max(abs(Avg))<0.001))||(uu>800)
        ll=ll+1
        imf(ll,:)=h1;
        break
       else
        h1=h;
       end
     end

    if (flag==1)
       flag=0;
       break
    else
     xtt=xtT-imf(ll,:);
     h1=xtt;
     xtT=h1;
    end
 end
 imf(ll+1,:)=h1;    %残余模态
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%极大极小值包络线 
%前后两端各添加一个极大极小值（采用平均的方法）
function [By,Cy]= Enovelope(bt,bg,ct,cg,t)
mxt1=t(1);
mxt2=t(end);
mit1=t(1);
mit2=t(end);
mxs1=bg(1);
mxs2=bg(end);
mis1=cg(1);
mis2=cg(end);
btt=[mxt1 bt mxt2];
ctt=[mit1 ct mit2];
bgg=[mxs1 bg mxs2];
cgg=[mis1 cg mis2];

By=spline(btt,bgg,t);              %极大值三次样条曲线
Cy=spline(ctt,cgg,t);              %极小值三次样条曲线

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function N1=ZeroNum(h1)
%--ZeroNum函数
%--求过零点个数
N1=0;
N=length(h1);
for i=1:N-1
    if h1(i)==0
        N1=N1+1;
    elseif h1(i)*h1(i+1)<0
       N1=N1+1;    %零点数   
    end
end
if h1(N)==0
    N1=N1+1;
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bg,bt,cg,ct,N2]=MaxMin(h1,t)
%--MaxMin函数
%--求局部极值点数
max=0;
max_No=0;
min=0;
min_No=0;
m=0;
k=0;
jj=0;
N=length(h1);
%  if h1(1)~=0
%        if (h1(1)>h1(2))
%            m=m+1;
%         bg(m)=h1(1);
%         bt(m)=t(1);
%        else 
%            k=k+1;
%         cg(k)=h1(1);
%         ct(k)=t(1);
%        end
%     end

for i=2:N-1
    if (h1(i)>=h1(i-1))&(h1(i)>=h1(i+1))
         m=m+1;
       bg(m)=h1(i);                   %极大值点
       bt(m)=t(i);
    end
    if (h1(i)<=h1(i-1))&(h1(i)<=h1(i+1))
         k=k+1;
       cg(k)=h1(i);                   %极小值点
       ct(k)=t(i);
    end
end
%   if h1(N)~=0
%     if h1(N)>h1(N)
%         m=m+1;
%         bg(m)=h1(N);
%         bt(m)=t(N);
%     else 
%         k=k+1;
%         cg(k)=h1(N);
%         ct(k)=t(N);
%     end
%    end
N2=k+m;                   %极值点数