function r = abpfeature(abp,OnsetTimes)
% ABPFEATURE  ABP waveform feature extractor.
%   r = ABPFEATURE(ABP,ONSETTIMES) extracts features from ABP waveform such
%   as systolic pressure, mean pressure, etc.
%
%   In:   ABP (125Hz sampled), 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
% 
%   Usage:
%   - OnsetTimes must be obtained using wabp.m
%
%   Written by James Sun (xinsun@mit.edu) on Nov 19, 2005.


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

%% P_sys, P_dias
Window    = 40;
OT        = OnsetTimes(1:end-1);
BeatQty   = length(OT);

[MinDomain,MaxDomain] = init(zeros(BeatQty,Window));
for i=1:Window
    MinDomain(:,i) = OT-i+1;
    MaxDomain(:,i) = OT+i-1;
end
MinDomain(MinDomain<1) = 1;  % Error protection
MaxDomain(MaxDomain<1) = 1;

[P_dias Dindex]  = min(abp(MinDomain),[],2);
[P_sys  Sindex]  = max(abp(MaxDomain),[],2);

DiasTime         = MinDomain(sub2ind(size(MinDomain),(1:BeatQty)',Dindex));
SysTime          = MaxDomain(sub2ind(size(MaxDomain),(1:BeatQty)',Sindex));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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));
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/125;  % RR time in seconds
sys_duration        = 0.3*sqrt(RR);
EndOfSys1           = round(OT + sys_duration*125);
SysArea1            = localfun_area(abp,OT,EndOfSys1,P_dias);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Systolic Area calculation using 'first minimum slope' method
SlopeWindow         = 35;
ST                  = EndOfSys1; % start 12 samples after P_sys

if ST(end) > (length(abp)-35)   % error protection
    ST(end) = length(abp)-35;
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; % Get rid of positive slopes

[MinSlope index]   = min(abs(Slope),[],2);
EndOfSys2          = SlopeDomain(sub2ind(size(SlopeDomain),(1:BeatQty)',index));
SysArea2           = localfun_area(abp,OT,EndOfSys2,P_dias);

% OUTPUT:
r = [SysTime'
     P_sys'
     DiasTime'
     P_dias'
     PP'
     MAP'
     BeatPeriod'
     mean_dyneg'
     EndOfSys1'
     SysArea1'
     EndOfSys2' 
     SysArea2']';
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

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