function  PeakIndex=peakdetection(wave,fs,dilation,thrhd,Loop)
%thrhd=0.5;
%Loop=5;

%输入参数说明
%wave为要检r波的信号
%fs为采样率
%dilation为膨胀系数，一般定义为0.1
%thrhd为阈值
%Loop为分段数，定义为5
%输出参数说明
%PeakIndex为差分信号的R波位置

NL=length(wave);
NLs=NL/Loop;
maxv=zeros(1,Loop);
for I=1:Loop;
    Indexs=round((I-1)*NLs+1:I*NLs);
    temp=wave(Indexs);
    maxv(I)=max(temp);
end
maxm=min(maxv);
Vector01=wave>=(maxm*thrhd);
 temp=[0,Vector01];
Vector10N1=Vector01-temp(1:NL);

%%%%%%%% dilation
II=1;
while 0<1;
   IndexN1=find(Vector10N1(II:end)==-1,1,'first');
   if isempty(IndexN1)
      break;
   else
      Index1=find(Vector10N1(II+IndexN1:end)==1,1,'first');
      if isempty(Index1)
          break;
      else
          if(Index1>fs*dilation)
              II=II+IndexN1+Index1;
              continue;
          else
             Vector10N1([II+IndexN1-1,II+IndexN1+Index1-1])=[0,0];
             II=II+IndexN1+Index1;
          end        
      end
   end
end

%%%%%%%%%%%%%%%%%%%%% Initial Peakdetection
PeakIndex=[];
II=1;
while 0<1;
   Index1=find(Vector10N1(II:end)==1,1,'first');
   if isempty(Index1)
      break;
   else
      IndexN1=find(Vector10N1(II+Index1:end)==-1,1,'first');
      if isempty(IndexN1)
          [M,temp]=max(wave(II+Index1-1:end));
          PeakIndex=[PeakIndex,II+Index1+temp-2];
          break;
      else
          [M,temp]=max(wave(II+Index1-1:II+Index1+IndexN1-1));
          PeakIndex=[PeakIndex,II+Index1+temp-2];
          II=II+Index1+IndexN1;
      end
   end
end
P2PInterval=PeakIndex(2:end)-PeakIndex(1:end-1);
MInterval=median(P2PInterval);
[M,MIndex]=min(abs(P2PInterval-MInterval));

%%%%%%%%%%%%%%%%%%补漏
PeaksRight0=PeakIndex(MIndex+1:end);
PeaksRightD=PeaksRight0;
II=1;
while 0<1;
    if length(PeaksRightD)>=2;
       if PeaksRightD(2)-PeaksRightD(1)<1.5*MInterval;
          PeaksRightD=PeaksRightD(2:end);                  
       else
          L=PeaksRightD(1)+round(0.8*MInterval);
          R=PeaksRightD(1)+round(1.2*MInterval);
          temp=wave(L:R);
          [M,Im]=max(temp);
          ImReal=L+Im-1;
          PeaksRightD(1)=ImReal;
          PeaksRight0=[PeaksRight0(1:II),ImReal,PeaksRight0(II+1:end)];          
       end
       II=II+1;
    else
        break;
    end

end
PeaksRight=PeaksRight0;
%%
PeaksLeft0=PeakIndex(1:MIndex);
PeaksRight0=fliplr(PeaksLeft0);
PeaksRightD=PeaksRight0;
II=1;
while 0<1;
    if length(PeaksRightD)>=2;
       if PeaksRightD(1)-PeaksRightD(2)<1.5*MInterval;
          PeaksRightD=PeaksRightD(2:end);          
       else
          L=PeaksRightD(1)-round(1.2*MInterval);
          R=PeaksRightD(1)-round(0.8*MInterval);
          temp=wave(L:R);
          [M,Im]=max(temp);
          ImReal=L+Im-1;
          PeaksRightD(1)=ImReal;
          PeaksRight0=[PeaksRight0(1:II),ImReal,PeaksRight0(II+1:end)];
       end
       II=II+1;
    else
       break;
    end
end
PeaksLeft=fliplr(PeaksRight0);
%%
PeakIndex=[PeaksLeft,PeaksRight];
P2PInterval=PeakIndex(2:end)-PeakIndex(1:end-1);
MInterval=median(P2PInterval);
[M,MIndex]=min(abs(P2PInterval-MInterval));

%%%%%%%%%%%%%%%%%%%%%%%%%%%去伪
PeaksRight0=PeakIndex(MIndex+1:end);
PeaksRightD=PeaksRight0;
II=1;
while 0<1;
    if length(PeaksRightD)>=3; 
       a=(PeaksRightD(3)-PeaksRightD(1))<=1.3*MInterval;
       b=(PeaksRightD(2)-PeaksRightD(1))<=0.8*MInterval;
       if a&b   
          Ix=find(PeaksRight0==PeaksRightD(2));
          PeaksRight0=[PeaksRight0(1:Ix-1),PeaksRight0(Ix+1:end)];
          PeaksRightD=PeaksRightD(3:end);
       else
          PeaksRightD=PeaksRightD(2:end); 
       end       
    else
       break;
    end
    II=II+1;
end
PeaksRight=PeaksRight0;

%%
PeaksLeft0=PeakIndex(1:MIndex);
PeaksRight0=fliplr(PeaksLeft0);
PeaksRightD=PeaksRight0;
II=1;
while 0<1;
    if length(PeaksRightD)>=3;
       a=PeaksRightD(1)-PeaksRightD(3)<=1.3*MInterval;
       b=PeaksRightD(1)-PeaksRightD(2)<=0.8*MInterval;
       if a&b   
          Ix=find(PeaksRight0==PeaksRightD(2));
          PeaksRight0=[PeaksRight0(1:Ix-1),PeaksRight0(Ix+1:end)];
          PeaksRightD=PeaksRightD(3:end);
       else
          PeaksRightD=PeaksRightD(2:end); 
       end       
    else
       break;
    end
    II=II+1;
end
PeaksLeft=fliplr(PeaksRight0);
%%
PeakIndex=[PeaksLeft,PeaksRight];
P2PInterval=PeakIndex(2:end)-PeakIndex(1:end-1);
MInterval=median(P2PInterval);
[M,MIndex]=min(abs(P2PInterval-MInterval));

% %%%%%%%%%%%%%%%%%%%%首尾  %%暂时不用8.3
temp=min(P2PInterval(2:end-1));
if P2PInterval(1)<=0.95*temp;
    PeakIndex=PeakIndex(2:end);
end
if P2PInterval(end)<=0.95*temp;
    PeakIndex=PeakIndex(1:end-1);
end


end