function [F1,Se,PPV,Nb] = run_sqi(refqrs,testqrs,thres,margin,windowlen,fs)
% sqi = run_sqi(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)
% fs: sampling frequency
%
% 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.
%
% When using this work, then please cite [1] and [2]:
% [1] Behar Joachim, Oster Julien, Qiao Li, Clifford Gari D. Signal Quality
% During Arrhythmia and its Application to False Alarm Reduction.
% IEEE Transactions on Biomedical Engineering. 60(6). 1660-6. 2013.
%
% [2] 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.
%
% PCinCC2014, version 1.0, June 2014
% Released under the GNU General Public License
%
% Copyright (C) 2014 Joachim Behar
% Oxford university, Intelligent Patient Monitoring Group - Oxford 2014
% joachim.behar@gmail.com
%
% Updates:
% 07-02-2014
% JB- testes with Octave -> running OK
%
% 02-10-2014
% Bug fix: Dealing with annotations close to the border of the search
% window. Lines 70 - 100
% David Springer
%
% 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;
if size(refqrs,1)>size(refqrs,2); refqrs=refqrs';end
if size(testqrs,1)>size(testqrs,2); testqrs=testqrs';end
start = margin*fs;
stop = (windowlen-margin)*fs;
refqrs = refqrs*fs; % convert into samples from time
testqrs = testqrs*fs; % convert into samples from time
try
refqrs = refqrs(refqrs>start & refqrs<stop)'; % reference annotations
testqrs = testqrs(testqrs>start & testqrs<stop)'; % test annotations
if ~isempty(refqrs)
NB_REF = length(refqrs);
NB_TEST = length(testqrs);
% == removing borders problems
indbord = find(refqrs<thres*fs | refqrs>(windowlen-thres)*fs); % reference QRS at the border
if ~isempty(indbord)
[IndMatchBord,DistQRSbord] = dsearchn(testqrs,refqrs(indbord));
Indeces_below_threshold = DistQRSbord<thres*fs; %Added line to find indices of annotation of interest (refqrs) below threshold (02-10-14)
IndMatchBord = IndMatchBord(Indeces_below_threshold);
NB_QRS_BORD = length(indbord); % QRS at the border of the window (0,1,2)
NB_MATCHING = length(IndMatchBord); % nb of corresponding mathing QRS (0,1,2)
if isempty(IndMatchBord);
refqrs(indbord) = [];
elseif NB_MATCHING<NB_QRS_BORD
refqrs(indbord(~Indeces_below_threshold)) = []; %Removing other indices not below threshold (02-10-14)
end
end
indbord = find(testqrs<thres*fs | testqrs>(windowlen-thres)*fs); % reference QRS at the border
if ~isempty(indbord)
[IndMatchBord,DistQRSbord] = dsearchn(refqrs,testqrs(indbord));
Indeces_below_threshold = DistQRSbord<thres*fs; %Added line to find indices of annotation of interest (testqrs) below threshold (02-10-14)
IndMatchBord = IndMatchBord(Indeces_below_threshold);
NB_QRS_BORD = length(indbord); % QRS at the border of the window (0,1,2)
NB_MATCHING = length(IndMatchBord); % nb of corresponding mathing QRS (0,1,2)
if isempty(IndMatchBord);
testqrs(indbord) = [];
elseif NB_MATCHING<NB_QRS_BORD
testqrs(indbord(~Indeces_below_threshold)) = []; %Removing other indices not below threshold (02-10-14)
end
end
% == 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
Nb.TP = TP;
Nb.FN = FN;
Nb.FP = FP;
else
F1=[];Se=[];PPV=[];Nb=[];
end
catch ME
F1=[];Se=[];PPV=[];Nb=[];
end
end