function out = EvalTimeDomainHRVstats(NN,tNN,sqi,HRVparams,windows_all)
% [NNmean,NNmedian,NNmode,NNvariance,NNskew,NNkurt,SDNN,NNiqr,RMSSD, pnn50
% btsdet,avgsqi] = EvalTimeDomainHRVstats(NN,tNN,sqi,HRVparams,windows_all)
%
% OVERVIEW: This function returns time domain HRV metrics calculated on
% input NN intervals.
%
% INPUT: MANDATORY:
% NN : a single row of NN (normal normal) interval
% data in seconds
% tNN : the time indices of the rr interval data
% (seconds)
% sqi : (Optional )Signal Quality Index; Requires
% a matrix with at least two columns. Column
% 1 should be timestamps of each sqi measure,
% and Column 2 should be SQI on a scale from 0 to 1.
% HRVparams : struct of settings for hrv_toolbox analysis
% windows_all : vector containing the starting time of each
% windows (in seconds)
%
% OUTPUT:
% out.NNmean :
% out.NNmedian :
% out.NNmode :
% out.NNvariace :
% out.NNskew : skewness
% out.NNkurt : kurtosis
% out.SDNN : standard deviation
% out.NNiqr :
% out.RMSSD :
% out.pnn50 : the fraction of consecutive beats that differ by
% more than a specified time.
% out.btsdet :
% out.avgsqi :
% out.tdflag : 2 - Not enough high SQI data
% 3 - Not enough data in the window to analyze
% 5 - Success
%
% REPO:
% https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
% ORIGINAL SOURCE AND AUTHORS:
% Written by Adriana N. Vest
% 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
%
% 10-25-2017 Modified by Giulia Da Poian, convert results in ms and
% removed code for writing results on file
%
% Verify input arguments
if nargin< 5
error('no input argument!!!!')
end
if isempty(sqi)
sqi(:,1) = tNN;
sqi(:,2) = ones(length(tNN),1);
end
% Set Defaults
windowlength = HRVparams.windowlength;
alpha = HRVparams.timedomain.alpha;
threshold1 = HRVparams.sqi.LowQualityThreshold;
threshold2 = HRVparams.RejectionThreshold;
% Preallocate arrays (all NaN) before entering the loop
out.NNmean = nan(1,length(windows_all));
out.NNmedian = nan(1,length(windows_all));
out.NNmode = nan(1,length(windows_all));
out.NNvariance = nan(1,length(windows_all));
out.NNskew = nan(1,length(windows_all));
out.NNkurt = nan(1,length(windows_all));
out.NNiqr = nan(1,length(windows_all));
out.SDNN = nan(1,length(windows_all));
out.RMSSD = nan(1,length(windows_all));
out.pnn50 = nan(1,length(windows_all));
out.btsdet = nan(1,length(windows_all));
out.avgsqi = nan(1,length(windows_all));
out.tdflag = nan(1,length(windows_all));
%Analyze by Window
% Loop through each window of RR data
for i_win = 1:length(windows_all)
% Check window for sufficient data
if ~isnan(windows_all(i_win))
% Isolate data in this window
idx_NN_in_win = find(tNN >= windows_all(i_win) & tNN < windows_all(i_win) + windowlength);
idx_sqi_win = find(sqi(:,1) >= windows_all(i_win) & sqi(:,1) < windows_all(i_win) + windowlength);
sqi_win = sqi(idx_sqi_win,:);
t_win = tNN(idx_NN_in_win);
nn_win = NN(idx_NN_in_win);
% Analysis of SQI for the window
lowqual_idx = find(sqi_win(:,2) < threshold1);
% If enough data has an adequate SQI, perform the calculations
if numel(lowqual_idx)/length(sqi_win(:,2)) < threshold2
out.NNmean(i_win) = mean(nn_win) * 1000; % compute and convert to ms
out.NNmedian(i_win) = median(nn_win)* 1000; % compute and convert to ms
out.NNmode(i_win) = mode(nn_win)* 1000; % compute and convert to ms
out.NNvariance(i_win) = var(nn_win)* 1000; % compute and convert to ms
out.NNskew(i_win) = skewness(nn_win);
out.NNkurt(i_win) = kurtosis(nn_win); %
out.NNiqr(i_win) = iqr(nn_win)* 1000; % compute and convert to ms
out.SDNN(i_win) = std(nn_win)* 1000; % compute and convert to ms % SDNN should only be done on longer data segments
% RMSSD
out.RMSSD(i_win) = runrmssd(nn_win)* 1000; % compute and convert to ms
% pNN50
out.pnn50(i_win) = pnna(nn_win, alpha); %
out.btsdet(i_win) = length(nn_win);
out.avgsqi(i_win) = mean(sqi_win(:,2));
out.tdflag(i_win) = 5; % 5 : 'sucess';
else
% 2: low SQI
out.tdflag(i_win) = 2;
end % end of conditional statements that run if SQI is above threshold2
else
% 3: Not enough data in the window to analyze;
out.tdflag(i_win) = 3;
end % end check for sufficient data
end % end of loop through window
end % end of function