function QRSref=QRSdetectorM(vadx,vdx,Fs,pth,RRts,pmQT)
% --------------------------------------------------------------------------------------------
% QRSdetectorM.m: QRS detector
%
% Input parameters:
% vadx : array of filtered absolute derivate values
% vdx : array of filtered derivate values
% Fs : sampling frequency
% pth : threshold on derivative
% RRts : RR (seconds) typical of the animal
% pmQt : fraction of QT length for QT maks
%
% Ouput parameters:
% qrsRef : array of QRS reference points (max signed derivative)
%
% Example:
% qrsM=QRSdetectorM(vadx,vdx,fs,0.4,0.86,1);
%
% NOTE:
% Decreasing the "pth" threshold leads to sensibility increasing but specificity decreasing.
% QT mask is proportional to the input typical RR.
% Bazzet formula, which involve sqrt(RRs), is a quite good approximation of QT length for humans,
% it is not valid across animal species
%
% 2002 Matlab translation of the "didactic" version of QRS detector developed in Mathcad.
% 2007 QRS fiducial point as max signed derivative
% 2008/11 Modified control row 119 (if(i < iinizio+QRSd || i<ifine+nsd))
%
% --------------------------------------------------------------------------------------------
% Author: Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% --------------------------------------------------------------------------------------------
% 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.
% --------------------------------------------------------------------------------------------
if (nargin < 3),
error('At least 3 parameters are required'); return
end
if (nargin < 4), pth=0.5; end
if (nargin < 5), RRts=0.86; end % 0.86=human % 0.25= mouse
if (nargin < 6), pmQT=1; end
sqRRts=sqrt(RRts);
QTlen=0.420*sqRRts; % normal Qt length (RR=1 => QT=0.42s, humans)
%QRSd=round(0.15*RRtc); % extension of the previous QRS detection is allowed in the interval (ifine, iinizio+QRSd)
QRSd=round((0.05+0.25*sqRRts)*Fs); % extension of the previous QRS detection is allowed in the interval (ifine, iinizio+QRSd)
maskQT=round((0.07+pmQT*QTlen)*Fs); % QT mask length (RR=1s => QT=0.42s, humans), QRS detection threshold decreases (linearly)
rthT=2.4; % from rthT*th if i=ifine to th if i=ifine+maskQt
% nsp=round(0.05*sqRRts*Fs); % max number of sample before threshold crossing (increasing)
nsp=round(0.15*sqRRts*Fs); % max number of sample before threshold crossing (increasing) 23-09-05
nsd=round(0.078*sqRRts*Fs); % max number of sample after threshold crossing (decreasing )
RRtc=Fs*RRts;
nqappr=fix(1.2*length(vadx)/RRtc); % approximate estimate of number of QRSs
QRSref=zeros(nqappr,1); inizio=zeros(nqappr,1);
vimaxd=zeros(nqappr,1); fine=zeros(nqappr,1);
isai=1.5*Fs; % index of first sample used in inizialization
fsai=min(length(vadx), floor(60*RRtc)); % index of last sample used in inizialization
w2=fix(2*RRtc); % wide windows containing at least one QRS
mD2=meanMaxSc(vadx(isai:fsai), w2, 1,1); % compute the average of maximum derivatives on windows of 2s
% (1% of minima and 1% of maxima are discard)
meaD=mD2;
th=pth*meaD;
% --- choose derivative signum for fiducial QRS pointer
%[minD, maxD] =mimaxsc(vdx(1:nsai),nsc,nsc);
[minD, maxD] =mimaxsc(vdx(isai:fsai),1,1);
if(maxD > -minD*1.1), vsdx=vdx; else vsdx=-vdx; end
jq=1;
QRS=0; % flag asserting a previous QRS threshold overcoming
maxd=0;
imaxd=1;
RRcm=RRtc;
iinizio=-maskQT;
ifine=iinizio;
for i=1:length(vadx) % main loop on derivative samples
vadxi=vadx(i);
if(QRS) % inside QRS interval
if(vadxi < th) % absolute derivative becomes lower than threshold
inizio(jq)=iinizio; % save the index of last threshold crossing
vimaxd(jq)=imaxd; % save the index of max derivative
ifine=i;
fine(jq)=i;
maxdc=min(maxd,th*4); % bounding of derivative maximum
% [maxs,imaxs]=max(vsdx(max(iinizio-nsp,1):min(ifine+nsd,length(vsdx))));
[maxs,imaxs]=max(vsdx(max(iinizio,1):min(ifine,length(vsdx))));
QRSref(jq)=max(iinizio,1)-1+imaxs;
meaD=0.97*meaD + (1-0.97)* maxdc; % updating the average of derivative maximum
th= pth*meaD; % updating the threshold on derivative
if(jq>1)
RRcj=QRSref(jq)-QRSref(jq-1);
RRcm=RRcm+ 0.97* sign(RRcj-RRcm)*min(abs(RRcj-RRcm), 0.1*RRcm) ;
RRsm=RRcm/Fs;
if(RRsm<RRts*.4 || RRsm>RRts*2.5), RRsm=RRts; RRcm=RRsm*Fs; end
sqRRsm=sqrt(RRsm);
QTlen=0.420*sqRRsm;
maskQT=round((0.07+pmQT*QTlen)*Fs);
nsp=round(0.15*sqRRsm*Fs);
nsd=round(0.078*sqRRsm*Fs);
end
QRS=0;
jq=jq+1;
elseif(vadxi > maxd)
imaxd=i; % save the index of max derivative
maxd=vadxi;
end
elseif(vadxi > th) % absolute derivative greater than threshold
if(i < iinizio+QRSd || i<ifine+nsd) % extend previous QRS detection
if(jq>1), jq=jq-1; end
ifine=i;
if(vadxi > maxd)
imaxd=i;
maxd=vadxi;
end
QRS=1; % set again QRS flag, previous QRS was not ended
% To avoid T wave the QRS detection threshold decreases (linearly) from "rthT*th" to "th" for
% i>=ifine and i<ifine+maskQt
elseif(vadxi > th*(rthT - (rthT-1)*(i-ifine)/maskQT))
imaxd=i;
iinizio=i; % save index of the first threshold crossing
maxd=vadxi;
QRS=1; % set flag, a new QRS is detected
end
end
end
if(jq>1)
QRSref(jq:end)=[];
else
QRSref=[];
end
end % = function ===============================================================