%% Estimate the quality of QRS complex detections
% Calculate the quality of QRS detections based on the the likelihood of
% measurements from RR interval series with respect to a trained model. See
% for further details. See 'A Pattern-Recognition Approach for
% Lead-Selection in Heartbeat Detection' CINC 2014.
%
% Example
%
% [ ratio, estimated_labs ] = CalcRRserieRatio(time_serie, ECG_header, start_end, ECG_type)
%
% Arguments:
% +time_serie: [cell] REQUIRED
%
% Cell array of size [nsig 1] with the detections
% performed for each lead.
%
% +ECG_header: [struct] OPTIONAL.
%
% Description of the ECG typically available in the
% ECG_header. Structure with fields:
%
% -freq: Sampling rate in Hz. (1)
%
% -nsig: Number of ECG leads. (size(ECG,2))
%
% -nsamp: Number of ECG samples. (size(ECG,1))
%
% +start_end: [numeric] OPTIONAL. Vector of size [2 1] with the start
% and end samples.
%
% +ECG_type: [numeric] OPTIONAL. Char indicating if you know something
% about the ECG recording type a priori. Possible
% values are: 'stress' 'arrhythmia' 'longterm' 'sinus'
%
% Output:
% + ratio: a measure of quality of the QRS detections
%
% + estimated_labs: a label for each heartbeat: {TP, FP and FN}
%
% See also ECGtask_QRS_detection
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate : 23/4/2013
% Copyright 2008-2015
function [ ratio, estimated_labs ] = CalcRRserieRatio(time_serie, ECG_header, start_end, ECG_type)
%% constants
min_pattern_separation = 350; % ms
max_pattern_separation = 1500; % ms
cECGtypes = {'stress' 'arrhythmia' 'longterm' 'sinus'};
%% start
if( nargin < 2 || isempty(time_serie) || isempty(ECG_header) )
ratio = 0;
return
end
if( nargin < 3 || isempty(start_end) )
start_end = [1 ECG_header.nsamp];
end
if( nargin < 4 || ~any(strcmpi(ECG_type, cECGtypes)) )
ECG_type = 'stress';
end
start_sample = start_end(1);
end_sample = start_end(2);
pb = progress_bar( 'Calculating ratios', 'Start' );
% co_ocurrences respect other leads
co_ocurrence = calc_co_ocurrences(time_serie);
if( isempty(co_ocurrence) )
return
end
dummy = prdataset([]);
dummy = prmapping([]);
% load the trained classifier.
aux_load = load([ 'QRS_q_' ECG_type '.mat']);
w_lablist = cellstr(getlabels(aux_load.wTrained_Classifier));
FN_lab = find(strcmpi(w_lablist, 'FN'));
FP_lab = find(strcmpi(w_lablist, 'FP'));
TP_lab = find(strcmpi(w_lablist, 'TP'));
lreferences = length(time_serie);
pb.Loops2Do = lreferences;
pb.checkpoint('Calculating ratios.');
ratio = nan(lreferences,1);
estimated_labs = cell(lreferences,1);
for ii = 1:lreferences
pb.start_loop();
this_ts = colvec(time_serie{ii});
if( length(this_ts) < (ECG_header.nsamp / max_pattern_separation) )
% warning( [mfilename ':few_anns'], 'Few detections in %d\n', ii );
ratio(ii) = 0;
else
pb.checkpoint([]);
k_gaps = calc_k_gaps(this_ts);
% build the feature matrix [ RR_i-1 RR_i RR_10 RR_60 co_ocurrences]
% RR_i-1
this_ts = round(this_ts * 1000 / ECG_header.freq);
pb.checkpoint([]);
aux_fm = [ calculate_RR_features(this_ts) colvec(co_ocurrence{ii}) ];
pb.checkpoint([]);
ds_result = prdataset(aux_fm) * aux_load.wTrained_Classifier;
estimated_labs{ii} = renumlab(ds_result * labeld, char(w_lablist) );
aux_val = sum(estimated_labs{ii} == TP_lab | estimated_labs{ii} == FN_lab);
if( aux_val == 0 )
this_se = 0;
else
this_se = sum(estimated_labs{ii} == TP_lab) / aux_val;
end
aux_val = sum(estimated_labs{ii} == TP_lab | estimated_labs{ii} == FP_lab);
if( aux_val == 0 )
this_pp = 0;
else
this_pp = sum(estimated_labs{ii} == TP_lab) / aux_val;
end
this_q = (2*this_se + this_pp)/3;
ratio(ii) = k_gaps * this_q;
end
pb.end_loop();
end
clear pb
function k = calc_k_gaps( time_serie )
time_serie = time_serie( time_serie >= start_sample & time_serie <= end_sample );
% add some samples in order to calculate the gaps among heartbeats (FN)
if( time_serie(1) ~= start_sample )
time_serie = [start_sample; time_serie];
end
if( time_serie(end) ~= end_sample )
time_serie = [time_serie; end_sample];
end
time_serie = round(time_serie * 1000 / ECG_header.freq);
RRserie = colvec(diff(time_serie));
RRserie = [0;RRserie];
gap_idx = find(RRserie > max_pattern_separation);
% gap_end_idx = gap_idx;
% gap_start_idx = gap_start_idx - 1;
% Percent of the time we have gaps.
k = 1 - (sum(RRserie(gap_idx)) / time_serie(end));
end
function fm = calculate_RR_features( time_serie )
RRserie = diff(time_serie);
RRserie = [RRserie(1); RRserie];
% build the feature matrix [ RR_i-1 RR_i RR_10 RR_60 ]
% RR_i-1
fm = [ RRserie(1); RRserie(1:end-1) ];
fm = [fm RRserie CalcFeatureRRx(time_serie, RRserie, 10000) CalcFeatureRRx(time_serie, RRserie, 60000)];
end
function start_end_aux = findStartEnd( bAux )
start_aux = find( bAux, 1, 'first' );
end_aux = find( bAux, 1, 'last' );
start_end_aux = [start_aux end_aux];
end
function aux_mean = CalcFeatureRRx(time_serie, RRserie, win_size)
% win_size in milliseconds
aux_seq = 1:length(RRserie);
aux_idx = arrayfun(@(a)( findStartEnd( time_serie >= (time_serie(a) - win_size) & ...
time_serie <= time_serie(a) )), ...
aux_seq, 'UniformOutput', false);
aux_mean = colvec(cellfun(@(a)(round(median(RRserie(a(1):a(2))))), aux_idx));
end
end