[m n]=size(data);
Lead=1;
d=abs(resample(diff(resample(data(:,Lead),1,10)),10,1));
hh=sort(d);
md=int32(m*5.5/6);
MeanD=hh(md);
j=0;
i=200;
while i<m-200
    while (d(i)<MeanD) && (i<m-200)
        i=i+1;
    end
    j=j+1;
    Rp(j)=i;
    i=i+200;
end
Rcount=j-1;

rint=80; %QRS interval length
Rstart=-20; %T-wave interval start
ModelR(1:rint,1)=0;
for j=1:5
    for ii=1:NrofLead
        for i=1:rint
           ModelR(i,1)=ModelR(i,1)+((data((Rp(j)+i+Rstart),ii)-data((Rp(j)+1+Rstart),ii))^2)/5;
        end;
    end;
end;

for j=1:Rcount
    CurrR(1:rint,1)=0;
    for ii=1:NrofLead
        for i=1:rint
           CurrR(i,1)=CurrR(i,1)+((data((Rp(j)+i+Rstart),ii)-data((Rp(j)+1+Rstart),ii))^2);
        end;
    end;
    Rtest(:,1)=CurrR(:,1);
    Rtest(:,2)=ModelR(:,1);
    v = xcov(CurrR,ModelR);
    [C,I]=max(v);
    R(j)=Rp(j)+I-rint;

    CurrR(1:rint,1)=0;
    for ii=1:NrofLead
        for i=1:rint
           CurrR(i,1)=CurrR(i,1)+((data((R(j)+i+Rstart),ii)-data((R(j)+1+Rstart),ii))^2);
        end;
    end;
    Rtest(:,3)=CurrR(:,1);
    ModelR=ModelR/2+CurrR/2;
    ModelRt(:,j)=ModelR;
       
end;