function [ac_results, dc_results, prsa_ac, prsa_dc] = prsa(rr, rri, HRVparams, sqi, WinStarIdxs)
% [ac_results, dc_results] = prsa(rr, rri, sqi, windows_all, HRVparams )
%
% OVERVIEW:
% Calculates acceleration and deceleration capacity values
% INPUT
% rr - rr intervals
% rri - index of rr intervals
% HRVparams - struct of settings for hrv_toolbox analysis
% sqi - 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.
% WinStarIdxs - Starting index of each windows to analyze
%
% OUTPUT
% ac_results - Acceleration Capacity value
% dc_results - Deceleration Capacity value
% prsa_ac - PRSA signal for AC anchor points
% prsa_dc - PRSA signal for DC anchor points
%
% AUTHORS
% L. Campana
% G. Clifford
% REFERENCE
% http://adsabs.harvard.edu/abs/2007Chaos..17a5112K
% http://h-r-t.org/prsa/en/
% REPO:
% https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
% DEPENDENCIES & LIBRARIES
% COPYRIGHT (C) 2016 AUTHORS (see above)
% This code (and all others in this repository) are under the GNU General Public License v3
% See LICENSE in this repo for details.
% Feb 13, 2017
% Code altered by Adriana N Vest for the HRV Toolbox. Most of code
% unchanged except for HRVparams and Initializations Section and the
% inclusion of windowing options.
% Oct 13, 2017
% Code altered by Giulia Da Poian to compute AC and DC value at different
% scale s
%
% These default parameters were a part of the original PRSA code developed
% by the authors above.
% tp = 20; thresh_per
% wl = 30; prsa_win_length
% fsd = 512; fs
%%
% Make vector a column
rr = rr(:);
if nargin < 5
error('no data provided')
end
if isempty(sqi)
sqi(:,1) = rri;
sqi(:,2) = ones(length(rri),1);
end
time_stamp = rri;
prsaWinLength = HRVparams.prsa.win_length;
s = HRVparams.prsa.scale;
PRSA_th = HRVparams.prsa.thresh_per;
plot_results = HRVparams.prsa.plot_results;
windowlength = HRVparams.windowlength;
SQI_th = HRVparams.sqi.LowQualityThreshold; % SQI threshold
WinQuality_th = HRVparams.RejectionThreshold; % Low quality windows threshold
Anchor_Low_th = 1-PRSA_th/100-.0001; % The lower limit for the AC anchor selection
Anchor_High_th = 1+PRSA_th/100; % The upper limit for the DC anchor selection
% Preallocation
ac_results = nan(length(WinStarIdxs),1);
dc_results = nan(length(WinStarIdxs),1);
prsa_ac = [];
prsa_dc = [];
% Run PRSA by Window
% Loop through each window of RR data
for i_win = 1:length(WinStarIdxs)
% Check window for sufficient data
if isnan(WinStarIdxs(i_win))
end
if ~isnan(WinStarIdxs(i_win))
% Isolate data in this window
% idx_NN_in_win = find(rri >= windows_all(i_win) & rri < 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( sqi(:,1) >= WinStarIdxs(i_win) & sqi(:,1) < WinStarIdxs(i_win) + windowlength,:);
nn_win = rr( rri >= WinStarIdxs(i_win) & rri < WinStarIdxs(i_win) + windowlength );
% Analysis of SQI for the window
lowqual_idx = find(sqi_win(:,2) < SQI_th);
% If enough data has an adequate SQI, perform the calculations
if numel(lowqual_idx)/length(sqi_win(:,2)) < WinQuality_th
% intialize
acm = [];
dcm = [];
drr_per = nn_win(2:end)./nn_win(1:end-1);
ac_anchor = (drr_per > Anchor_Low_th) & (drr_per <= .9999); % defines ac anchors, no changes greater then 5%
dc_anchor = (drr_per > 1) & (drr_per <= Anchor_High_th);
ac_anchor = [0; ac_anchor(:)];
dc_anchor = [0; dc_anchor(:)];
ac_anchor(1:prsaWinLength) = 0; % sets first window to 0
ac_anchor(length(ac_anchor)-prsaWinLength+1:length(ac_anchor)) = 0; % sets last window to 0
dc_anchor(1:prsaWinLength) = 0; % sets first window to 0
dc_anchor(length(dc_anchor)-prsaWinLength+1:length(dc_anchor)) = 0; % sets last window to 0
ac_ind = find(ac_anchor);
dc_ind = find(dc_anchor);
for i = 1:length(dc_ind)
dcm(i,:) = 1000*nn_win(dc_ind(i)-prsaWinLength:dc_ind(i)+prsaWinLength-1); % convert from sec to msec
end
for i = 1:length(ac_ind)
acm(i,:) = 1000*nn_win(ac_ind(i)-prsaWinLength:ac_ind(i)+prsaWinLength-1); % convert from sec to msec
end
prsa_ac = mean(acm);
prsa_dc = mean(dcm);
% Edited by Adriana Vest: Added the following if statements to
% deal with scenarios when ac or dc is not computable for a
% particular window. An example scenario is when the RR
% intervals are flat.
if ~isnan(prsa_dc) % Edited by Giulia Da Poian
dc = (sum(prsa_dc(prsaWinLength+1:prsaWinLength+s)) - sum(prsa_dc(prsaWinLength-(s-1):prsaWinLength))) ./ (2*s);
dc_results(i_win,1) = dc; % assign output of window
end
if ~isnan(prsa_ac) % Edited by Giulia Da Poian
ac = (sum(prsa_ac(prsaWinLength+1:prsaWinLength+s)) - sum(prsa_ac(prsaWinLength-(s-1):prsaWinLength))) ./ (2*s);
ac_results(i_win,1) = ac; % assign output of window
end
% Load custom colors
custom_colors.red = [72 11 11] / 100;
custom_colors.green = [30 69 31] / 100;
% Plot results
if plot_results == 1
figure(1);
plot(time_stamp,nn_win,'k+');
hold on;
figure(2)
subplot(2,1,2)
plot(dcm,'k--');
legend('Deceleration');
hold on;
figure(2)
subplot(2,1,1)
plot(acm,'k--');
legend('Acceleration');
hold on;
figure(3)
plot(time_stamp,nn_win,'k-+');
hold on;
plot(time_stamp(ac_ind),nn_win(ac_ind), 'color', custom_colors.green, 'marker', '+');
plot(time_stamp(dc_ind),nn_win(dc_ind), 'color', custom_colors.red, 'marker', '+');
legend('non-anchor','ac anchor','dc anchor');
title('RR anchors');
figure(4);
subplot(2,1,1)
plot([-2:2*prsaWinLength-2],dcm,'k--')
title('dc matrix')
hold on
plot([-2:2*prsaWinLength-2],prsa_dc,'r');
subplot(2,1,2)
plot([-2:2*prsaWinLength-2],acm,'k--')
hold on
plot([-2:2*prsaWinLength-2],prsa_ac,'r');
title('ac matrix')
end % end of plotting condition
else % else, if SQI is not adequate
end % end of conditional statements run when SQI is adequate
else % else, if window is NaN
end % end of check for sufficient data
end % end of loop through windows
end % end of function