[a, at2, at1]=A_duration(p, Fs, flag)

function [a, at2, at1]=A_duration(p, Fs, flag)
% % A_duration: Calculates the A-duration for impulsive noise
% % 
% % Syntax:  
% %
% % [a, at2, at1]=A_duration(p, Fs, flag);
% %
% % **********************************************************************
% %
% % Description 
% % 
% % This program calculates the A-duration for impulsive noise analysis
% % 
% % **********************************************************************
% % 
% % Input variables
% % 
% % p is the sound pressure in Pa for a single channel data array
% %                 the default value is randn(1, 50000);
% %
% % Fs sampling rate in Hz.  default value is 100000 Hz.  
% % 
% % flag=0; % Controls whether to use p, -p or abs(p) for identifying the
% %         % global peak pressure. The A-duration is calculated using the 
% %         % zero crossings before an after the global peak. 
% %         % 
% %         % 0 is for ICP microphones
% %         % 1 is for 200 Volt polarized microphones
% %         % 
% %
% % **********************************************************************
% %
% % Output variables
% %
% % a is the A-duration in seconds
% %
% % at2 is the time of first zero crossing after the peak in seconds
% %
% % at1 is the time of first zero crossing before the peak in seconds
% %
% % **********************************************************************
% 
% 
% 
% Example='1';
%
% % Example impulsive data with background noise
%
% Fs=50000; fc=200; td=1; tau=0.1; delay=0.1; A1=3; A=20;
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
% % p               % Pa sound pressure, single channel data array.
%                   % p should have only one impulse.
% Fs=50000;         % Hz sample rate frequency
% flag=0;
% 
% [a, at2, at1]=A_duration(p, Fs, flag);
% 
% % **********************************************************************
% % 
% % Reference: Guido F. Smoorenburg, "Damage Risk Criteria for Impulsive
% %            Noise," New Perspectives on Noise Induced Hearing Loss, 
% %            Raven Press, New York, pages(471-490) 1982
% % 
% % 
% %
% % **********************************************************************
% %
% % A_duration.m was originally developed by Chucri Kardous.  
% %
% % This implementation of A_duration was written by Edward L. Zechmann  
% % 
% %     date 11 December    2007
% % 
% % modified 17 December    2007    Added Comments 
% % 
% % modified 13 August      2008    Updated Comments
% % 
% % modified 21 September   2008    Check output for being empty
% %                                 Updated Comments
% % 
% % modified  6 August      2010    Added additional error handling
% %                                 Updated Comments
% % 
% % 
% % 
% % 
% % **********************************************************************
% %
% % Please feel free to modify this code.
% % 
% % See also: B_mil_1474D_duration, B_Duration, C_Duration, D_Duration
% % 


if (nargin < 1 || isempty(p)) || ~isnumeric(p)
    p=randn(1, 50000);
end

if (nargin < 2 || isempty(Fs)) || ~isnumeric(Fs)
    Fs=50000;
end

if (nargin < 3 || isempty(flag)) || ~isnumeric(flag)
    flag=0;
end



% A-duration

if isequal(flag, 0)
    p2=p;
elseif isequal(flag, 1)
    p2=-p;
else
    p2=abs(p);
end

[maxp maxp_index]=max(p2);

p_max=p(maxp_index);

if p_max < 0;
    p=-p;
end


% find start time for A-duration, which is the first zero-crossing before
% the peak pressure
flag1=0;
if (maxp_index > 1) && (maxp_index <= length(p))
    e1=(maxp_index);
else
    e1=1;
end

while (flag1 == 0) && (e1 > 1)
    e1 = e1-1;
    if p(e1) <= 0
        flag1=1;
    end
end

% interpolate to find better begin time
if abs(p(e1)-p(e1+1)) >= 10^-12 && ~isequal(flag1, 0)
    at1 = (0-p(e1))/(p(e1+1)-p(e1))+e1;
else
    at1=e1;
end

% find end time for A-duration, i.e. the last zero-crossing after peak
flag1=0;

if (maxp_index > 1) && (maxp_index <= length(p))
    e1=(maxp_index);
else
    e1=2;
end

while (flag1 == 0) && (e1 < length(p))
    e1 = e1+1;
    if p(e1) <= 0
        flag1=1;
    end
end

% interpolate to find better end time
if abs(p(e1)-p(e1-1)) >= 10^-12 && ~isequal(flag1, 0)
    at2 = (0-p(e1-1))/(p(e1)-p(e1-1))+e1-1;
else
    at2=e1;
end

% A-duration in indices
% end time minus beginning time
at2=1/Fs*at2;
at1=1/Fs*at1;

a = (at2-at1);

if isempty(a);
    a=-1;
end