[b_min, b, p2]=E_duration(p, Fs)

function [b_min, b, p2]=E_duration(p, Fs)
% % E_duration: Calculates the E-duration for impulsive noise
% % 
% % Syntax:  [b_min, b, p2]=E_duration(p, Fs);
% %
% % **********************************************************************
% %
% % Description 
% % 
% % This program calculates the E-duration for impulsive noise analysis
% % 
% % 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
% % 
% % **********************************************************************
% % 
% % 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.  
% % 
% % **********************************************************************
% %
% % Output variables
% %
% % b_min is the smallest estimate of the E-duration 
% % It is recommended to use b_min as the estimated E-duration.
% % 
% % b is the array of estimated E-durations for each threshold value
% % thresholds form 1 dB to 20 dB are used to make the esimate more robust.
% % Some thresholds typically will have anomalies and overestimate the 
% % E-duration.  
% %
% % p2 is the processed data array.
% % The processing involvse the running average of the absolute value of 
% % the hilbert transform.
% % 
% % This quantity, p2, is the time average instantaneous sound pressure.
% %
% % **********************************************************************
% %
% Example;
%
% % 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
% 
% [b_min, b, p2]=E_duration(p, Fs);
% 
% % 
% % **********************************************************************
% % 
% % 
% % Subprograms
% % 
% % 
% % 
% % List of Dependent Subprograms for 
% % E_duration
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % 
% % Program Name   Author   FEX ID#
% % 1) convert_double		Edward L. Zechmann			
% % 2) sub_mean2		Edward L. Zechmann						
% % 
% % 
% % **********************************************************************
% %
% % E_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 22 January     2009    Updated Comments
% % 
% % modified 19 January     2010    Changed name to E-Duration 
% %                                 from B-Duration.  Updated Comments
% % 
% % modified 11 February    2011    Updated Comments
% %                                 
% % 
% % 
% % 
% % **********************************************************************
% %
% % Please feel free to modify this code.
% % 
% % See also: A_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

[buf, p2]=sub_mean2(abs(hilbert(p)), Fs, 2000);

% array of threshold values in dB 
threshold_a=[1, 2, 4, 8, 10, 15, 20]; % dB
num_thres=length(threshold_a);

% E-duration (assuming only one impulse)

[maxp maxp_index]=max(abs(p2));
b=zeros(num_thres, 1);
num_data=length(p);

for e2=1:num_thres;
    thres=10^(-threshold_a(e2)/20);
    maxp_thres=maxp*thres;
    x_thres_index=find(p2 > maxp_thres);

    % find start time for E-duration, which is the first 20 dB threshold-crossing before
    % the peak pressure
    flag1=0;
    e1=x_thres_index(1)+1;

    if (x_thres_index(1)+1) > 1 && (x_thres_index(1)+1) < length(p2)
        e1=x_thres_index(1)+1;
    else
        e1=1;
    end

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

    % interpolate to find better begin time
    if abs(p2(e1)-p2(e1+1)) >= 10^-12
        bt1 = (maxp_thres-p2(e1))/(p2(e1+1)-p2(e1))+e1;
    else
        bt1=e1;
    end

    % find end time for E-duration, i.e. the last 20 dB threshold-crossing after peak
    flag1=0;

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

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

    % interpolate to find better end time
    if abs(p2(e1)-p2(e1-1)) >= 10^-12
        bt2 = (maxp_thres-p2(e1-1))/(p2(e1)-p2(e1-1))+e1-1;
    else
        bt2=e1;
    end

    % E-duration in indices
    b(e2) = (bt2-bt1);

end

% normaize b by the 20 dB threshold  and convert to seconds
b=1/Fs*(20./threshold_a)'.*b;

% Recommended to use the minimum of the estimated E-durations
% 
% Background noise inceases the E-duration 
% Flat or jumpy subsequent peaks can increase the E-duration 
% 
% Both of the above affects excessively increase the E-duration

% Make sure that the E-duration is positive
b_min=min(b(b > 0.00000001));

% Make sure that E-duration is not longer than the data array
if b_min > 1/Fs*num_data
    b_min=1/Fs*num_data;
end

if isempty(b)
    b=-1;
    b_min=-1;
end