function [TO, TS,nPVCs, TOsingle, TSsingle , avgTachogram,PVCidxs] = HRT_Analysis(RRInts,Labels,BeatsBefore, BeatsAfter, filtMethod, GraphOn)
% [TO, TS, TOsingle, TSsingle,avgTachogram] = HRT_Analysis(RRInts,Labels,BeatsBefore, BeatsAfter, GraphOn)
% OVERVIEW:
% This function return TO and TS, i.e., the basic parameters of
% heart rate turbulence (HRT) used to quantify the return to
% equilibrium of heart rate after a premature ventricular
% contraction (PVC).
% TO is the Turbulence Onset, and it is defined as the percentage
% difference between average value of the first two normal RR
% following the PVC and the last two normal intervals preciding the
% PVC
% TS is the Turbulence Slope, and it is calculated by constructing an
% average ectopy-centered time series and determining the steepest
% slope each possible five consecutive normal intervals in the
% post-PVC
%
% INPUTS:
% RRInts : Vector containing RR intervals data (in seconds)
% Labels : Vector containing annotations of the RR data at
% each point indicating the type of the beat (see
% https://www.physionet.org/physiobank/annotations.shtml)
% BeatsBefore : number of RR intervals to consider before the PVC
% (default 5)
% BeatsBefore : number of RR intervals to consider after the PVC
% (default 16)
% filtMethod : (optional) select 'mean5before' to filter RR intervals
% that change by more than 20% with respect to the
% mean of the 5 last sinus intervals (the reference interval)
% otherwise it will only filter intervals 15% shorter or
% longer than the preceding interval
% GraphOn : (optional) set GraphOn = 1 for a plot of the HRT analys
%
% OUTPUTS:
% TO : average turbulence onset (TO)
% TS : turbulence slop (TS) of the average tachogram
% nPVCs : number of PVCs used to compute the average tachogram and average TO
% TOsingle : a vector containing the turbulence onset (TO) of each PVC
% TSsingle : a vector containing the turbulence slop (TS) of each PVC
% avgTachogram : average PVC tachogram
% PVCidxs : indexes of the PVCs used in the analysis (relative to the RR position)
%
% DEPENDENCIES & LIBRARIES:
% PhysioNet Cardiovascular Signal Toolbox
% https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%
% REFERENCE:
% Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
% Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%
% REPO:
% https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
% ORIGINAL SOURCE AND AUTHORS:
% Giulia Da Poian
% COPYRIGHT (C) 2018
% LICENSE:
% This software is offered freely and without warranty under
% the GNU (v3 or later) public license. See license file for
% more information
if nargin<2
error('')
end
if nargin<3
BeatsBefore = 5;
end
if nargin<4
BeatsAfter = 16;
end
if nargin<5
filtMethod = 'other';
end
if nargin<6
GraphOn = 0;
end
% Convert to ms
RRInts = RRInts *1000;
% Find PVcs
PVCs = find(Labels=='V');
% Only single ventricular beats are considered
PVCs(diff(PVCs)<=BeatsAfter) = [];
TO = NaN;
TS = NaN;
TOsingle = NaN;
TSsingle = NaN;
avgTachogram = nan(1,2+BeatsBefore+BeatsAfter);
PVCidxs = NaN;
nPVCs = 0;
% Pre-process the original RR intervals to get tachograms time series from
% which unphysiological beats and not Normal 'N' beats are remove, also
% computes the average tachogram
if isempty(PVCs)
return;
end
[tachograms, avgTachogram,PVCidxs, nPVCs] = HRT_preprocessing(RRInts,Labels, PVCs, BeatsBefore, BeatsAfter,filtMethod);
% Analyze tachograms and avgTachogram to get HRT values
[TO, TS, TOsingle, TSsingle] = computeHRT(tachograms, avgTachogram,BeatsBefore);
% Plots
if GraphOn
if isempty(avgTachogram)
return;
end
PlotHRT(tachograms, avgTachogram, BeatsBefore,TO,TS)
end
end % end of Eval_HRT function
function [CleanTachos, AvgTacho,out_idxs_pvcs,nPVCs] = HRT_preprocessing(RRInts,Labels, PVCs, BeatsBefore, BeatsAfter,filtMethod)
% Function that performs a complete HRT preprocessing steps that includes:
% 1. Find all VPC-tachograms (pre_pvc_beats+PVC+CP+post_pvc_beats)
% 2. Filter all VPC-tachograms; i.e. find all the availables VPC-tachograms
% to be used on the HRT analysis
% REFERENCE: Clifford et al. 'Advanced Methods and Tools for ECG Data
% Analysis' (pp. 89-92)
% Initialize variables for:
Tachos = []; % all possible tachograms
TachoLabels = []; % and their labels
CleanTachos = [];
idxPVC = BeatsBefore+1;
idxCP = idxPVC+1;
% Mark all beats after a PVC as C (temporary variable)
tmpLab = Labels;
tmpLab(find(Labels == 'V')+1) = 'C';
tmpLab = tmpLab(1:length(Labels));
% Do not consider RR intervals < 300ms or > 2000ms, unless they are
% PVCs
RRInts( (RRInts<300 | RRInts>2000) & tmpLab~='V' & tmpLab~='C') = NaN;
Labels( (RRInts<300 | RRInts>2000) & tmpLab~='V' & tmpLab~='C') = 'r';
% Mark all RR interval related to 'N' beats that change by more
% than 20% with respect to the median of the five last sinus intervals
% If a PVC tachogrm will have a NaN marked beat it will not be used
idxN = find(tmpLab == 'N'); % use tmpLab to exclude also the CP
if strcmp(filtMethod,'mean5before')
B = 1/5*[1 1 1 1 1];
avgNN = filter(B,1,RRInts(idxN));
refNNinterval = avgNN(5:end-1);
rri_avg = [nan(5,1);RRInts(idxN(6:end))./refNNinterval]; % nan(1,5) for alignment with RRInts indexes
notUse = find(( rri_avg> 1.2) | (rri_avg< 0.8));
else
tmpRR = RRInts(idxN);
rri_rel = [NaN ; tmpRR(2:end)./tmpRR(1:end-1)];
notUse = find((rri_rel < 0.85)| (rri_rel > 1.15) );
end
RRInts(idxN(notUse)) = NaN;
% Guarantee last VPC position allows post_pvc_beats beats after it.
if PVCs(end) +1 + BeatsAfter > length(RRInts) % +1 for the CP after the pvc
PVCs(end) = [];
end
idxs_pvcs = [];
for m = 1: length(PVCs)
idxTachos = PVCs(m)-BeatsBefore : PVCs(m)+1+BeatsAfter;
if idxTachos(end) >= length(RRInts) || idxTachos(1) <= 0 || idxTachos(2) <= 0
continue
end
Tachos = [Tachos RRInts(idxTachos)];
TachoLabels = [TachoLabels Labels(idxTachos)];
idxs_pvcs = [idxs_pvcs PVCs(m)];
end
% Every possible tachogram is going to be "filtered" to determine
% if it is a valid tachogram to be used in HRT analysis
% Indices used to check condition (excluding the PVC and the one after the PVC)
idxSinus = [1:BeatsBefore idxCP+1:idxCP+BeatsAfter];
out_idxs_pvcs = [];
for idx = 1:size(Tachos,2)
tmpTacho = Tachos(:,idx);
tmpLabs = TachoLabels(:,idx);
% compute all the conditions that a VPC tachogram must fulfill
% All beats, except for VPC, must be 'N'
if ~CheckForNormalBeats(tmpLabs, BeatsBefore)
continue;
end
% Do not consider RR intervals where || RR(n-1) - RR(n)|| > 200ms
if ~isempty(find(abs(diff(tmpTacho(1:BeatsBefore)))>=200,1)) || ...
~isempty(find(abs(diff(tmpTacho(BeatsBefore+3:end)))>=200,1))
continue;
end
% Remove tacho with RR intervals that change by more than 20% with
% respect to the mean of the 5 last sinus intervals (the reference interval)
% and also do not consider RR intervals < 300ms or > 2000ms
% Look at RRInts mark as NaN (PVC and CP after PVC can be NaN)
if isnan(sum(tmpTacho(idxSinus)))
continue;
end
% Only use PVCs with a minimum prematurity of 20%
if tmpTacho(BeatsBefore+1) >= 0.8*tmpTacho(BeatsBefore)
continue;
end
% Exclude extrasystolic pause which is at least 20% longer than the normal interval
if tmpTacho(BeatsBefore+2) <= 1.2*tmpTacho(BeatsBefore)
continue;
end
CleanTachos = [CleanTachos tmpTacho];
out_idxs_pvcs = [out_idxs_pvcs idxs_pvcs(idx)];
end
if isempty(CleanTachos)
% there is no tachograms that fulfill all conditions
AvgTacho = [];
nPVCs = 0;
else
AvgTacho = mean(CleanTachos,2);
nPVCs = length(out_idxs_pvcs);
end
end
function validTach = CheckForNormalBeats(tachLab, rri_before)
aux = tachLab([1:rri_before rri_before+2:end]);
cond1 = tachLab(rri_before+1) == 'V';
OtherBeats = find(aux ~= 'N',1);
if isempty(OtherBeats) && cond1
validTach = 1;
else
validTach = 0;
end
end
function [TS, iTS, p]= TurbulenceSlope(Tach,BeatsBefore)
% make sure Tach is an Nx1 vector
if size(Tach,1)<size(Tach,2)
Tach = Tach';
end
%Computes the Turbulence Slope on the tachogram given by parameter
SegLeng = 5;
posPC = BeatsBefore+1;
posFin = length(Tach) - (SegLeng -1) ;
slopes = zeros(1,length(Tach));
idx = 1;
for m = posPC+2:posFin
seg = Tach(m:m+SegLeng-1);
p(idx,:) = polyfit(1:5, seg', 1);
slopes(idx) = p(idx,1);
idx =idx+1;
end
[TS, iTS] = max(slopes);
end
function [TO, num, den] = TurbulenceOnset(Tach,PVCidx)
% Computes the Turbulence Onset on the tachogram given by parameter
num = (Tach(PVCidx+2)+Tach(PVCidx+3))-(Tach(PVCidx-1)+Tach(PVCidx-2));
den = Tach(PVCidx-1)+Tach(PVCidx-2);
TO = (num/den)*100;
end
function [TO, TS, TOsingle, TSsingle ] = computeHRT(CleanTachos, AvgTacho, BeatsBefore)
%Computes the TS and TO on all conditions tachograms (tachograms_ok) and on the mean tachogram
TSsingle = nan(1,size(CleanTachos,2));
TOsingle = nan(1,size(CleanTachos,2));
if isempty(AvgTacho)
% if there is no tachograms that fulfill all conditions, the value of the mean tachogram is nan,
% so the TS and TO values will be nan too
TS = NaN;
TO = NaN;
else
% compute TS and TO for each VPC tachogram
for tac = 1 : size(CleanTachos,2)
TSsingle(tac) = TurbulenceSlope(CleanTachos(:,tac),BeatsBefore);
TOsingle(tac) = TurbulenceOnset(CleanTachos(:,tac),BeatsBefore+1); % BeatsBefore+1 is the PVC position
end
% compute the TS and TO from the average VPC tachogram
TS = TurbulenceSlope(AvgTacho,BeatsBefore);
TO = mean(TOsingle);
end
end
function PlotHRT(tachograms, avgTachogram, pre_pvc_beats,TO,TS)
N = size(tachograms,2);
figure1 = figure('Position',[150 300 1400 600]);
subplot1 = subplot(1,2,1,'Parent',figure1);
plot(tachograms,'Parent',subplot1,'Marker','x');
xlabel('# of RR Interval');ylabel('RR Intervals [ms]');
title({['RR interval profiles of ' num2str(N) ' ventricular premature'], [' complexes (VPC) aligned by VPC position']})
set(subplot1,'FontSize',18);
% Create subplot
subplot2 = subplot(1,2,2,'Parent',figure1);
hold(subplot2,'on');
plot(avgTachogram,'Marker','o','LineWidth',1.5)
xlabel('# of RR Interval');ylabel('RR Intervals [ms]');
set(subplot2,'FontSize',18);
% Plot TO
l1 = (avgTachogram(pre_pvc_beats-1)+ avgTachogram(pre_pvc_beats))/2;
x1 = pre_pvc_beats-2 : pre_pvc_beats+3;
line(x1, l1*ones(length(x1)),'LineWidth',2,'LineStyle','--',...
'Color',[0.929411768913269 0.694117665290833 0.125490203499794]);
l2 = (avgTachogram(pre_pvc_beats+3)+ avgTachogram(pre_pvc_beats+4))/2;
x2 = pre_pvc_beats : pre_pvc_beats+5;
line(x2, l2*ones(length(x2)),'LineWidth',2,'LineStyle','--',...
'Color',[0.929411768913269 0.694117665290833 0.125490203499794]);
% Plot TS
[TS, iTS,p] = TurbulenceSlope(avgTachogram,pre_pvc_beats);
f = polyval(p(iTS,:),1:5);
plot(pre_pvc_beats+1+iTS:pre_pvc_beats+1+iTS+4,f,'LineWidth',3,...
'Color',[0.850980401039124 0.325490206480026 0.0980392172932625])
title('Average Ectopy-Centered Time Series')
% Add TS and TO values
descr = {['Turbulence Onset TO:' num2str(TO) '%'] ;
' ';
['Turbulence Slope TS:' num2str(TS)]};
% Create textbox
annotation(figure1,'textbox',...
[0.74 0.13 0.13 0.2],'Color',[1 1 1],...
'String',descr,'LineWidth',2,'FitBoxToText','off','EdgeColor',[0.08 0.17 0.56],...
'BackgroundColor',[0 0.45 0.75],'FontSize',18);
end