function [cleanNN, cleantNN, flagged_beats] = RRIntervalPreprocess(rr,time,annotations,HRVparams)
%
% [cleanNN, cleantNN, flagged_beats] = RRIntervalPreprocess(rr, time, annotations, HRVparams)
%
% OVERVIEW: This function preprocesses RR interval data in preparation
% for running HRV analyses. Noise and non-normal beats are
% removed from the RR interval data and replaced with
% interpolated data using a interpolation method of choice.
%
% INPUT: rr - a single row of rr interval data in seconds
% time - the time indices of the rr interval data
% (seconds)
% annotations - (optional) annotations of the RR data at each
% point indicating the quality of the beat
% HRVparams - struct of settings for hrv_toolbox analysis
%
% OUTPUT: cleanNN - normal normal interval data
% cleantNN - time indices of NN data
% flagged_beats - percent of original data removed
% REFERENCE:
% Clifford, G. (2002). "Characterizing Artefact in the Normal
% Human 24-Hour RR Time Series to Aid Identification and Artificial
% Replication of Circadian Variations in Human Beat to Beat Heart
% Rate using a Simple Threshold."
% REPO:
% https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
% ORIGINAL SOURCE AND AUTHORS:
% Script written by Adriana N. Vest
% Dependent scripts written by various authors
% (see functions for details)
%
% 09-013-2017
% Edited by Giulia Da Poian
% -Updated the method for measuring changes in the current RR interval
% from the last N (N=5) intervals and excluding intervals that
% change by more than a certain percentage define by
% HRVparams.preprocess.per_limit
% - Removed SQI, windows containing low quality SQI segments are removed
% in EvalTimeDomainHRVstats and EvalFrequencyDomainHRVstats
%
% COPYRIGHT (C) 2016
% LICENSE:
% This software is offered freely and without warranty under
% the GNU (v3 or later) public license. See license file for
% more information
%
% CINC 2002 - Cite for rationale for % good data for cutoff
% 2004 paper - var goes up with numbers of beats dropped out
% don't put phantom beats in for the lomb
% check input
if isempty(time)
time = cumsum(rr);
end
if isempty(annotations)
annotations = repmat('N',[length(rr) 1]);
end
if nargin < 4
error('Not enough input arguments!')
end
figures = HRVparams.gen_figs;
% 1. Identify data that is too close together and Remove
% These are not counted towards the total signal removed
idx_remove = find(diff(time) < 1/HRVparams.Fs);
% could make this the refractory period - and have the variable in the settings
% document
rr(idx_remove+1) = [];
time(idx_remove) = [];
annotations(idx_remove) = [];
clear idx_remove;
% 2. Find Artifact Annotations and Remove Those Points
% These are not counted towards the total signal removed
idx_remove = strcmp('|', annotations); % find artifacts
rr(idx_remove) = [];
annotations(idx_remove) = [];
time(idx_remove) = [];
idx_remove2 = strcmp('~', annotations); % find artifacts
rr(idx_remove2) = [];
annotations(idx_remove2) = [];
time(idx_remove2) = [];
clear idx_remove idx_remove2
% 3. Remove Large Gaps in Data At Beginning of Dataset
%if t(1) > settings.preprocess.forward_gap
% t = t - t(1);
%end
% 4. Remove Large RR intervals Caused by Gaps
% These are not counted towards the total signal removed
idx_remove = find(rr >= HRVparams.preprocess.gaplimit);
rr(idx_remove) = [];
time(idx_remove) = [];
annotations(idx_remove) = [];
clear idx_remove;
% Detect PVCs and Implement Phantom Beat
% Look for premature ectopic beats
% replace value of ectopic beat with a linear interpolated point
% replace sample number
% now next beat becomes normal
% with more than two Vent ectopic beats, cut off data.
% pvcthresh = .3;
% rrdiff = diff(rr);
% idx = find(rrdiff > pvcthresh);
% for i = 1:length(idx)
% ppvc = rr(idx(i)+1);
% if rr(idx(i)+1+1) > rr(idx(i)+1)
% if (rr(idx(i)+1+2) - rr(idx(i))) < pvcthresh
% % If all is true, we can assume we have a PVC
% gap = rr(idx(i)+1+2)-rr(idx(i)+1-1);
% even = gap/3;
% rr(i) = NaN;
% rr(i+1) = NaN;
% avgrr = mean(rr(i-10:i-1));
% rr(i) = even + rr(i-1);
% rr(i+2) = even + rr(i);
% end
% end
% end
%hold on; plot(t,rr);
% 5. Find Non-'N' Annotations
[~,~,ann_outliers] = AnnotationConversion(annotations);
if ~(length(annotations) == length(rr))
ann_outliers = zeros(length(rr),1);
end
% 6. Find RR Over Given Percentage Change
perLimit = HRVparams.preprocess.per_limit;
idxRRtoBeRemoved = FindSpikesInRR(rr, perLimit);
% 7. Find Long Running Outliers
% 8. Combine Annotations and Percentage Outliers
% Combination of methods
outliers_combined = idxRRtoBeRemoved(:) + ann_outliers(:);
outliers = logical(outliers_combined);
% remove extra points that may have been introduced at the end of the file
if (length(outliers) > length(ann_outliers))
fprintf('Too many Outliers - Removing last Point\n');
z = length(outliers);
outliers(z) = [];
end
% timestamp = [find(diff([-1; outliers; -1]) ~= 0)]; % where does outliers change?
% runlength = diff(timestamp); % how long is each segment
% % to split them for 0's and 1's
% runlength0 = runlength(1+(outliers(1)==1):2:end);
% runlength1 = runlength(1+(outliers(1)==0):2:end);
%
% r0_rmv_small = (find(runlength0 > 10));
% r1_rmv_small = (find(runlength1 > 10));
% idx_groups_outliers = find(runlength1 > 10);
% 9. Remove or Interpolate Outliers
idx_outliers = find(outliers == 1);
% Keep count of outliers
numOutliers = length(idx_outliers);
rr_original = rr;
rr(idx_outliers) = NaN;
switch HRVparams.preprocess.method_outliers
case 'cub'
NN_Outliers = interp1(time,rr,time,'spline','extrap');
t_Outliers = time;
case 'pchip'
NN_Outliers = interp1(time,rr,time,'pchip');
t_Outliers = time;
case 'lin'
NN_Outliers = interp1(time,rr,time,'linear','extrap');
t_Outliers = time;
case 'rem'
NN_Outliers = rr;
NN_Outliers(idx_outliers) = [];
t_Outliers = time;
t_Outliers(idx_outliers) = [];
otherwise % By default remove outliers
NN_Outliers = rr;
NN_Outliers(idx_outliers) = [];
t_Outliers = time;
t_Outliers(idx_outliers) = [];
end
if figures
figure;
plot(time,rr_original,t_Outliers,NN_Outliers);
legend('raw','interp1(after outliers removed)')
end
% 10. Identify Non-physiologic Beats
toohigh = NN_Outliers > HRVparams.preprocess.upperphysiolim; % equivalent to RR = 2
toolow = NN_Outliers < HRVparams.preprocess.lowerphysiolim; % equivalent to RR = .375
idx_toolow = find(toolow == 1);
NN_NonPhysBeats = NN_Outliers;
NN_NonPhysBeats(idx_toolow) = NaN;
numOutliers = numOutliers + length(idx_toolow);
switch HRVparams.preprocess.method_unphysio
case 'cub'
NN_NonPhysBeats = interp1(t_Outliers,NN_NonPhysBeats,t_Outliers,'spline','extrap');
t_NonPhysBeats = t_Outliers;
flagged_beats = logical(outliers(:) + toohigh(:)+ toolow(:));
case 'pchip'
NN_NonPhysBeats = interp1(t_Outliers,NN_NonPhysBeats,t_Outliers,'pchip');
t_NonPhysBeats = t_Outliers;
flagged_beats = logical(outliers(:) + toohigh(:)+ toolow(:));
case 'lin'
NN_NonPhysBeats = interp1(t_Outliers,NN_NonPhysBeats,t_Outliers,'linear','extrap');
t_NonPhysBeats = t_Outliers;
flagged_beats = logical(outliers(:) + toohigh(:)+ toolow(:));
case 'rem'
NN_NonPhysBeats(idx_toolow) = [];
t_NonPhysBeats = t_Outliers;
t_NonPhysBeats(idx_toolow) = []; % Review this line of code for improvement
otherwise % use cubic spline interpoletion as default
NN_NonPhysBeats = interp1(t_Outliers,NN_NonPhysBeats,t_Outliers,'pchip');
t_NonPhysBeats = t_Outliers;
end
if figures
hold on;
plot(t_NonPhysBeats,NN_NonPhysBeats+.01);
hold on; plot(t_NonPhysBeats,toolow,'o')
legend('raw','interp1(after outliers removed)',...
'interp2(after too low)','toolow')
end
% 11. Interpolate Through Beats that are Too Fast
toohigh = NN_NonPhysBeats > HRVparams.preprocess.upperphysiolim; % equivalent to RR = 2
idx_outliers_2ndPass = find(logical(toohigh(:)) ~= 0);
NN_TooFastBeats = NN_NonPhysBeats;
NN_TooFastBeats(idx_outliers_2ndPass) = NaN;
numOutliers = numOutliers + length(idx_outliers_2ndPass);
if strcmp(HRVparams.preprocess.method_unphysio,'rem')
flagged_beats = numOutliers;
end
switch HRVparams.preprocess.method_outliers
case 'cub'
NN_TooFastBeats = interp1(t_NonPhysBeats,NN_TooFastBeats,t_NonPhysBeats,'spline','extrap');
t_TooFasyBeats = t_NonPhysBeats;
case 'pchip'
NN_TooFastBeats = interp1(t_NonPhysBeats,NN_TooFastBeats,t_NonPhysBeats,'pchip');
t_TooFasyBeats = t_NonPhysBeats;
case 'lin'
NN_TooFastBeats = interp1(t_NonPhysBeats,NN_TooFastBeats,t_NonPhysBeats,'linear','extrap');
t_TooFasyBeats = t_NonPhysBeats;
case 'rem'
NN_TooFastBeats(idx_outliers_2ndPass) = [];
t_TooFasyBeats = t_NonPhysBeats;
t_TooFasyBeats(idx_outliers_2ndPass) = []; % Review this line of code for improvement
otherwise % USe cubic spline interpoletion as default
NN_TooFastBeats = interp1(t_NonPhysBeats,NN_TooFastBeats,t_NonPhysBeats,'spline','extrap');
t_TooFasyBeats = t_NonPhysBeats;
end
if figures
hold on;
plot(t_TooFasyBeats,NN_TooFastBeats);
legend('raw','interp1(after outliers removed)',...
'interp2(after too low)','toolow','interp3 (after too fast removed)')
end
% 12. Remove erroneous data at the end of a record
% (i.e. a un-physiologic point caused by removing data at the end of
% a record)
while NN_TooFastBeats(end) > HRVparams.preprocess.upperphysiolim % equivalent to RR = 2
NN_TooFastBeats(end) = [];
t_TooFasyBeats(end) = [];
end
cleanNN = NN_TooFastBeats;
cleantNN = t_TooFasyBeats;
end