[b, pp, sf, Lpt2, Lpt1, Baseline]=B_mil_1474D_duration(p, Fs, flag, make_plot)

function   [b, pp, sf, Lpt2, Lpt1, Baseline]=B_mil_1474D_duration(p, Fs, flag, make_plot)
% % B_mil_1474D_duration: Calculates the B-duration according to MIL-STD-1474D Requirement 4.
% %
% % Syntax:  
% % 
% % [b, pp, sf, Lpt2, Lpt1]=B_mil_1474D_duration(p, Fs, flag, make_plot);
% %
% % ********************************************************************
% %
% % Description
% %
% % This program calculates the B-duration for impulsive noise analysis
% %
% % Reference: mil standard 1474D pages 43-4712 February 1997
% %
% % ********************************************************************
% %
% % Input variables
% %
% % p=randn(1, 10000)   % Sound pressure time record in Pa for a single channel data array
% %                     % the default value is randn(1, 50000);
% %
% % Fs=100000;  % Sampling rate in Hz.  
% %             % default value is 100000 Hz.
% %
% % flag=1;     % Multiplies the time record p by -1 to correct for the 
% %             % voltage inversion of polarized microphones.  
% %             % 0 does not alter the time record, p. 
% %             % default value is flag=0;
% % 
% % make_plot=1;    % 1 Plots the processing data of the peaks.  
% %                 % Otherwise no plots are made. 
% %                 % 
% %                 % default is make_plot=0;  
% %
% %
% % ********************************************************************
% %
% % Output variables
% %
% % b is the estimated B-duration in seconds. 
% % 
% % pp is the primary portion of the durationin seconds.
% % 
% % sf is the duration secondary fluctuationsin seconds.
% % 
% % Lpt2 is the index of the first decent below Lplus after the global 
% %      positive peak.
% % 
% % Lpt1 is the index of the first rise above Lplus before the global 
% %      positive peak.
% % 
% % 
% % ********************************************************************
% 
% 
% 
% Example='1';
%
% % Example impulsive data with high background noise
% % The Impulse does not exceed the background noise enough to limit the
% % duration
% 
% Fs=50000; fc=200; td=1; tau=0.1; delay=0.1; A1=3; A2=20;
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
%
% flag=0;           % do not invert the p time record
% make_plot=1;      % plot the simulation
%
% [b, pp, sf]=B_mil_1474D_duration(p, Fs, flag, make_plot);
%
% 
% 
% Example='2';
%
% % Example impulsive data with low background noise
% % The Impulse does not exceed the background noise enough to limit the
% % duration
% 
% Fs=50000; fc=200; td=1; tau=0.1; delay=0.1; A1=0.1; A2=20;
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
%
% flag=0;           % do not invert the p time record
% make_plot=1;      % plot the simulation
%
% [b, pp, sf]=B_mil_1474D_duration(p, Fs, flag, make_plot);
%
% 
% Example='3';
%
% % Example impulsive data with low background noise
% % The impulse is shorter with subsequent flucuating random noise.
% 
% Fs=50000; fc=200; td=0.3; tau=0.01; delay=0.1; A1=0.1; A2=20;
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
%
% flag=0;           % do not invert the p time record
% make_plot=1;      % plot the simulation
% p=[p randn(1, 100) 0.5*randn(1, 1000) 2*randn(1, 100) 0.5*randn(1, 1000) randn(1, 100) 0.1*randn(1, 20000)];
% [b, pp, sf]=B_mil_1474D_duration(p, Fs, flag, make_plot);
%
%
% 
% Example='4';
%
% % Example impulsive data with low background noise
% % The impulse is shorter with subsequent flucuating random noise.
% 
% Fs=50000; fc=2000; td=0.3; tau=0.01; delay=0.1; A1=0.1; A2=20;
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
%
% flag=0;           % do not invert the p time record
% make_plot=1;      % plot the simulation
% p=[p randn(1, 10000) 0.5*randn(1, 1000) 2*randn(1, 100) 0.5*randn(1, 1000) randn(1, 100) 0.1*randn(1, 20000)];
% p(6500)=4;
% [b, pp, sf]=B_mil_1474D_duration(p, Fs, flag, make_plot);
% 
% 
% 
% Example='5';
%
% % Example single impulse data with low background noise
% % The impulse has a very long a-duration with no 
% % subsequent flucuating random noise.
% 
% Fs=50000; fc=200; td=0.3; tau=0.001; delay=0.1; A1=0.1; A2=20;
% [p, t]=analytic_impulse(Fs, fc, td, tau, delay, A1, A2);
%
% flag=0;           % do not invert the p time record
% make_plot=1;      % plot the simulation
%
% [b, pp, sf]=B_mil_1474D_duration(p, Fs, flag, make_plot);
% 
% 
% 
% % ********************************************************************
% % 
% % References: 
% %            
% %            Mil Standard 1474D, DEPARTMENT OF DEFENSE
% %            DESIGN CRITERIA STANDARDm, NOISE LIMITS, 12 February 1997
% %            www.silencertests.com/docs/mil-std-1474d.pdf 
% %            
% %            Guido F. Smoorenburg, "Damage Risk Criteria for Impulsive
% %            Noise," New Perspectives on Noise Induced Hearing Loss, 
% %            Raven Press, New York, pages(471-490) 1982
% % 
% % 
% % ********************************************************************
% %
% % Sub Programs
% % 
% % 
% % List of Dependent Subprograms for 
% % B_mil_1474D_duration
% % 
% % 
% % Program Name   Author   FEX ID#
% % 1) A_duration		
% % 2) B_Duration_second_flucts		
% % 3) convert_double		
% % 4) sub_mean	
% %
% %
% % ********************************************************************
% %
% %
% % B_mil_1474D_duration was written by Edward L. Zechmann
% %
% %     date 22 January     2009    Modified from code from William Murphy
% %
% % modified 31 January     2009    Added plotting of B-duration
% %
% % modified  2 February    2009    Debugged the A_duration
% %
% % modified  8 October     2010    Removed sub_mean2. sub_mean2 was 
% %                                 subtracting the running average.
% %                                 
% %
% % ********************************************************************
% %
% % Please feel free to modify this code.
% %
% % See also: A_Duration, B_Duration, B_mil_1474D_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

if (nargin < 4 || isempty(make_plot)) || ~isnumeric(make_plot)
    make_plot=0;
end

% Make sure that p is a positive pressure corrected for polarization
% effects
if ~isequal(flag, 0)
    p=-p;
end

[p]=convert_double(p);

[m1 n1]=size(p);

if m1 > n1
    p=p';
    [m1 n1]=size(p);
end

% make sure p has only one channel
p=p(1, :);


% Find the positive global peak value
[pmax pos_gl_pk_ix] = max(p);

% Find the negative global peak value
[pmin neg_ix] = min(p);


% Calculate the A-duration
[A, at2, at1]=A_duration(p, Fs, 0);

JA=floor(Fs*A);
Jat2=ceil(Fs*at2);
Jat1=floor(Fs*at1);





% ***********************************************************************
%
% Find the Baseline
%

% Bplus is 5 percent of the difference between the positive and negative
% global peak levels.
Bplus = 0.05*abs(pmax-pmin);

% Calculate number of points in 5 milliseconds of time record
num_pts_5ms = floor(0.005*Fs);

% Calculate number of points in 2 milliseconds of time record
num_pts_2ms = floor(0.002*Fs);

% Caculate the maximum number of steps marching backwards with a 5 ms
% time interval using a 2 ms step
max_iterations=floor((pos_gl_pk_ix-num_pts_5ms)/num_pts_2ms);



if max_iterations < 1
    % not enough pretrigger data to calculate a baseline
    Baseline=0;
else

    % initilize control variables
    e1=0;
    flag3=1;

    % March backward in time using 2 millisecond steps
    % until the baseline is found or the program runs
    % out of data points.
    %
    % This ensures that the 5 ms baseline interval does not precede the
    % first excursion above L+ or below L- by more than 2 ms.
    while (e1 < max_iterations) && isequal(flag3, 1)

        e1=e1+1;
        % ll is the lower limit
        ll=pos_gl_pk_ix-num_pts_5ms-(e1-1)*num_pts_2ms;

        % ul is the upper limit
        ul=pos_gl_pk_ix-(e1-1)*num_pts_2ms;

        if ll < 1
            ll=1;
            flag3=0;
        end

        % p5ms is the time record of the 5 ms baseline interval
        p5ms=p(1, ll:ul);

        % Find the positive local peak value
        [pmax_lo pos_lo_pk_ix] = max(p5ms);

        % Find the negative local peak value
        [pmin_lo neg_lo_pk_ix] = min(p5ms);

        Baseline=mean(p5ms);

        if ((abs(pmax_lo-Baseline) < Bplus) &&  (logical(abs(pmin_lo+Baseline) < Bplus))) || isequal(flag3, 0)
            flag3=0;
            Baseline=mean(p5ms);
        end

    end
end


% Lplus is 20 dB below the positive global peak level
Lplus = pmax*0.1;

% Lplus is 20 dB below the positive global peak level
Lminus = -Lplus + Baseline;



% Find the point before first excursion above L+ before the 
% positive global peak pressure

e1=find(p(1:pos_gl_pk_ix) > Lplus, 1)-1;

if e1 < 1
    e1=1;
end

% Interpolate to find better Lpt1 index.
%
% Lpt1 is the index of the last excursion above L+ before the positive
% global peak pressure.
if abs(p(e1)-p(e1+1)) >= 10^-12
    Lpt1 = (Lplus-p(e1))/(p(e1+1)-p(e1))+e1;
else
    Lpt1=e1;
end


% Find the first point after positive global peak pressure below L+ 
% 

bigend=max([Jat2, length(p)-pos_gl_pk_ix]);

if bigend > length(p)
    bigend=length(p);
end

if pos_gl_pk_ix > length(p)
    pos_gl_pk_ix=length(p);
    bigend=length(p);
end

e1=pos_gl_pk_ix-1+find(p(pos_gl_pk_ix:bigend) < Lplus, 1);

if e1 > n1-1
    e1=n1-1;
end

Lpt22=e1;

% Interpolate to find better Lpt2 index.
%
% Lpt2 is the index of the first excursion below L+ after the positive
% global peak pressure.
if abs(p(e1)-p(e1+1)) >= 10^-12
    Lpt2 = -(Lplus-p(e1))/(p(e1-1)-p(e1))+e1;
else
    Lpt2=e1;
end



% Calculate T1
N=5;
T1=N*A;

% Make sure that T1 is less than 30 ms
if T1 > 0.030
    T1=0.030;
end

% Update the number of data points
t1=floor(Fs*T1);

% **********************************************************************
% 
% Locate Point P
%

% Select all data points after the first crosing of Lplus after the 
% global maximum peak 

buf=p( (Lpt22+1):end);

% Find the points outside the inteval Lplus to Lminus
buf2=find( buf > Lplus);
buf22=find(buf < Lminus);
buf2=union(buf2, buf22);

% Determine if there are any intervals of duration T1 which remain with the 
% range Lplus to Lminus.

% Initialize the flag and Point P and the primary portion
flag4=0;
pp=0;

if length(buf2) > 1
    
    % Determine the lengths of the intervals between data points that 
    % are outside the range Lplus to Lminus
    buf3=diff(buf2);

    % Find the intervals greater than length t1.
    ints_gt_t1=find(buf3 > t1, 1);

    % Determine the first point of the interval.
    if ~isempty(ints_gt_t1)
        P=Lpt22+buf2(ints_gt_t1(1));
    else
        % no intervals can be found
        % Return the last point above Lplus or below Lminus as the B-duration
        
        P=Lpt22+buf2(end);
        pp=1/Fs*(P-Lpt1);
        flag4=1;
    end
    
    while (P-Lpt22 < t1) && ~isempty(ints_gt_t1) && logical(ints_gt_t1(1) < t1) && logical(N > 1)

        % Calculate T1
        N=N-1;
        T1=N*A;

        % Make sure that T1 is less than 30 ms
        if T1 > 0.030
            T1=0.030;
        end

        % Update the number of data points
        t1=floor(Fs*T1);
        
        % Find the first interval larger than t1
        ints_gt_t1=find(buf3 > t1, 1);
        
        % Determine the first point of the interval.
        if ~isempty(ints_gt_t1)
            P=Lpt22+buf2(ints_gt_t1);
        end
        
    end
    
else
    % no intervals can be found
    % Return the end of array p as the B-duration
    P=Lpt22;
    pp=1/Fs*(P-Lpt1);
    
end


flag5=0;
flag6=0;
Q=P;
sf=0;


if isequal(flag4, 0)

    % Quantity, T is defined as
    T=1/Fs*(P-Lpt1);

    R=0.3*T;
    tr1=Fs*R;

    % The primary portion, pp is defined as
    if R <= T1
        pp=T;
    else
        
        % Select all data points after the point P.   
        buf=p( P:end);
        if ~isempty(buf)
            
            % Find the points outside the inteval Lplus to Lminus
            buf2=find( buf > Lplus);
            buf22=find(buf < Lminus );
            buf2=union(buf2, buf22);
            
            
            % Determine the lengths of the intervals between data points 
            % that are outside the range Lplus to Lminus.  
            buf3=diff(buf2);
            
            % Find the intervals greater than length tr1.
            ints_gt_tr1=find(buf3 > tr1, 1);
    
            % Find point Q
            if buf2(ints_gt_tr1) > t1
                Q=P+buf2(ints_gt_tr1);
                flag5=1;
            else
                Q=P;
            end
            
        
        else
            Q=n1;
        end
        
        pp=1/Fs*(Q-Lpt1);

    end


    if Q < n1 && logical(length(p) > 10)
        
        flag6=1;
        % Find durations of the secondary fluctuations
        make_plot2=0;
        [sf, sf2, ctime, spa, mpa, error_cond]=B_Duration_second_flucts(p, Fs, Q, Lplus, Lminus, make_plot2);
        
    end
    

end


if P > n1
    P=n1;
end

if P < 1
    P=1;
end

if sf > 0.1*pp
    b=pp+sf;
else
    b=pp;
end

if isempty(b)
    b=-1;
end


if isequal(make_plot, 1)

    Lt=1;
    
    t=1/Fs*(0:(n1-1));
    figure;
    plot(t, p);
    uno=ones(1,2);
    t2=[t(1) t(end)];
    
    
    % Plot Lplus 
    hold on;
    plot(t2, Lplus*uno,    'k', 'LineWidth', Lt);
    text( 0.9*t(end), 1.1*Lplus,'L_{plus}', 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Bottom');
    
    % Plot Lminus
    plot(t2, Lminus*uno,   'k', 'LineWidth', Lt);
    text( 0.9*t(end), 1.1*Lminus,'L_{minus}', 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Top');
    
    % Plot Baseline
    plot(t2, Baseline*uno, 'r', 'LineWidth', Lt);
    text( 0.9*t(end), Baseline,'Baseline', 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Middle', 'BackgroundColor', [1 1 1]);
    
    % Plot Global Max Peak
    plot( (pos_gl_pk_ix-1)/Fs, p(pos_gl_pk_ix), '^k', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
    text( (pos_gl_pk_ix-1)/Fs, p(pos_gl_pk_ix),'  Global Maximum', 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Bottom');
    
    % Plot Global Min Peak
    plot( (neg_ix-1)/Fs, p(neg_ix), 'vk', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
    text( (neg_ix-1)/Fs, p(neg_ix),'  Global Minimum', 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Bottom');

    % Plot point Lpt1, (Last Point Before excursion above Lplus    
    plot( (Lpt1-1)/Fs, Lplus,  'xk', 'linestyle', 'none', 'markersize', 11, 'LineWidth', Lt);
    text( (Lpt1-1)/Fs-0.001, 1.1*Lplus, 'Point S', 'HorizontalAlignment', 'Right', 'VerticalAlignment', 'Bottom');
    
    % Plot point P
    plot( (P-1)/Fs, p(P), 'ok', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
    if p(P) > 0
        vertical_align='Bottom';
    else
        vertical_align='Top';
    end
    text((P-1)/Fs, 1.1*p(P),'  Point P', 'HorizontalAlignment', 'Left', 'VerticalAlignment', vertical_align);
    
    % Plot point Q
    if isequal(flag5, 1)
        plot( (Q-2)/Fs, p(Q-1), 'sc', 'linestyle', 'none', 'markersize', 7, 'LineWidth', Lt);
        text( (Q-2)/Fs, p(Q-1),'  Point Q', 'HorizontalAlignment', 'Left', 'VerticalAlignment', 'Bottom');
    end
    
    % Plot subsequent fluctuations
    if isequal(flag6, 1)
        
        % put a green circle for each point above the threshold
        for e1=1:length(ctime); 
            plot( (Q-1+ctime(e1))/Fs, p(Q+ctime(e1)), 'og', 'linestyle', 'none', 'markersize', 7);
        end
        
        % put a red circle for each single point above the threshold
        for e1=1:length(spa); 
            plot( (Q-1+spa(e1))/Fs, p(Q+spa(e1)), 'or', 'linestyle', 'none', 'markersize', 7);
        end
        
        % put a black circle for the beginning and end point of each series of multiple points above the threshold
        for e1=1:length(mpa); 
            plot( (Q-1+mpa(e1))/Fs, p(Q+mpa(e1)), 'ok', 'linestyle', 'none', 'markersize', 7);
        end
    end
    
    
    title({'Calculation of B-duration'}, 'Fontsize', 24);
    ylabel('Amplitude (Pa)', 'Fontsize', 18);
    xlabel('Time seconds', 'Fontsize', 18);
    legend_str=cell(8,1);
    legend_str(1:8, 1)={'Signal', 'L_{plus}', 'L_{minus}', 'Baseline', 'Global Max Peak', 'Global Min Peak', 'Point S', 'Point P'};
    if isequal(flag5, 1)
        legend_str{9, 1}='Point Q';
    end
    
    legend(legend_str, 'Location', 'NorthEast'); 
    
    
    rr=(pmax-pmin);
    ylim([0.1*round(10*(pmin-0.1*rr)) 0.1*round(10*(pmax+0.1*rr))]);
    xlim([t(1) t(end)]);
    set(gca, 'Fontsize', 16);
   
   
end