% -------------------------------------------------------------------------------------------------
% QRSdetector.m: QRS detector
%   Translation to Matlab of the Mathcad "didactic" version
%
%   Input parameters:
%	 vadx : array of filtered absolute derivate values  
%	 Fs   : sampling frequency (in Hertz).  
%	 pth  : threshold 
%
%   Ouput parameters:
%   qrsM  : matrix of 3 rows:
%           QRS supThDer
%           QRS maxDer, 
%           QRS infThDer
%   
%	Example: 
%	 qrsV=QRSdetector(vadx,fs,0.6); 
%
% Copyright (C) 2004 Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
%
% This program is free software; you can redistribute it and/or modify it under the terms
% of the GNU General Public License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version.
%
% This program is distributed "as is" and "as available" in the hope that it will be useful,
% but WITHOUT ANY WARRANTY of any kind; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program;
% if not, write to the Free Software Foundation, Inc.,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% -------------------------------------------------------------------------------------------------

function qrsM=QRSdetector(vadx,Fs,pth)

if (nargin < 2),
    error ( 'At least 2 parameters are required' ); 
elseif (nargin == 2)
    pth=0.5; 
end

QRSd=0.15*Fs;
maskQT=0.4*Fs;

w02=fix(0.2*Fs);
mD02=meanMaxSc(vadx, w02, 1,1);

w2=fix(2*Fs);
mD2=meanMaxSc(vadx, w2, 1,1);

meaD=mD2;
th=pth*meaD;
j=1;
QRS=0;
maxd=0;
imaxd=1;
iinizio=1;
%ifine=maskQT;
ifine=iinizio;
for i=1:length(vadx)
    value=vadx(i);
    if(QRS)
        if(value < th)
            inizio(j)=iinizio;
            vimaxd(j)=imaxd;
            ifine=i;
            fine(j)=i;
            maxdc=min(maxd,th*4);
            meaD=0.97*meaD + (1-0.97)* maxdc;
            th= pth*meaD;
            QRS=0;
            j=j+1;
        elseif(value > maxd)
            imaxd=i;
            maxd=value;
        end
    elseif(value > th)
        if(i > ifine+maskQT)|(value > th*(2-(i-(ifine+QRSd))/maskQT))
            imaxd=i;
            iinizio=i;
            maxd=value;
            QRS=1;
        elseif(i < ifine+QRSd)
            if(j>1) j=j-1; end
            if(value > maxd)
                imaxd=i;
                maxd=value;
            end
            QRS=1;
        end
    end
end
if(j>1)
    qrsM=[inizio; vimaxd; fine];
else qrsM=[];
end