function  [rt_array2, lv_array2, lv_array_BG2]=rt_script(filenames_in, Fc, pbs, pab, outfile, num_x_filter, t_high, t_low, Room_Name, plot_all_mics, smd, num_std, method)
% % Syntax;
% % [rt_array2, lv_array2, lv_array_BG2, ptsa, nptsa]=rt_script(filenames_in, Fc, pbs, pab, outfile, num_x_filter, t_high, t_low, Room_Name, plot_all_mics, smd, num_std, method)
% % 
% % This program calculates the reverberation time in one third octave 
% % bands.  This program supports two measuremnt methods: 1) the speaker 
% % on speaker off method and 2) the balloon pop method.   The first method
% % using loud speakers is considered to be mmore accurate since it excites
% % the room modes more effectively.  The first method generally has a 
% % better signal to noise ratio than the second method.  However somtimes
% % it is difficult to transport and setup the equipment for the speaker on
% % speakre off method.  It is usually much easier anf more fun to pop
% % balloons.  The balloon pop method may be sufficiently accurate 
% % and is usually much more convenient.  
% % 
% % The data is processed by filtering the data into a third octave bands.
% % The instantaneous sound pressure level is calculated using the 
% % hilbert transform.  
% % 
% % A moving time averaged mean value with a 40 ms sample rate for the 
% % speaker on-off method and 10 ms sample rate for the balloon pop method
% % is used to smooth the instantaneous sound pressure level.  
% % 
% % For the speaker on-off method, the signal high level is calculated.   
% % For the balloon pop method the maximum data point is found.
% % 
% % The background level is estimated from the last t_low seconds of data.
% % 
% % Cutoff times from the signal high to the signal low are calculated.  
% % The signal cutoff are after the speaker was shut off or after the 
% % balloon popped.  The background cutoff is before the time record 
% % reaches the background level.  
% % 
% % The instantaneous sound pressure level data between the cutoffs is 
% % the transient reverberant decay.  
% % 
% % A decay line is robustly fit to the transient reverberant decay.
% % 
% % The slope of the transient instantaneous amplitude is used to calculate 
% % the T60 reverberation time.
% % 
% % The program outputs a pdf file for each center frequency and one tab-
% % delimited text file which serve as documentation of the reverberation 
% % time calculations.
% % 
% % The program automatically loads the matlab files from the current
% % working directory.  Each data file must contain a data variable which 
% % is the time record of the sound pressure.  The units must be linear.  
% % The program is setup for units of Pa and a size eqaual to 
% % (num_mics, num_data_points).
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %  
% % See reverb_time.m program for a baloon pop example. 
% % Here is a speaker on-off example and a description of the input values.
% % 
% Example=''; % Simulated reverberation data 
%             % For the Speaker on Speaker off method
% 
% % This example is for a four microphone with each microphone recording at
% % a different level.  The reverberation time is set to a different for each
% % microphone to illustrate the decay line fitting.
%
% filenames_in={'speaker_on_example.mat'};   
%                           % cell array of filenames
%
% Fc=400;                   % Center Frequency array default has 31 bands
%                           % from 20 Hz to 20 KHz.
%                           %
%                           % There is a one-to-one correspondence between
%                           % The array Fc and the filenames
%                           %
%                           % Each file is assigned the next
%                           % center frequency in the array Fc
%                           
% pbs=15;                   % Percentage of the total range in dB below the average signal
%                           % level to discard from the curve fit.  Increase pbs to avoid
%                           % spurious high level sound data at the speaker cutoff time
%                           %
% pab=40;                   % Percentage of the range above the average background level
%                           % to discard from the curve fit. ncrease pab to avoid spurious background
%                           % noise, where the reveberant decay approaches the background level.
%                           % This is more of a problem in environments with very low background levels;
%                           % for example laboratory hemi-anechoic chambers;
%                           
% pab=0;                    % 0 is a good number for acoustic environments with high continuous background levels;
%                           % for example animal shelters with dogs barking in the surrounding rooms.
%                           
% outfile='Reverb_data';    % Filename for the data for each microphone.  The data includes
%                           % the reverberation times, high sound level at the beginning of the recording,
%                           % and the background sound pressure level at the end of the recording.
%                           % If more than 1 microphone is used then
%                           % standard deviation, confidence interval and outlier data are also printed.
%                           
% num_x_filter=2;           % 2 is a typical value.  Number of times to filter the data.
%                           % Minimum value is 1.
%                           % Typically a value of 2 to 10 is good.  At low
%                           % frequencies (Fc < 100), num_x_filter=10 has a significant phase shift.
%                           
% t_high=0.5;               % (s) duration of signal on at the beginning of the time record.
%                           % t_high is used to calculate the average
%                           % signal level.  The speaker cutoff time is
%                           % estimated using the average signal level.
%                           % A longer t_high time is better when there is filter phase delay
%                           % especially at low frequencies (Fc < 100) or if the speaker turns on slowly.
%                           
% t_low=0.5;                % (s) duration of quiet time at the end of the time record
%                           % used to calculate the average background
%                           % level, which is used to estimate the background cutoff time.
%                           % Background cutoff occurs when the reverberant decay approaches the
%                           % background noise level.  A longer t_low time helps to avoid
%                           % the dips and peaks in the background level. Dips and peaks
%                           % can be caused by sudden quiet and background dog barks etc...
%                           
% Room_Name='Auditorium';   % The room name appears on the pdf plot
%                           
% plot_all_mics=1;          % 1 plot all microphone time records
%                           % 0 plots only the microphone time record for the
%                           % microphone with the reverberation time
%                           % closest to the mean value
%                           
% smd=0;                    % acronym for "surpress microphone data" from the pdf plot legend.
%                           % smd=0; Allow all mics to have a level, reverberation time, and
%                           % quality marker on the plot legend.
%                           % smd=1; Allow only 1 mic to have a level, reverberation time, and
%                           % quality marker on the plot legend
%                           %
% num_std=2;                % Keep all reverberation times within 
%                           % 2 standard deviations of the geometric mean
%                           % num_std is the number of standard deviations 
%                           % to accept data.  1 is good for data with high
%                           % background noise, 2 is better for laboratory
%                           % measurements.  
%
% method=1;                 % 1 is the speaker on speaker off method.  
%                           % otherwise is the balloon pop method
%
% % Now create the file
% 
% tau=[2, 2.3, 3.0, 2.4];   % time constants for exponential decay
% 
% tb=5;                     % seconds, duration of decay signal
% 
% Fs=50000;                 % sampling frequency rate
% 
% t=1/Fs*(1:(tb*Fs));       % seconds time array for calculating
%                           % decay signal
% 
% A=[100, 110, 115, 130];   % RMS Amplitude in (Pa)
% decay1=A(1).*exp(-1./tau(1).*t).*randn(1,length(t));
% decay2=A(2).*exp(-1./tau(2).*t).*randn(1,length(t));
% decay3=A(3).*exp(-1./tau(3).*t).*randn(1,length(t));
% decay4=A(4).*exp(-1./tau(4).*t).*randn(1,length(t));
% 
% % The test signal has
% %     1) nominal one second of signal high
% %     2) reverberant decay
% %     3) nominal one second of background noise
% %
% SP=zeros(4, 2*Fs+length(t));
% SP(1, :)=[A(1)*randn(1,Fs)       decay1 A(1).*exp(-1./tau(1).*tb)*randn(1, Fs);];
% SP(2, :)=[A(2)*randn(1,1.5*Fs)   decay2 A(2).*exp(-1./tau(2).*tb)*randn(1, 0.5*Fs);];
% SP(3, :)=[A(3)*randn(1,1.2*Fs)   decay2 A(3).*exp(-1./tau(3).*tb)*randn(1, 0.8*Fs);];
% SP(4, :)=[A(4)*randn(1,0.8*Fs)   decay2 A(4).*exp(-1./tau(4).*tb)*randn(1, 1.2*Fs);]; 
%
% filenames_in={'speaker_on_example.mat'};
% t_SP=1/Fs*(1:length(SP)); 
% save(filenames_in{1}, 'SP', 't_SP');
%
% [rt_array2, lv_array2, lv_array_BG2]=rt_script(filenames_in, Fc, pbs, ...
% pab, outfile, num_x_filter, t_high, t_low, Room_Name, plot_all_mics, smd, num_std, method);
% 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % Output Variable Description
% % 
% % rt_array2       cell array of reverberation times
% % 
% % lv_array2       cell array of A-weighted average signal levels during t_high
% % 
% % lv_array_BG2    cell array of A-weighted average background levels during t_low
% %
% % The size of the ouput cell arrays are {num_files, 1} 
% % each cell has an array size (num_mics, 1)
% %
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% % Uses the Signal Processing Toolbox
% % Uses the following progrmas which can be found on matlab central
% % Uses the octave program set from Christophe Couvreur 
% % Uses programs from Alexandros Leontitsis
% % Uses the modified programs from Alexandros Leontitsis
% % 
% % List of Subprograms
% % 
% % data_loader2
% %     file_extension
% %     convert_double
% % reverb_time
% %     oct3dsgn    From Christophe Couvreur found on matlab central
% %     sub_mean
% %     LMTSreg     A modified program original from Alexandros Leontitsis
% %                     found on matlab central
% %         LMSloc  From Alexandros Leontitsis found on matlab central
% %         LMS_trim
% %             rand_int
% %     LMTSregor   A modified program original from Alexandros Leontitsis
% %                     found on matlab central
% % adsgn           From Christophe Couvreur found on matlab central
% % cdsgn           From Christophe Couvreur found on matlab central
% % Leq_all_calc
% %     Aweight_time_filter
% %     Cweight_time_filter
% % call_ci
% % ci              from Raymond Reynolds found on matlab central
% % 
% % save_a_plot_reverb_time
% % file_extension
% %
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% % Program was written by Edward L. Zechmann
% % 
% %      date  16 October   2007
% %
% %  modified   9 November  2007    revision 1 revised comments
% %                                 (grammar and clarity)
% %
% %  modified   2 January   2008    updated comments
% %
% %  modified  13 January   2008    Added the data_loader2
% %
% %  modified  10 February  2008    Revised the balloon pop method
% %
% %  modified  13 February  2008    Further revised the balloon pop method
% %                                 Added an example.  Updated comments    
% %
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Please feel free to modify this code as you please.
% %

figure(1);
delete(1);

% program assumes all matlab files in directory are recordings

num_files=length(filenames_in);

if isempty(Fc)
    Fc=[20 25 31 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000 20000 ];
end

if length(Fc) < num_files
    fprintf(1, '%s\r', 'Not enough Center Frequencies for input Fc');
    fprintf(1, '%s\r', 'Adding more Center Frequencies for input Fc');
    num1=ceil( (num_files-length(Fc))/31);
    for e1=1:num1;
        Fc=[Fc [20 25 31 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000 20000] ];
    end
    fprintf(1, '%s\r', 'New List of Center frequencies');
    fprintf(1, '%d\t', Fc);
    fprintf(1, '\r');
end

% The program processes the first file with the center frequency of Fc(1)


count=0;
rt_array2=cell(num_files, 1);    % Reverberation time cell array
lv_array2=cell(num_files, 1);    % Reverberation Level cell array
lv_array_BG2=cell(num_files, 1); % Background Level cell array
meanlevel=0;
mic_num_array=zeros(1,num_files);
ptsa=cell(1,num_files);
nptsa=cell(1,num_files);

default_wav_settings_in={};
default_mat_config_in={};

if num_x_filter < 1
    num_x_filter=1;
end

for e1=1:num_files;

    %reset filter coefficients 
    % B and A are for A-weighted time record
    Ba=1;
    Aa=1;
    % Bc and Ac are for C-weighted time record
    Bc=1;
    Ac=1;
    
    % load the data file
    [SP_var, vibs_var, Fs_SP_var, Fs_vibs_var, default_wav_settings_out, default_mat_config_out]=data_loader2(filenames_in{e1}, default_wav_settings_in, default_mat_config_in);
    
    % set the default wav file setttings
    % for configuring the channels
    default_wav_settings_in=default_wav_settings_out;
    default_mat_config_in=default_mat_config_out;
    
    num_vars=length(SP_var);
    
    for e3=1:num_vars;
        
        % Initialze the time data variable
        SP=[];
        
        if iscell(SP_var)
            
            if length(SP_var) >= e3
                SP=SP_var{e3};
            end

            if length(Fs_SP_var) >= e3
                Fs_SP=Fs_SP_var{e3};
            end

            % currently there is not any code for reverberation time of 
            % vibrations
            %if length(vibs_var) >= e7
            %    vibs=vibs_var{e7};
            %end
            %
            %if length(Fs_vibs_var) >= e7
            %    Fs_vibs=Fs_vibs_var{e7};
            %end
        else
            SP=SP_var;
            Fs_SP=Fs_SP_var;
            
            % currently there is not any code for reverberation time of
            % vibrations
            %vibs=vibs_var;
            %Fs_vibs=Fs_vibs_var;
        end
        
        [num_mics, num_pts]=size(SP);
        if num_mics > num_pts
            SP=SP';
            [num_mics, num_pts]=size(SP);
        end
        
        % Use the Fc array to input the center frequency of each recording file.
        fc=Fc(e1);

        dt=1/Fs_SP;
        t_SP=dt*(1:num_pts);
        
        title_str={Room_Name, 'For each Microphone', 'Level (dBA), T_6_0 (s), Quality', };
        rt_array=zeros(num_mics,1);     % buffer variable for the reverberation times
        
        lv_array=zeros(num_mics,1);     % buffer variable for the A-weighted signal level
        lv_array_BG=zeros(num_mics,1);  % buffer variable for the A-weighted background level
            
        if logical(num_mics > 0) && logical(num_pts > 0)
            
            count=count+1;
            mic_num_array(count)=num_mics;     % keep track of the number of mics in each file
            
            figure(1);
            delete(1);
            color_a={[0 0 1], [0 1 0], [1 0 0], [1 0 1], [0 1 1], [0.75 0.5 1], [0.4 0.2 0.3], [0.1 0.5 0.1],  [1 0.5 1], 0*[1 1 1], 0.4*[1 1 1], 0.7*[1 1 1], 0.9*[1 1 1]};
            Lya=zeros(num_mics, num_pts);   %
            rtpa=zeros(num_mics, num_pts);  %
            legend_str=cell(num_mics, 1);   % all the strings for the legend
            
            t_r=max(t_SP)-min(t_SP);
            
            for e2=1:num_mics;
                
                [rt, Ly, rtp]=reverb_time(SP(e2, :), Fs_SP, fc, t_high, t_low, num_x_filter, pbs, pab, method);
                
                Lya(e2, :)=Ly;
                rtpa(e2, :)=rtp;
                rt_array(e2)=rt;
                
                if method == 1
                    % Calculate the A-weighted signal level
                    
                    if isequal(Ba, 1) && isequal(Aa, 1)
                        [LeqA, LeqA8, LeqC, LeqC8, Leq, Leq8, Ba, Aa, Bc, Ac]=Leq_all_calc(SP(e2, 1:floor(Fs_SP*t_high)), Fs_SP, 1);
                    else
                        [LeqA, LeqA8, LeqC, LeqC8, Leq, Leq8, Ba, Aa, Bc, Ac]=Leq_all_calc(SP(e2, 1:floor(Fs_SP*t_high)), Fs_SP, 1, Ba, Aa, Bc, Ac);
                    end

                else
                    LeqA=20*log10(1./0.00002.*max(abs(SP(e2, :))));
                end
                
                lv_array(e2)=LeqA;

                % Calculate the A-weighted background level
                if isequal(Ba, 1) && isequal(Aa, 1)
                    [LeqA_BG, LeqA8_BG, LeqC_BG, LeqC8_BG, Leq_BG, Leq8_BG, Ba, Aa, Bc, Ac]=Leq_all_calc(SP(e2, [(num_pts-floor(Fs_SP*t_low)):end]), Fs_SP, 1);
                else
                    [LeqA_BG, LeqA8_BG, LeqC_BG, LeqC8_BG, Leq_BG, Leq8_BG, Ba, Aa, Bc, Ac]=Leq_all_calc(SP(e2, [(num_pts-floor(Fs_SP*t_low)):end]), Fs_SP, 1, Ba, Aa, Bc, Ac);
                end
                
                lv_array_BG(e2)=LeqA_BG;
                
                % Concatentate the legend string parts
                if rt < 0.001
                    legend_str{e2}=[ sprintf('NAN')];
                else
                    legend_str{e2}=[ sprintf('%5.1f', 0.1*round(10*LeqA)), '  ', sprintf('%-6.3f', 1/1000*round(1000*rt))];
                end

                if e2 == num_mics

                    % removal of outliers
                    rt_array_cal=rt_array(rt_array > 0);
                    gm=geomean2(rt_array_cal);  % calculate the geometric mean
                    stdrt=std(rt_array_cal);	% calculate standard deviation
                    
                    % find all values within num_std standard deviations of the geometric mean
                    if length(rt_array) > 3
                        min_lim=max([gm - num_std*stdrt 0]);
                        pts=intersect( find(rt_array > min_lim), find(rt_array < gm + num_std*stdrt));
                        
                        npts=setdiff( 1:num_mics, pts); % find the indices of the outlying data
                    else
                        pts=1:num_mics; % small data set keep all data
                        npts=[];        % small data no outliers
                    end
                    
                    outliers=zeros(size(rt_array));
                    if num_mics > 1
                        for e4=1:num_mics;
                            outliers(e4)=ismember(e4, npts);
                            if smd ~= 1 || num_mics > 1
                                if isequal(outliers(e4), 1)
                                    legend_str{e4}=[legend_str{e4}, ' \oslash '];
                                else
                                    legend_str{e4}=[legend_str{e4}, ' \surd '];
                                end
                            end
                        end
                    end

                    ptsa{count}=pts;   % Cell array keeps track of all the indices of acceptable reverberation times
                    nptsa{count}=npts; % Cell array keeps track of all the indices of outlying reverberation times

                    rts=rt_array(pts); % The set of accepted reverberation times

                    % Calculate the geometric mean, standard deviation, and
                    % 95% confidence interval for the accepted reverberation times
                    mn_rt=geomean2(rts); % calculate the geometric mean
                    stdrt=std(rts);  % calculate standard deviation
                    ci_int=call_ci(rts,95);  % calculate 95% confidence interval of the standard error of the t-distribution with a two-sided test

                    % calculate the mean level of the signal during t_high
                    meanlevel=round(10*log10(mean(10.^(lv_array(pts, 1)./10))));

                    % Find the index of the median reverberation time for plotting
                    % purposes
                    [mbuf ix]=min(abs(rt_array-mn_rt));

                    if smd == 1
                        legend_str={legend_str{ix}};
                    end

                    figure(1);
                    % round to the nearest factor of 5 for plotting
                    % purposes
                    if plot_all_mics == 1
                        ymax=5*ceil(max(max(0.2*Lya)));
                        ymin=5*floor(min(min(0.2*Lya)));
                    else
                        ymax=5*ceil(1+max(max(0.2*Lya(ix, :))));
                        ymin=5*floor(-1+min(min(0.2*Lya(ix, :))));
                    end

                    ylim([ymin ymax]);
                    hv=0.05*(ymax-ymin);

                    if plot_all_mics == 1
                        % Plot the processed time averaged level of the time
                        % record.
                        for e5=1:num_mics;
                            e4=mod(e5-1, 12)+1;
                            plot(t_SP, Lya(e5, :),  'color', color_a{e4}, 'Linewidth', 1);
                            hold on;
                        end
                        % plot the solution of the reverberation time calculation

                        for e5=1:num_mics;
                            e4=mod(e5-1, 12)+1;
                            plot(t_SP, rtpa(e5, :), 'color', color_a{e4}, 'Linewidth', 2);
                        end

                        if smd == 1
                            % plot the first row after the header
                            plot(min(t_SP)+t_r*[0.62 0.65], (ymax-0.20*(ymax-ymin)-hv*(1-1))*[1 1], 'color', color_a{1}, 'Linewidth', 2);
                        else
                            for e5=1:num_mics; % plot a whie background for the colors of the lines
                                fill(min(t_SP)+t_r*[0.61 0.66 0.66 0.61 0.61], (ymax-0.20*(ymax-ymin)-hv*(e5-1))+[1-hv/2 1-hv/2 1+hv/2 1+hv/2 1-hv/2], [1 1 1], 'EdgeColor', 'none');
                            end
                            
                            for e5=1:num_mics; % plot the colors of the lines
                                plot(min(t_SP)+t_r*[0.62 0.65], (ymax-0.20*(ymax-ymin)-hv*(e5-1))*[1 1], 'color', color_a{e5}, 'Linewidth', 2);
                            end
                        end
                    else
                        if smd == 1
                            % plot the median reverberation time
                            % add only the median mic data
                            plot(t_SP, Lya(ix, :), 'color', color_a{1}, 'Linewidth', 1);
                            hold on;
                            plot(t_SP, rtpa(ix, :), 'color', color_a{1}, 'Linewidth', 2);

                            fill(min(t_SP)+t_r*[0.61 0.66 0.66 0.61 0.61], (ymax-0.20*(ymax-ymin)-hv*(1-1))+[1-hv/2 1-hv/2 1+hv/2 1+hv/2 1-hv/2],[1 1 1], 'EdgeColor', 'none');
                            plot(min(t_SP)+t_r*[0.62 0.65], (ymax-0.20*(ymax-ymin)-hv*(1-1))*[1 1], 'color', color_a{1}, 'Linewidth', 2);
                        else
                            % plot the median reverberation time only add all of
                            % the mic data
                            plot(t_SP, Lya(ix, :), 'color', color_a{1}, 'Linewidth', 1);
                            hold on;
                            % plot the solution of the reverberation time calculations
                            plot(t_SP, rtpa(ix, :), 'color', color_a{1}, 'Linewidth', 2);
                            for e5=1:num_mics;
                                fill(min(t_SP)+t_r*[0.61 0.66 0.66 0.61 0.61], (ymax-0.20*(ymax-ymin)-hv*(e5-1))+[1-hv/2 1-hv/2 1+hv/2 1+hv/2 1-hv/2], [1 1 1], 'EdgeColor', 'none');
                            end
                            plot(min(t_SP)+t_r*[0.62 0.65], (ymax-0.20*(ymax-ymin)-hv*(ix-1))*[1 1], 'color', color_a{1}, 'Linewidth', 2);
                        end
                    end

                    % Give the pdf plot a title
                    % the time includes an estiamte of the reverberation time and
                    % if there are two or more mics then there will also be a two-sided 95%
                    % confidence interval
                    text(min(t_SP)+0.78*t_r, ymax-0.17*(ymax-ymin), title_str, 'Fontsize', 20, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'BackgroundColor', 'w' );

                    % Give the pdf plot a legend
                    if smd == 1
                        text(min(t_SP)+0.66*t_r, ymax-0.20*(ymax-ymin)-hv*(1-1), legend_str{1}, 'Interpreter','Tex', 'Fontsize', 19, 'FontName','FixedWidth', 'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle', 'BackgroundColor', 'w');
                    else
                        for e4=length(legend_str):(-1):1;
                            text(min(t_SP)+0.66*t_r, ymax-0.20*(ymax-ymin)-hv*(e4-1), legend_str{e4}, 'Interpreter','Tex', 'Fontsize', 19, 'FontName','FixedWidth', 'HorizontalAlignment', 'left', 'VerticalAlignment', 'middle', 'BackgroundColor', 'w');
                        end
                    end

                    set(gca, 'FontSize', 20);

                    xlabel('Time (Seconds)', 'Fontsize', 20);
                    ylabel('Instantaneous Amplitude ( dB ref. 20 \mu Pa)', 'Fontsize', 20);

                    if method == 1
                        title1=['Mean Level ', num2str(meanlevel), ' dBA, Freq = '];
                    else
                        title1=['Peak Level ', num2str(meanlevel), ' dB, Freq = '];
                    end
                    
                    if num_mics  > 1
                        title([title1, num2str(Fc(e1)), ' Hz, Reverb. Time = ', sprintf('%5.3f +- %5.3f (s)', [mn_rt ci_int])], 'Fontsize', 22);
                    else
                        title([title1, num2str(Fc(e1)), ' Hz, Reverb. Time = ', sprintf('%5.3f (s)', mn_rt)], 'Fontsize', 22);
                    end
                    
                    rt_array2{count}=rt_array;          % Collect the reverberation times
                    lv_array2{count}=lv_array;          % Collect the A-weighted signal levels
                    lv_array_BG2{count}=lv_array_BG;    % Collect the A-weighted background levels
                    cum_var_counter(count)=e1;          % Cumulative variable Counter
                    hold on;
                    ylim([ymin ymax]);

                    % Give the PDF plot a file name with the format
                    % Room_Name_frequency_Hz_Var#
                    if e3 > 1
                        filename1=[ Room_Name, '_', num2str(floor(Fc(e1))), '_Hz', '_', num2str(e3) ];
                    else
                        filename1=[ Room_Name, '_', num2str(floor(Fc(e1))), '_Hz' ];
                    end
                    % Save the PDF plot
                    save_a_plot_reverb_time(1, filename1);
                end
            end
        end
    end
end



% Delete the PDF Plot from the screen.
% Speeds up processing and cleans up the screen.
delete(1);

% Open Output file to save the Reverberation Times, A-weighted signal 
% Levels, and A-weighted background levels.
% Make sure outfile has a .txt extension
[filename_base, ext]=file_extension(outfile);
fid=fopen([filename_base '.txt'], 'w');

% % *********************************************************************
%Print Reverberation Times;
% % *********************************************************************
fprintf(fid, 'Reverberation Times\r\n');
fprintf(fid, 'Center Frequency\tMicrophone Number');

max_num_mics=max(mic_num_array);

tabs1='';
for e1=1:(max_num_mics);
    tabs1=[tabs1, '\t'];
end

fprintf(fid, tabs1);

if num_mics  > 1
    fprintf(fid, 'Mean\tSTD\t95%%CI\tOutliers\r\n');
else
    fprintf(fid, 'Mean\r\n');
end

fprintf(fid, 'Hz\t');

nums=1:(max_num_mics);
fprintf(fid, '%i\t', nums);
if num_mics  > 1
    fprintf(fid, '(s)\t(s)\t(s)\r\n');
else
    fprintf(fid, '(s)\r\n');
end


for e3=1:count;
    e1=cum_var_counter(count);
    pts=ptsa{e3};
    npts=nptsa{e3};
    rts=rt_array2{e3}(pts);
    rts_cal=rts(rts>0);
    mn_rt=geomean2(rts_cal);
    std_rt=std(rts_cal);
    ci_int=call_ci(rts_cal,95);

    fprintf(fid, '%i\t', floor(Fc(e3)));
    fprintf(fid, '%s', sprintf('%6.3f\t', 0.001*round(1000*[rt_array2{e3}])) );

    tabs1='';
    for e2=1:(max_num_mics-mic_num_array(e3));
        tabs1=[tabs1, '\t'];
    end

    fprintf(fid, tabs1);
    if num_mics  > 1
        fprintf(fid, '%s', sprintf('%6.3f\t', 0.001*round(1000*[mn_rt, std_rt, ci_int])) );
        for e2=1:length(npts);
            fprintf(fid, '%d\t', npts(e2) );
        end
    else
        fprintf(fid, '%s', sprintf('%6.3f\t', 0.001*round(1000*mn_rt)) );
    end

    fprintf(fid, '\r\n');

end

fprintf(fid, '\r\n\r\n\r\n');
% % *********************************************************************
%Print Signal High Levels;
% % *********************************************************************
if method == 1
    fprintf(fid, 'Time Averaged Speaker on Sound Pressure Levels (dBA)\r\n');
else
    fprintf(fid, 'Peak Sound Pressure Levels (dB)\r\n');
end

fprintf(fid, 'Center Frequency\tMicrophone Number');

max_num_mics=max(mic_num_array);

tabs1='';
for e1=1:(max_num_mics);
    tabs1=[tabs1, '\t'];
end

fprintf(fid, tabs1);
if num_mics  > 1
    fprintf(fid, 'Mean\tSTD\t95%%CI\tOutliers\r\n');
else
    fprintf(fid, 'Mean\r\n');
end

fprintf(fid, 'Hz\t');

nums=1:(max_num_mics);
fprintf(fid, '%i\t', nums);
if method == 1
    if num_mics  > 1
        fprintf(fid, '(dBA)\t(dBA)\t(dBA)\r\n');
    else
        fprintf(fid, '(dBA)\r\n');
    end
else
    if num_mics  > 1
        fprintf(fid, '(dB)\t(dB)\t(dB)\r\n');
    else
        fprintf(fid, '(dB)\r\n');
    end
end


for e3=1:count;
    e1=cum_var_counter(count);
    pts=ptsa{e3};
    npts=nptsa{e3};
    mn_rt=10*log10(mean(10.^(lv_array2{e3}(pts)./10)));
    std_rt=std(lv_array2{e3}(pts));
    ci_int=call_ci(lv_array2{e3}(pts),95);
    
    fprintf(fid, '%i\t', floor(Fc(e3)));
    fprintf(fid, '%s', sprintf('%6.1f\t', 0.1*round(10*[lv_array2{e3}])) );

    tabs1='';
    for e2=1:(max_num_mics-mic_num_array(e3));
        tabs1=[tabs1, '\t'];
    end

    fprintf(fid, tabs1);

    if num_mics  > 1
        fprintf(fid, '%s', sprintf('%6.1f\t', 0.1*round(10*[mn_rt, std_rt, ci_int])) );
        for e2=1:length(npts);
            fprintf(fid, '%d\t', npts(e2) );
        end
    else
        fprintf(fid, '%s', sprintf('%6.1f\t', 0.1*round(10*mn_rt)) );
    end

    fprintf(fid, '\r\n');

end

fprintf(fid, '\r\n\r\n\r\n');
% % *********************************************************************
%Print Background Levels;
% % *********************************************************************
fprintf(fid, 'Time Averaged Background Sound Pressure Levels dBA\r\n');
fprintf(fid, 'Center Frequency\tMicrophone Number');

max_num_mics=max(mic_num_array);

tabs1='';
for e1=1:(max_num_mics);
    tabs1=[tabs1, '\t'];
end

fprintf(fid, tabs1);
if num_mics  > 1
    fprintf(fid, 'Mean\tSTD\t95%%CI\tOutliers\r\n');
else
    fprintf(fid, 'Mean\r\n');
end

fprintf(fid, 'Hz\t');

nums=1:(max_num_mics);
fprintf(fid, '%i\t', nums);
if num_mics  > 1
    fprintf(fid, '(dBA)\t(dBA)\t(dBA)\r\n');
else
    fprintf(fid, '(dBA)\r\n');
end

for e3=1:count;
    e1=cum_var_counter(count);
    pts=ptsa{e3};
    npts=nptsa{e3};
    mn_rt=10*log10(mean(10.^(lv_array_BG2{e3}(pts)./10)));
    std_rt=std(lv_array_BG2{e3}(pts));
    ci_int=call_ci(lv_array_BG2{e3}(pts),95);
    
    fprintf(fid, '%i\t', floor(Fc(e3)));
    fprintf(fid, '%s', sprintf('%6.1f\t', 0.1*round(10*[lv_array_BG2{e3}])) );

    tabs1='';
    for e2=1:(max_num_mics-mic_num_array(e3));
        tabs1=[tabs1, '\t'];
    end

    fprintf(fid, tabs1);

    if num_mics  > 1
        fprintf(fid, '%s', sprintf('%6.1f\t', 0.1*round(10*[mn_rt, std_rt, ci_int])) );
        for e2=1:length(npts);
            fprintf(fid, '%d\t', npts(e2) );
        end
    else
        fprintf(fid, '%s', sprintf('%6.1f\t', 0.1*round(10*mn_rt)) );
    end

    fprintf(fid, '\r\n');

end

% % *********************************************************************
% % Close the text file
% % *********************************************************************
fclose(fid);
fclose('all');
