function [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low, num_x_filter, pbs, pab, method)
% % Syntax:
% % [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low, num_x_filter,pbs, pab, method);
% % 
% % This program supports two methods of calculating the reverberation
% % time: the speaker on speaker off method and the balloon pop method.
% % The speaker on speaker off method is considered to be much more
% % accurate.  It is highly recommended to use the speaker on speaker off
% % method if possible.
% % 
% % This program calculates reverberation times for 1/3 octave broad 
% % band noise for multiple microphones.  See the description of the
% % variable fc.
% %   
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Description of the Input Variables
% % 
% % SP is the time record in units Pa and has size (num_mics,
% %     num_samples).  SP is automatically rearranged so that the number  
% %     of elements in the second dimension are greater than the first  
% %     dimension. 
% % 
% % Fs (Hz) is the samping rate.
% % 
% % fc (Hz) is the third octave band center frequency.  
% % 
% %     Filters have anomalies and program returns meaningless data 
% %     for fc > Fs/3.  
% %     
% %     It is difficult to excite source below 31 Hz with enough energy   
% %     for a meaningful reverberation time.
% %     
% %     For the balloon pop method, the impulsive reverberant decay must 
% %     reach the background level.  
% %     
% %     For the Balloon pop method measurement times of 10 seconds or
% %     longer are recommended for fc < 200 Hz.  
% %       
% % 
% % t_high (s) is used in the speaker on spekaer off method.
% %     t_high is the duration of signal on at the beginning of the time 
% %     record.  This is used to calculate the average signal level.
% % 
% % t_low (s) is used by both methods.  t_low is the duration of quiet 
% %     time at the end of the time record.  t_low is used to calculate the 
% %     average background level.  
% % 
% % num_x_filter is the number of times to filter the time record. 
% %     The minimum value of num_x_filter is 1.
% %     Typically a value of 2 is sufficient to remove the contributions 
% %     from other third octave bands.  
% % 
% %     If there is spurious background noise, espcially in the low
% %     frequencies fc < 100 Hz, then additional filtering may be
% %     beneficial.  num_x_filter=10 can remove up to 120 dB of noise.  
% %     Unfortunately the third octave band filter changes the calculated 
% %     reverberation times, so be careful.  
% % 
% %     Increasing num_x_filter does not necessarily improve accuracy.
% %     Increasing the input signal speaker level or balloon pop level 
% %     may improve the accuracy for both methods.
% %     
% % pbs is an acronym for the percentage of the total range in dB below 
% %     the average signal level.  pbs is used to discard spurious 
% %     high level sound data at the speaker cutoff. pbs=20 usually works.
% % 
% % pab is an acronym for the percentage of the range above the average  
% %     background level.  pab is used to discard spurious background 
% %     noise where the reveberant decay approaches the background level.
% %     This is more of a problem in environments with very low background 
% %     For example, measuring reverberation time in a hemi-anechoic 
% %     chamber.
% % 
% % method selects the speaker on speaker off method or the balloon pop
% %     method.  1 for speaker, otherwise balloon pop.
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % Description of the Output Variables
% %
% % rta     % The array of the reverberation times
% %
% % Lya     % The array of processed time records
% %
% % rtpa    % The array of time records of the solution of the 
% %         % reverberation time calculation
% %
% % SP2     % Resampled and filtered time record
% % 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% 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.
% 
% 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];
% 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);];
% 
% figure(1); for e1=1:4; subplot(4, 1, e1); plot(SP(e1, :)); hold on; end;
%
% Fs=50000;         % (Hz) Sampling rate
% 
% fc=400;           % (Hz) Center frequency of the third-octave band
% 
% t_high=0.5;       % (s) duration of loud time at the beginning of the time record
%                   % used to calculate the average signal level,
%                   % which is usead to estimate the speaker cutoff time
%                   % A longer t_high time is better when there is filter phase delay
%                   % especially at low frequencies (Fc < 100) or
%                   % if the speaker turns onslowly
% 
% t_low=0.5;        % (s) duration of quiet time at the end of the time record
%                   % used to calculate the average background level for
%                   % calculating the time when the background noise level is reached
%                   % a longer time gets the background level above
%                   % the dips and peaks in the background level from
%                   % sudden quiet and background dog barks etc...
% 
% num_x_filter=2;   % Number of times to filter the data.  Minimum value is 1
%                   % typically a value of 2 is sufficient to remove noise.
%                   % num_x_filter=10 can remove 120 dB of noise.
%                   % Increasing num_x_filter changes the cacualted 
%                   % reverberation time.
% 
% pbs=15;   % Percentage of the total range in dB below the average signal
%           % level to discard from fitting to a line.  This avoids
%           % spurious high level sound data at the speaker cutoff.
% 
% pab=40;   % Percentage of the range above the average background level
%           % to discard.  Increases the background cutoff level above 
%           % spurious background noise where the reveberation decay
%           % approaches the background level.  A value of 40 worked in
%           % a laboratory hemi-anechoic chamber and may work for 
%           % acoustical environments with very low background noise levels.
% 
% pab=15;   % A value of 15 worked for an animal shelter with dogs barking
%           % in the surrounding rooms.  
% 
% method=1; % 1 is the speaker on speaker off method.  
%           % otherwise is the balloon pop method
% 
% [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low,num_x_filter,pbs,pab, method);
% 
% % Exact reverberation times are
% t_60=log(1000).*tau;
%
% % The percent error is
% Percent_error=100.*(1-(t_60./rta'));
% %
% %
% 
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% Example=''; % Simulated reverberation data 
%             % For the Balloon Pop 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.
% 
% 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];
% 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-tenth of a second of background noise
% %     2) reverberant decay
% %     3) nominal one second of background noise
% 
% SP=zeros(4, 2*Fs+length(t));
% SP(1, :)=[0.01*A(1)*randn(1,0.1*Fs) decay1 A(1).*exp(-1./tau(1).*tb)*randn(1, Fs);];
% SP(2, :)=[0.01*A(2)*randn(1,0.1*Fs) decay2 A(2).*exp(-1./tau(2).*tb)*randn(1, Fs);];
% SP(3, :)=[0.01*A(3)*randn(1,0.1*Fs) decay2 A(3).*exp(-1./tau(3).*tb)*randn(1, Fs);];
% SP(4, :)=[0.01*A(4)*randn(1,0.1*Fs) decay2 A(4).*exp(-1./tau(4).*tb)*randn(1, Fs);];
% 
% Fs=50000;        
% fc=400;           
% t_high=0.5;  
% t_low=0.5;  
% num_x_filter=2;   
% pbs=15;   
% pab=15;   
% method=2; 
% 
% [rta, Lya, rtpa]=reverb_time(SP, Fs, fc, t_high, t_low,num_x_filter,pbs,pab, method);
% 
% % Exact reverberation times are
% t_60=log(1000).*tau;
%
% % The percent error is
% Percent_error=100.*(1-(t_60./rta'));
% 
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% % 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
% % findextrema From Schuberth Schuberth found on matlab central
% % 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
% %
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% % Program was written by Edward L. Zechmann
% % 
% %      date 16 October   2007
% %
% %  modified 13 January   2008     updated comments
% %                                 changed from filter.m to filtfilt.m
% %                                 Added impulsive nosie method
% % 
% %  modified 14 January   2008     Reduced number of points allowed for
% %                                 the robust line fit
% %                                 updated comments
% %
% % modified 13 February   2008     Added the balloon pop method and
% %                                 example
% %                                 Made the fit adaptive for both methods
% %                                 Added a leverage point to the 
% %                                 Speaker on speaker off method when the 
% %                                 decay line fit does not properly
% %                                 intersect the signal high time record
% % 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% % Please feel free to modify this code.
% %

if nargin < 9
    method=1;
end

fc=round(fc);
[SP]=convert_double(SP);

[m1, n1]=size(SP);

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

if nargin < 6
    num_x_filter=1;
else
    if num_x_filter < 1
        num_x_filter=1;
    end
end

dt=1/Fs;

%Resample timeseries so that filter has better performance

p=1;
q=1;
flag1=0;
Fsn=Fs;

if fc < Fs/50
    Fsn=40*fc;
    flag1=1;
elseif fc > Fs/10
    Fsn=10*fc;
    flag1=1;
end

if flag1 == 1

    dtn=1/Fsn;
    [p,q] = rat(dt/dtn);
    buf=resample(SP(1, :), p, q);
    n2=length(buf);
    SP2=zeros(m1,n2);
    SP2(1, :)=buf;
    clear('buf');
    if m1 > 1
        for e1=2:m1;
            SP2(e1, :)=resample(SP(e1, :), p, q);
        end
    end
    Fs2=Fs*p/q;
else
    Fs2=Fs;
    SP2=SP;
end
% Calculate the filter coefficients for the third octave filter
[Bc,Ac] = oct3dsgn(fc,Fs2,3);
clear('SP');

for e1=1:m1;
    for e2=1:num_x_filter;
        SP2(e1, :)=filtfilt(Bc,Ac,SP2(e1, :)); %Process data with a 1/3 octave bandpass filter
    end
end

if flag1 == 1

    buf=resample(SP2(1, :), q, p);
    SP=zeros(m1,length(buf));
    SP(1, :)=buf;
    clear('buf');

    if m1 > 1
        for e1=2:m1;
            SP(e1, :)=resample(SP2(e1, :), q, p);
        end
    end

else
    SP=SP2;
    Fs=Fs2;
end

clear('SP2');

% return SP to size [m1, n1]
% check size of SP zero pad if necessary
[m2 n2]=size(SP);


if n1 ~= n2
    if n2 > n1
        SP=SP(1:m1, 1:n1);
    else
        a=min(min(abs(SP)));
        for e1=1:m1;
            for e2=1:n1;
                if e2 <= n2
                    SP(e1, e2)=SP(e1, e2);
                else
                    SP(e1, e2)=a;
                end
            end
        end
    end
end

[m1 n1]=size(SP);
SP=SP-mean(SP, 2)*ones(1, n1);
a=1:m1;


SP = abs(hilbert(SP')');

I_high2=floor(t_high/dt);
I_low2=floor(t_low/dt);

rta=zeros(m1,1);
Lya=zeros(m1,n1);
rtpa=zeros(m1,n1);

for e1=1:m1;

    if isequal(method, 1)
        % Process data for the speaker on speaker off method
        [buf1, y22 ]=sub_mean(SP(e1, 1:n1), Fs, 25);
    else
        % Process data for the balloon pop method
        [buf1, y22 ]=sub_mean(SP(e1, 1:n1), Fs, 100);
    end

    Ly=20*log10(y22(1, :)./(20*10^(-6)));
    
    % reset the values for I_high and I_low adn the other quantities
    I_high=I_high2;
    I_low=I_low2;
    IX_est=0;
    IX3=n1;
    count=0;
    t1=0;
    t2=n1-I_low;
    IX=1;
    flag=0;
    
    while (logical(abs(IX3-IX_est)/IX3 > 0.1) || logical(abs(I_high-t1)/I_high > 0.1) || logical(abs(I_low-(n1-t2))/I_low > 0.1)) && logical(count < 3)
                
        if isequal(count, 0)
            IX3=n1-I_low;
        end
        count=count+1;
        I_high=max([I_high, IX]);
        buf1=[n1-t2, n1-IX3];
        buf1=buf1(buf1>0);
        I_low=min(buf1);
        
        
        if I_low > n1
            I_low=n1-1;
        end
        
        if I_low < 1
            I_low=1;
        end
        
        if I_high > n1
            I_high=n1-1;
        end
        
        if I_high < 1
            I_high=1;
        end
        
        if isequal(method, 1)
            % Level at beginning
            Lb=sqrt(sum(y22(1, 1:I_high)'.^2)./I_high);
        else
            % Level of the balloon pop
            Lb=max(y22(1, :));
        end

        SPLb=20*log10(Lb./(20*10^(-6)));

        % Level at the end
        Le=sqrt(sum(y22(1, (n1-I_low):n1)'.^2)./I_low);
        SPLe=20*log10(Le./(20*10^(-6)));

        % Find the beginning and ending of the decay
        % Set cutoff level pbs percent of dB range below the signal level
        % Set cutoff level pab percent of dB range above the background level

        % Find last data points where the level is below the signal level or
        % below the balloon pop level
        if isequal(method, 1)
            IX=find( Ly > SPLb-pbs/100*(SPLb-SPLe), 1, 'last' );

            % find data points where Ly is below the background level
            IX2=find( Ly < SPLe+pab/100*(SPLb-SPLe) );

            % find the first point below the balloon pop level, but above the
            % background level
            IX3=IX2(find(IX2 > IX, 1 ));
        else
            [max_LY IX]=max(Ly);
            IX2=find( Ly < SPLe+pab/100*(SPLb-SPLe) );
            IX3=IX2(find(IX2 > IX, 1 ));
        end

        % Make sure there are plenty of data points between IX and IX3
        if IX3 - IX < 1000;
            if IX3+1000 > n1
                IX=IX-1000;
            else
                IX3=IX3+1000;
            end
        end
        
        
        if length(IX:IX3) > 5

            if isequal(method, 1)
                % SPeaker On Speaker Off Method
                IX3=min([max([IX3, IX_est]), n1]);
                
                if IX3 > n1
                    IX3=n1-1;
                end
                
                if IX3 < 1
                    IX3=1;
                end
                
                rtd=Ly(1, IX:IX3);
                tt=dt*(IX:IX3);
                [LMSout,blms,Rsq]=LMTSreg(rtd',tt', 1000, 100000);
                p=[blms(2) blms(1)];
                t2=floor(Fs*(SPLe-p(2))/p(1));
                IX_est=t2;

            else
                
                % Balloon Pop Method
                IX3=max(IX3, IX_est);
                
                if IX3 > n1
                    IX3=n1-1;
                end
                
                if IX3 < 1
                    IX3=1;
                end
                
                rtd=Ly(1, IX:IX3);
                tt=dt*(IX:IX3);
                num_pts=length(rtd);
                num_bins=min([num_pts/5, 1000]);
                rbp=floor(num_pts/num_bins);
                bins=(1:rbp:num_pts);
                epts=zeros((length(bins)-1), 1);
                
                for e2=1:(length(bins)-1);
                    [buf, eptsix]=max(rtd(bins(e2):bins(e2+1)));
                    epts(e2)=bins(e2)-1+eptsix;
                end
                
                [ma,mi]=findextrema(rtd);
                rtd_n1=length(rtd);
                ma=unique(round(ma));
                ma=intersect(ma, epts);
                
                if length(ma) > 10
                    if ma(1) < 1
                        ma(1)=1;
                    end

                    if ma(end) > rtd_n1
                        ma(end)=rtd_n1;
                    end

                    ma=unique(ma);

                    [LMSout,blms,Rsq]=LMTSregor(-max_LY+rtd(ma)',-tt(1)+tt(ma)', 1000, 100000);
                    p=[blms(1) max_LY];
                else
                    [LMSout,blms,Rsq]=LMTSregor(-max_LY+rtd',-tt(1)+tt', 1000, 100000);
                    p=[blms(1) max_LY];

                end

                t2=floor(Fs*(SPLe-p(2))/p(1));
                IX_est=IX+t2;

            end
        else
            count=1000000;
            p(1)=60;
            p(2)=SPLb;
            rtp=SPLb*ones(1, n1);
        end
        
        t1=ceil(Fs*(SPLb-p(2))/p(1)); % first intersection point
        t2=floor(Fs*(SPLe-p(2))/p(1)); % second intersection point
        
        if ((-t1+IX)/IX > 0.3 || logical(p(1) > 0)) && isequal(method, 1)
            % Method failed
            % a leverage point is necessary to ensure a good fit
            % The intesection point of the signal high cutoff limit and 
            % the instantaneous sound level data will be the leverage point
            flag=1;
            count=4;
            [LMSout,blms,Rsq]=LMTSregor(-Ly(IX)+rtd',-dt*(IX)+tt', 1000, 100000);
            p=[blms(1) Ly(IX)];
            
            t1=ceil(Fs*(SPLb-p(2))/p(1)); % first intersection point
            t2=floor(Fs*(SPLe-p(2))/p(1)); % second intersection point
        
        end      
        
        rty = polyval(p,dt*(t1:t2));

        if isequal(method, 1) && ~isequal(flag, 1)
            rtp=[ SPLb*ones(1, t1-1) rty SPLe*ones(1, (n1-t2))];
        else
            rtp=[ SPLb*ones(1, IX+t1-1) rty SPLe*ones(1, (n1-IX-t2))];
        end
        
    end
    
    rt=-60/p(1);                % Estimated Reveberation Time
    rta(e1)=rt;                 % Array of estimated reverberation time measurements
    Lya(e1, :)=Ly;              % Time record of processed time record
    rtpa(e1, :)=rtp(1, 1:n1);   % Time record of reverberation time solution
end




