function QRSref=QRSdetectorF2(vadx,vdx,Fs,pth,RRst,RRsi,lensi,fb,pmQT)
% --------------------------------------------------------------------------------------------
% QRSdetectorF2.m: QRS detector
% QRS fiducial point as max signed derivative
% It needs to start from a controlled signal interval
%
% Input parameters:
% vadx : array of filtered absolute derivate values
% vdx : array of filtered derivate values
% Fs : sampling frequency
% pth : threshold on derivative
% RRst : RR (seconds) typical of the animal
% RRsi : RR (seconds) initial value
% lensi: initialization interval length (seconds)
% fb : forward (0) - backward (1) direction flag
% pmQt : fraction of QT length for QT mask (NOT used)
%
% Ouput parameters:
% QRSref (reference point, max signed derivative )
%
% Example:
% qrsM=QRSdetectorF2(vadx,vdx,fs,0.4,0.86,1);
% --------------------------------------------------------------------------------------------
% 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), RRst=0.86; end % 0.86=humans % 0.25= mice
if (nargin < 6), RRsi=RRst; end
if (nargin < 7), lensi=RRsi*10; end
if (nargin < 8), fb=0; end % 0=forward, 1=backward direction
if (nargin < 9), pmQT=1; end
RRtc=Fs*RRst;
RRci=Fs*RRsi;
lenci=lensi*Fs;
% sqRRst=sqrt(RRst);
% QTlen=0.420*sqRRst; % normal Qt length (RR=1 => QT=0.42s, humans)
m=5; mu=0.01; alp=mu/10; % LMS algorithm parameters (RR adaptive predictive model)
nqrsMap=ceil(length(vdx)/RRci)*2; % used only for vector preallocation
QRSref=zeros(nqrsMap,1);
RRcVz=zeros(nqrsMap,1);
%maskQT=round((0.07+pmQT*QTlen)*Fs); % QT mask length (RR=1s => QT=0.42s, humans)
wleftAmi=0.2; wrightAmi=0.2;
wleft= round(0.5*RRst*Fs); wright= round(1.2*RRst*Fs);
wleftc= round(0.15*RRst*Fs); wrightc= round(0.15*RRst*Fs);
isai=1; % index of first sample used in inizialization
fsai=min(length(vadx),round(max(4*RRci, lenci))); % index of last sample used in inizialization
w2=fix(2*RRci); % length of window containing at least one QRS
mD2=meanMaxSc(vadx(isai:fsai), w2, 0.5,0.5); % compute the average of maximum derivatives on windows of 2*RRci
% (0.5% of minima and 0.5% of maxima are discarded)
meaD=mD2;
th=pth*meaD; % threshold on derivative
% --- choose derivative signum for fiducial QRS pointer
[minD, maxD] =mimaxsc(vdx(isai:fsai),1,1);
if(fb)
if(minD < -maxD*1.1), vsdx=-vdx; else vsdx=vdx; end % backward
else
if(maxD > -minD*1.1), vsdx=vdx; else vsdx=-vdx; end % forward
end
j=1;
RRsm=RRsi; RRcm=RRsm*Fs; RRest=RRcm;
QRSrefj=1;
while QRSrefj<=length(vdx) % main loop on derivative samples
QRSrefjlmi=max(QRSrefj-wleft,1); QRSrefjrma=min(QRSrefj+wright,length(vsdx));
if(j>1)
%---
weightw=ones(wleft+wright+1,1);
for k=wrightc:wright
if(weightw(wleft+k)-0.02>wrightAmi), weightw(wleft+k+1)= weightw(wleft+k)- 0.002;
else weightw(wleft+k+1)= weightw(wleft+k); end
end
for k=wleftc:wleft
if(weightw(wleft-k+2)-0.02>wleftAmi), weightw(wleft-k+1)= weightw(wleft-k+2)- 0.002;
else weightw(wleft-k+1)= weightw(wleft-k+2); end
end
weightw=weightw(1:QRSrefjrma-QRSrefjlmi+1);
%----
[maxs,imaxs]=max(weightw.*vsdx(QRSrefjlmi: QRSrefjrma)); % max of derivative
else
%QRSrefjrma=min(QRSrefj+round(RRcm),length(vsdx));
weightw=ones(wleft+wright+1,1);
for k=wrightc:wright
if(weightw(wleft+k)-0.02>wrightAmi), weightw(wleft+k+1)= weightw(wleft+k)- 0.002;
else weightw(wleft+k+1)= weightw(wleft+k); end
end
weightw=weightw(1:QRSrefjrma-QRSrefjlmi+1);
[maxs,imaxs]=max(weightw.*vsdx(QRSrefjlmi: QRSrefjrma)); % max of derivative
end
if(maxs>th*weightw(imaxs)) % comparison of max of derivative with weighted threshold
QRSrefj=QRSrefjlmi-1+imaxs; % position of the detected QRS
QRSref(j)=QRSrefj; % QRS pointer saving
if(j>1)
RRcj=QRSref(j)-QRSref(j-1);
RRcmp=RRcm; RRczj=RRcj-RRcmp;
RRcVz(j,1)=RRczj;
RRcm=RRcmp+ 0.05* sign(RRczj)*min(abs(RRczj), 0.05*RRcmp);
RRsm=RRcm/Fs;
if(RRsm<RRst*.5 || RRsm>RRst*2.5), RRsm=RRst; 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);
wleft = round(0.3*RRsm*Fs); wright = round(0.8*RRsm*Fs);
% wleftc= round(0.15*RRsm*Fs); wrightc= round(0.15*RRsm*Fs);
wleftc= round(0.02*RRsm*Fs); wrightc= round(0.02*RRsm*Fs);
RRest=RRcm;
if(j==m+1)
RRcVzj=RRcVz(2:m+1);
a = lpc(RRcVzj,m); w=-a(end:-1:2)'; % initializing predictive filter weigths
RRczest=w'*RRcVzj; % zero mean estimated RR value
RRest=RRcmp+RRczest; % estimated RR value
sz2=RRcVzj'*RRcVzj/m; % initializing average power
%RRest = filter([0 -a(2:end)],1,RRcVzj);
elseif(j>m+1)
sz2 = sz2 + alp*sign(RRczj*RRczj-sz2)*min(RRczj*RRczj-sz2, sz2); % average power updating
cp=0.01*sz2;
er=RRczj-RRczest; % predictive error
RRcVzj=RRcVz(j-m+1:j);
if(er*er < 0.5*sz2)
w = w + (2*mu/(cp+RRcVzj'*RRcVzj))*er*RRcVzj; % weigth updating by LMS algorithm
else
w=0.97*w;
end
RRczest=w'*RRcVzj; % zero mean estimated RR value
if(abs(RRczest)>0.03*RRcmp)
RRczest=sign(RRczest)*0.03*RRcmp;
end
RRest=RRcmp+RRczest; % estimated RR value
end
end
QRSrefj=QRSrefj+round(RRest); % estimated next QRS position
j=j+1; % increasing the index of qrs pointer vector
else
QRSrefj=QRSrefjlmi+imaxs+wleft;
end
end
if(j>1)
QRSref=QRSref(1:j-1);
else
QRSref=[];
end
end %== function ================================================================
%