function sqi = bsqi(refqrs,testqrs,thres,margin,windowlen,fs)
% sqi = bsqi(refqrs,testqrs,thres,margin,windowlen,fs)
% compare two sets of annotation with one as the reference (refqrs) and one
% as the test (testqrs)
%
% inputs
% refqrs: reference qrs annotation (in sec)
% testqrs: test qrs annotations (in sec)
% thres: threshold (in sec,default 0.05s)
% margin: margin time not include in comparison (in sec,default 2s)
% windowlen: length of the comparison window (in sec,default 60s)
%
% output
% sqi: match proportion according to some criteria you can change
% depending on what you are looking for (can be Se, PPV or F1 measure).
% See at the end of the function.
%
% Reference: inspired from
% Li, Qiao, Roger G. Mark, and Gari D. Clifford. "Robust heart rate estimation
% from multiple asynchronous noisy sources using signal quality indices and
% a Kalman filter." Physiological measurement 29.1 (2008): 15.
%
%
% FECG extraction toolbox, version 1.0, Sept 2013
% Released under the GNU General Public License
%
% Copyright (C) 2013 Joachim Behar
% Oxford univseristy, Intelligent Patient Monitoring Group - Oxford 2013
% joachim.behar@eng.ox.ac.uk
%
% Last updated : 02-08-2013
%
% 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 in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.
% == managing inputs
if nargin<2; error('bsqi: wrong number of input arguments \n'); end;
if nargin<3; thres=0.05; end;
if nargin<4; margin=2; end;
if nargin<5; windowlen=60; end;
if nargin<6; fs=1000; end;
start = margin*fs;
stop = (windowlen-margin)*fs;
refqrs = refqrs*fs;
testqrs = testqrs*fs;
try
refqrs = refqrs(refqrs>start & refqrs<stop)'; % reference annotations
testqrs = testqrs(testqrs>start & testqrs<stop)'; % test annotations
NB_REF = length(refqrs);
NB_TEST = length(testqrs);
% == core function
[IndMatch,Dist] = dsearchn(refqrs,testqrs); % closest ref for each point in test qrs
IndMatchInWindow = IndMatch(Dist<thres*fs); % keep only the ones within a certain window
NB_MATCH_UNIQUE = length(unique(IndMatchInWindow)); % how many unique matching
TP = NB_MATCH_UNIQUE; % number of identified ref QRS
FN = NB_REF-TP; % number of missed ref QRS
FP = NB_TEST-TP; % how many extra detection?
Se = TP/(TP+FN);
PPV = TP/(FP+TP);
F1 = 2*Se*PPV/(Se+PPV); % accuracy measure
sqi = PPV; % EDIT ME: for my application PPV better catches what we want but you can build your bsqi based on Se,PPV or F1 measure
catch ME
for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
sqi = 0;
end
end