[sf, sf2, ctime, spa, mpa, error_cond]=B_Duration_second_flucts(p, Fs, Q, Lplus, Lminus, make_plot)

function [sf, sf2, ctime, spa, mpa, error_cond]=B_Duration_second_flucts(p, Fs, Q, Lplus, Lminus, make_plot)
% % B_Duration_second_flucts: Calculates the secondary fluctuations for the B-duration of impulsive noise
% %
% % Syntax:
% %
% % [sf, sf2, error_cond]=B_Duration_second_flucts(p, Fs, Q, Lplus, Lminus, make_plot);
% %
% % **********************************************************************
% %
% % Description
% %
% % This program calculates the duration of the subsequent fluctuations
% % according to MIL-STD-1474D Requirement 4.
% %
% % Two methods are used to calculate the secondary fluctuation durations.
% % Mtehod 1 uses interpolation which is generaly more accurate but may
% % have bugs.  Method 2 is the sum of the points above the threchshold
% % which is less accurate but a quick check for the interpolation method.
% %
% % Using two methods provides for an error condition to calcualte whether
% % there is a mistake.
% %
% % **********************************************************************
% %
% % Input variables
% %
% % p=radnn(1, 50000);  % Sound pressure time record in Pa for a single
% %                     % channel data array the
% %                     % default value is randn(1, 50000);
% %
% % Fs=50000;           % Sampling rate in Hz.
% %                     % default value is 100000 Hz.
% %
% % Q=1;                % Data point that the subsequent fluctuations
% %                     % proceed.
% %                     % default value is [buf Q]=max(abs(p));
% %
% % Lplus=1;            %
% %                     % default value is Lplus=0.1*max(abs(p));
% %
% % Lminus=-1;          %
% %                     % default value is Lminus=-0.1*max(abs(p));
% %
% % make_plot=1;        % 1 Plots the time records and places a circle at
% %                     % each chosen peak. Otherwise no plots are made.
% %                     % default value is 0. (No plotting)
% %
% % **********************************************************************
% %
% % Output variables
% %
% % sf is the interpolated subsequent fluctuation duration in seconds
% %
% % sf2 is the rough estimate of the subsequent fluctuation duration by
% %      counting points without interpolation.
% %
% % ctime array of indices of all points above Lplus or below Lminus.
% %
% % spa array of indices of all single consecutive points above Lplus
% %      or below Lminus.
% %
% % mpa array of indices of the first and last points of multiple
% %      consecutive points above Lplus or below Lminus.
% %
% % error_cond is 1 if the absolute value of the difference betweeen
% %      c and c2 is greater than the sum of the multiple points and the
% %      single points above the threshold.
% %
% % **********************************************************************
%
%
% Example='1';
%
% % Example impulsive data with background noise
%
% 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);
% % p               % Pa sound pressure, single channel data array.
%                   % p should have only one impulse.
% Fs=50000;         % Hz sampling rate, 50000 Hz is the default
% make_plot=0;      % 0 for no plot
%                   % 1 for make the plot
%                   % Default is to not make any plots.
%                   % The plot is the absolute value of the data array.
%                   % blue line for the 10 dB threshold.
%                   % green circle for each point above the threshold.
%                   % red circle for each single point above the
%                   % threshold.
%                   % black circle for each first and last point of a
%                   % multiple point series above the threshold.
%
% [sf, sf2, ctime, spa, mpa, error_cond]=B_Duration_second_flucts(p, Fs, Q, Lplus, Lminus, 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
% %
% %
% % **********************************************************************
% %
% % B_Duration_second_flucts is written by Edward L. Zechmann
% %
% %     date 27 January     2009
% %
% % modified 31 January     2009    Added plotting of B-duration
% %
% % **********************************************************************
% %
% % Please feel free to modify this code.
% %
% % See also: B_mil_1474D_duration, A_Duration, B_Duration, C_Duration, D_Duration
% %


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

% If sampling rate is not specified, then
% assume a data acquisition rate of 50 KHz;
if (nargin < 2 || isempty(Fs)) || ~isnumeric(Fs)
    Fs=100000;
end

if (nargin < 3 || isempty(Q)) || ~isnumeric(Q)
    [buf Q]=max(abs(p));
end

Q=round(Q(1));
Q(Q<1)=1;

if (nargin < 4 || isempty(Lplus)) || ~isnumeric(Lplus)
    Lplus=0.1*max(abs(p));
end

if (nargin < 5 || isempty(Lminus)) || ~isnumeric(Lminus)
    Lminus=-0.1*max(abs(p));
end

% The default is to not make any plots
if (nargin < 6 || isempty(make_plot)) || ~isnumeric(make_plot)
    make_plot=0;
end

% new_zero is the midpoint
new_zero=0.5*(Lplus+Lminus);

% new_thres is the half_length
new_thres=0.5*(Lplus-Lminus);

p3=p((Q+1):end);

p2=abs(p3-new_zero);

% B_Duration_second_flucts
ctime=find(p2 > new_thres);

c=0;
ct1=0;
ct2=0;
c2=length(ctime);

dct=diff(ctime);

flag1=0;
mpa=[];
spa=[];

if isempty(ctime) || logical(length(p2) < 10)

    sf=0;
    sf2=0;
    ctime=0;
    spa=0;
    mpa=0;
    error_cond=0;
    t_array=0;
    t_array=0;

else

    if isempty(dct)


        % C-duration for a single very large impulse
        thres=new_thres/max(max(p3));
        [maxp maxp_index]=max(p2);
        maxp_thres=maxp*thres;
        x_thres_index=find(p2 > maxp_thres);

        if ~isempty(x_thres_index)
            % find start time for C-duration, which is the first 10 dB
            % threshold-crossing before the peak pressure
            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
            flag1=0;
            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
                ct1 = (maxp_thres-p2(e1))/(p2(e1+1)-p2(e1))+e1;
            else
                ct1=e1;
            end

            % find end time for C-duration, i.e. the last 10 dB threshold-crossing
            % after peak
            if (x_thres_index(end)-1) > 1 && ((x_thres_index(end)-1) <= length(p2))
                e1=(x_thres_index(end)-1);
            else
                e1=2;
            end

            flag1=0;
            while (flag1 == 0) && (e1 < length(p2))
                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
                ct2 = (maxp_thres-p2(e1-1))/(p2(e1)-p2(e1-1))+e1-1;
            else
                ct2=e1;
            end

        else
            ct2=1;
            ct1=1;

        end

        % C-duration in indices
        c = (ct2-ct1);
        c2=c;
        error_cond=0;

        t_array=(ct2-ct1);

    else
        % mpa multiple points above threshold
        % spa single point above threshold

        % C-duration for a single very large impulse
        thres=new_thres/max(max(p2));
        [maxp maxp_index]=max(p2);
        maxp_thres=maxp*thres;


        mpa=zeros(c2, 1);
        spa=zeros(c2, 1);

        mpc=0; % multiple point counter
        spc=0; % single point counter

        for e1=1:(length(dct));

            if dct(e1) > 1
                if flag1 == 1
                    mpc=mpc+1;
                    mpa(mpc)=ctime(e1);
                else
                    spc=spc+1;
                    spa(spc)=ctime(e1);
                end
                flag1=0;
            else
                if flag1 == 0
                    mpc=mpc+1;
                    mpa(mpc)=ctime(e1);
                end
                flag1=1;
            end
        end

        mpa=mpa(1:mpc);
        spa=spa(1:spc);

        t_array=[];

        % Every pair of multiple points above threshold
        % must have an end point
        % if no last end point, then the last
        % point becomes the last end point
        if mod(length(mpa), 2) > 0
            mpc=mpc+1;
            mpa(mpc)=ctime(end);
        end

        nmpa=floor(length(mpa)/2);
        t_array=zeros(nmpa+spc, 1);

        for e1=1:nmpa;

            % get height of first data point before crossing threshold
            if (mpa(2*e1-1)-1) > 1
                pt1=p2(mpa(2*e1-1)-1);
            else
                pt1=p2(1);
            end

            % get height of data point after crossing threshold
            pt2=p2(mpa(2*e1-1));

            % interpolate for data point before crossing threshold
            if abs(pt2-pt1) >= 10^-12
                t3L=1*(pt2-maxp_thres)/(pt2-pt1)+0;
            else
                t3L=0;
            end


            % get height of data point
            pt22=p2(mpa(2*e1));

            % get height of next data point
            if (mpa(2*e1)+1) < length(p2)
                pt3=p2(mpa(2*e1)+1);
            else
                pt3=p2(end);
            end

            % interpolate for data point after crossing threshold
            if abs(pt22-pt3) >= 10^-12
                t3U=1*(pt22-maxp_thres)/(pt22-pt3)+0;
            else
                t3U=0;
            end

            t_array(e1)=t3U+t3L+mpa(2*e1)-mpa(2*e1-1);

        end

        for e1=1:spc;

            % get height of previous data point
            if (spa(e1)-1) > 1
                pt1=p2(spa(e1)-1);
            else
                pt1=p2(1);
            end

            % get height of data point
            pt2=p2(spa(e1));

            % interpolate for data point before crossing threshold
            if abs(pt2-pt1) >= 10^-12
                t3L=1*(pt2-maxp_thres)/(pt2-pt1)+0;
            else
                t3L=0;
            end

            % get height of next data point
            if (spa(e1)+1) < length(p2)
                pt3=p2(spa(e1)+1);
            else
                pt3=p2(end);
            end


            % interpolate for data point after crossing threshold
            if abs(pt2-pt3) >= 10^-12
                t3U=1*(pt2-maxp_thres)/(pt2-pt3)+0;
            else
                t3U=0;
            end

            t_array(e1+nmpa)=t3U+t3L;

        end

        % calculate the interpolated sum
        c=sum(t_array);

        % calculate the rough estimated sum
        c2=length(ctime);

        % Calculate if the maximum error is satisfied
        error_cond=0;
        max_error=spc + nmpa;
        if abs(c2 - c) > max_error
            error_cond=1;
        end

        if make_plot == 1
            close all;

            figure(2);
            % plot the absolute value of the data array
            plot(1/Fs*(0:length(p)-1), abs(p));
            hold on;

            % draw the threshold line
            plot( 1/Fs*[0 length(p)-1], maxp_thres*[1 1], 'k');

            % put a green circle for each point above the threshold
            for e1=1:length(ctime);
                plot( (Q-1+ctime(e1))/Fs, p2(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, p2(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, p2(mpa(e1)), 'ok', 'linestyle', 'none', 'markersize', 7);
            end

        end

    end

end

% convert c and c2 to seconds
sf=1/Fs*c;
sf2=1/Fs*c2;

if isempty(sf)
    sf=-1;
end

if isempty(sf2)
    sf2=-1;
end