function [SMI,hr] = assess_regularity(QRS,secDer,fs,CI,segL,debug)
% compute smoothness of hrv or of hr time series given a confidence interval (CI).
% The underlying assumption for using this function is that the more smooth
% the QRS time series is the most likely it is to be a meaningful QRS time series.
%
% inputs
% QRS: QRS fiducials (required, in data points number)
% secDer: use second derivative? (i.e. HRV instead of HR to
% compute SMI, default: 1)
% fs: sampling frequency (default: 1kHz)
% CI: confidence interval (default: 0.96)
% segL: length of the ecg segment in seconds (default: 60sec)
%
% outputs
% SMI: regularity measure which corresponds to the standard deviation
% of hrv OR hr (depending on secDer option)
% hr: heart rate
%
%
% FECG extraction toolbox, version 1.0, Sept 2013
% Released under the GNU General Public License
%
% Copyright (C) 2013 Joachim Behar
% Oxford university, 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.
% == manage inputs
if nargin<1; error('assess_regularity: wrong number of input arguments \n'); end;
if nargin<2; secDer = 1; end;
if nargin<3; fs = 1000; end;
if nargin<4; CI = 1; end;
if nargin<5; segL = 60; end;
if nargin<6; debug = 0; end;
% == core function
try
% == compute variability
hr = 60./(diff(QRS)/fs);
if secDer
% rather than looking at the distribution of hr this option is looking at
% the distribution of hrv. This makes more sense because this way
% we are looking into high variation in deltaHR from a measure to the
% following one rather than the variability of absolute value of HR (which
% might be high if the foetus HR is changing significantly)
% use the second derivative of QRS location i.e. ~diff(HR) i.e. hrv
% because we do not know how many QRS have been detected we are going
% to resample at a defined number of datapoints xi
xi = fs*(0.1:0.45:segL-1);
% hr is resampled using a spline function to smooth it. So the output
% is a smoothed version of the hr (yi) at a defined number of
% datapoints (xi)
yi = interp1(QRS(1:end-1),hr,xi,'spline');
% now taking the derivative of smoothed hear rate. This will give the
% hrv
hrv = sort(diff(yi));
hrv_N = length(hrv);
% plot(hr); hold on, plot(yi,'r');
% we tolerate some mistakes using a confidence interval
if CI~=1
CI_sup = ceil(hrv_N*(CI+(1-CI)/2));
CI_inf = ceil(hrv_N*((1-CI)/2));
hrv_CI = hrv(CI_inf:CI_sup);
else
hrv_CI = hrv;
end
% output the std of the hrv
% SMI = std(hrv_CI); % OLD version
% new version (25-08-2013)
SMI = length(find(hrv_CI>30 | hrv_CI<-30));
% this looks at the absolute number of outliers in hrv_CI with >
% 30bpm drop or increase from a point to the next.
else
% use first derivative of RR time series i.e. hr
hr_sorted = sort(hr);
hr_N = length(hr_sorted);
if CI~=1
CI_sup = ceil(hr_N*(CI+(1-CI)/2));
CI_inf = ceil(hr_N*((1-CI)/2));
hr_CI = hr_sorted(CI_inf:CI_sup);
else
hr_CI = hr_sorted;
end
SMI = std(hr_CI);
end
catch ME
for enb=1:length(ME.stack); disp(ME.stack(enb)); end;
SMI = 100; hr = [];
end
% == plots
if debug
hist(hrv_CI,40); xlabel('hr or hrv histogram');
title(['assess regularity in term of hr or hrv, REGULARITY:' SMI]);
set(findall(gcf,'type','text'),'fontSize',14,'fontWeight','bold');
end
end