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