function r = abpfeature(abp,OnsetTimes, Fs)
% ABPFEATURE  ABP waveform feature extractor.
%   r = ABPFEATURE(ABP,ONSETTIMES) extracts features from ABP waveform such
%   as systolic pressure, mean pressure, etc.
%
%   In:     ABP = signal (default 125Hz sampled)
%           ONSETTIMES = times of onset (in samples)
%   Out:    Beat-to-beat ABP features
%           Col 1:  Time of systole   [samples]
%               2:  Systolic BP       [mmHg]
%               3:  Time of diastole  [samples]
%               4:  Diastolic BP      [mmHg]
%               5:  Pulse pressure    [mmHg]
%               6:  Mean pressure     [mmHg]
%               7:  Beat Period       [samples]
%               8:  mean_dyneg
%               9:  End of systole time  0.3*sqrt(RR)  method
%              10:  Area under systole   0.3*sqrt(RR)  method
%              11:  End of systole time  1st min-slope method
%              12:  Area under systole   1st min-slope method
%              13:  Pulse             [samples]
% 
%   Usage:
%   - OnsetTimes must be obtained using wabp.m
%
%   Written by James Sun (xinsun@mit.edu) on Nov 19, 2005.
%   Updated by Alistair Johnson, 2014.
%
%   LICENSE:    
%       This software is offered freely and without warranty under 
%       the GNU (v3 or later) public license. See license file for
%       more information




if length(OnsetTimes)<30 % don't process anything if too little onsets
    r = [];
    return
end

%% P_sys, P_dias
if nargin<3
Fs        = 125; % Sampling Frequency
end
Window    = ceil(0.32*Fs);
OT        = OnsetTimes(1:end-1);
BeatQty   = length(OT);

% [MinDomain,MaxDomain] = init(zeros(BeatQty,Window));

for i=1:Window % Vectorized version
    MinDomain(:,i) = OT-i+1;
	MaxDomain(:,i) = OT+i-1;
end

%MinDomain = bsxfun(@minus,OT,(0:Window-1)); 
%MaxDomain = bsxfun(@plus,OT,(0:Window-1)); 
MinDomain(MinDomain<1) = 1;  % Error protection
MaxDomain(MaxDomain<1) = 1;
MinDomain(MinDomain>length(abp)) = length(abp);
MaxDomain(MaxDomain>length(abp)) = length(abp);

[P_dias Dindex]  = min(abp(MinDomain),[],2);% Get lowest value across 'Window' samples before beat onset
[P_sys  Sindex]  = max(abp(MaxDomain),[],2);% Get highest value across 'Window' samples after beat onset

DiasTime         = MinDomain(sub2ind(size(MinDomain),(1:BeatQty)',Dindex)); % Map offset indices Dindex to original indices
SysTime          = MaxDomain(sub2ind(size(MaxDomain),(1:BeatQty)',Sindex)); % Map offset indices Sindex to original indices

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pulse Pressure [mmHg]
PP         = P_sys - P_dias;

% Beat Period [samples]
BeatPeriod = diff(OnsetTimes);

% Mean,StdDev, Median Deriv- (noise detector)
dyneg          = diff(abp);
dyneg(dyneg>0) = 0;

[MAP,stddev,mean_dyneg] = init(zeros(BeatQty,1));
if OnsetTimes(end)==numel(abp)
    OnsetTimes(end) = numel(abp) - 1;
end

for i=1:BeatQty
    interval       = abp(OnsetTimes(i):OnsetTimes(i+1));
    MAP(i)         = mean(interval);
    stddev(i)      = std(interval);

    dyneg_interval = dyneg(OnsetTimes(i):OnsetTimes(i+1));
    dyneg_interval(dyneg_interval==0) = [];
    if min(size(dyneg_interval))==0
        dyneg_interval = 0;
    end
    mean_dyneg(i)  = mean(dyneg_interval);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Systolic Area calculation using 0.3*sqrt(RR)
RR                  = BeatPeriod/Fs;  % RR time in seconds
sys_duration        = 0.3*sqrt(RR);
EndOfSys1           = round(OT + sys_duration*Fs);
SysArea1            = localfun_area(abp,OT,EndOfSys1',P_dias);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Systolic Area calculation using 'first minimum slope' method
SlopeWindow         = ceil(0.28*Fs);
ST                  = EndOfSys1; % start 12 samples after P_sys

if ST(end) > (length(abp)-SlopeWindow)   % error protection
    ST(end) = length(abp)-SlopeWindow;
end

SlopeDomain = zeros(BeatQty,SlopeWindow);
for i=1:SlopeWindow
    SlopeDomain(:,i) = ST+i-1;
end
% y[n] = x[n]-x[n-1]
Slope              = diff(abp(SlopeDomain),1,2);
Slope(Slope>0)     = 0; % Set positive slopes to zero

[MinSlope index]   = min(abs(Slope),[],2); % Find first positive slope
EndOfSys2          = SlopeDomain(sub2ind(size(SlopeDomain),(1:BeatQty)',index));
SysArea2           = localfun_area(abp,OT,EndOfSys2,P_dias);
Pulse              = 60./BeatPeriod;


% OUTPUT:
% Ensure that there is no concatenation error by using
% (:) indexing to force column vectors
r = [SysTime(:),P_sys(:),DiasTime(:),P_dias(:),PP(:),...
    MAP(:),BeatPeriod(:),mean_dyneg(:),EndOfSys1(:),SysArea1(:),...
    EndOfSys2(:),SysArea2(:),Pulse(:)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function:

function SysArea = localfun_area(abp,onset,EndSys,P_dias)
%% Input:  abp,
%%         onset, P_dias, end of systole, beat duration in unit of samples, same length
%% Output: systolic area, warner correction factor

BeatQty = length(onset);
SysArea = init(zeros(BeatQty,1));
for i=1:BeatQty
    SysArea(i) = sum(abp(onset(i):EndSys(i))); % faster than trapz below
    %b(i) = trapz(abp(onset(i):EndSys(i)));
end
EndSys = EndSys(:); onset = onset(:); % force col vectors
SysPeriod  = EndSys     - onset;

% Time scale and subtract the diastolic area under each systolic interval
SysArea = (SysArea - P_dias.*SysPeriod)/125; % Area [mmHg*sec]
end

% for initializing a variable
function varargout = init(z)
for i=1:nargout
    varargout(i) = {z};
end
end