function [T, fo, X, P, E, I] = iso3382(ximp, Fs, k, d)
    // Reverberation time by Schroeder's reverse integration
    // method specified in International Standard ISO 3382
    //
    // Usage:
    //       [T, fo, X, P, E, I] = iso3382(ximp, Fs, k, d)  
    //
    // where
    //       ximp: Impulse response as a row vector 
    //       Fs:   Sampling rate in Hz (default is 44100 Hz)
    //       k:    Octave fraction denominator (default is 1)
    //       d:    Useful dynamic range of decay in dB (default
    //             is 30)
    //       SN:   Signal-to-noise ratio in dB, either a scalar,  
    //             a vector with a value for each band's frequency
    //             or a sample of about 1 s of background noise
    //             (default is 45)
    //       T:    Vector with octave-fraction band reverberation
    //             times in s 
    //       fo:   Octave-fraction band center frequency
    //       X:    Matrix whose rows represent octave-fraction 
    //             band-filtered versions of ximp
    //       P:    Matrix whose rows contain the coefficients of
    //             the linear regression of the -5 dB to -35 dB 
    //             decay of the corresponding rows of E  
    //       E:    Matrix whose rows represent reverse-integrated
    //             square filtered impulse response
    //       I:    Matrix whose rows contain the indices where
    //             the corresponding rows of E reach -5 dB and
    //             -35 dB below their maximum (or (-5 + d) dB)
    //
    // NOTES: 
    // 1  Requires fftfilter and bandfilter
    // 2  Zero padding is added at the beginning to prevent
    //    signal loss due to non-causal transient response of
    //    the FFT band filtering 
    // 3  The script converts automatically to row vector the
    //    input if it is presented as column vector, which is 
    //    the case when ximp has been imported by wavread
    // 4  k is optional. Possible values are 0 (full band),
    //    1, 2 and 3 (octave, half octave and one third octave). 
    //    Default: 1
    // 5  SN is optional, its default value being 45. The
    //    function features smart detection of data type
    // 6  d is optional, its default value being 30. It may be
    //    reduced to 20 and still comply with ISO 3382. This
    //    is useful when there is too much background noise in
    //    one or more bands.  
    // 7  Sampling rate is optional, its default value being
    //    44100
    // 8  All output arguments except T are optional. 
    // 
    // REFERENCES
    // Schroeder, M. (1965) New method of measuring reverberation
    //    time. JASA, 37, 409–412 (1965).
    //    https://doi.org/10.1121/1.1909343
    // ISO 3382-1 (2009) Acoustics — Measurement of room acoustic 
    //    parameters. Part 1: Performance spaces
    //
    // ----------------------------------------
    // Author: Federico Miyara
    // Date:   2013-10-17
    //         2017-06-06
    //         2020-03-26
    //         2024-08-14 

    // Number of input arguments
    rhs = argn(2)

    // Default input arguments if some are not provided
    if rhs<5
        SN = 45;
    end
    if rhs<4
        d = 30;
    end
    if rhs<3
        k = 1;
    end
    if rhs<2
        Fs = 44100;
    end

    // Load band-filtering function
    exec("bandfilter.sci");

    // Conversion to row vector if necessary
    ximp = ximp(:)';

    // Zero padding at the beginning to prevent transient 
    // issues especially at low-frequency bands
    x = [zeros(1,32768), ximp];
    t = (0:length(x)-1)/Fs;

    // Type of spectral band
    select k
    case 0
        // Full audio band without filtering
        fo = Fs/2;
    case 1
        // Octave band frequencies 
        fo = [63  125  250  500  1000  2000  4000  8000  16000]';
    case 2
        // Half octave band frequencies
        fo = [22.4  31.5  45  63  90  125  180  250  355  500  710 ...
              1000  1400  2000  2800  4000  5600  8000  11200  16000]';
    case 3
        // One third octave band frequencies
        fo = [ 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]';
    end

    // Preallocation of filtered signal matrix
    X = zeros(length(fo), length(x));

    // Band filtering of impulse response
    for h=1:length(fo)
        if k==0
            X = x;
        else
            if k==0
                X(h,:) = x;
            else
                X(h,:) = bandfilter(x, fo(h), k, Fs);
            end
        end 
    end

    // Reverberation time for each band
    for h=1:length(fo)
        [TT, pp, EE, ii] = schroeder(X(h,:), Fs, d);
        T(h,:) = TT;
        P(h,:) = pp;
        E(h,:) = EE;
        I(h,:) = ii;
    end

endfunction


function [T, p, E, i] = schroeder(x, Fs, d)
    // Reverse integration, linear fitting and 
    // reverberation time computacion
    // 
    // Usage:
    //       [T, p, E] = schroeder(x, Fs, d)
    // where
    //       x:  Impulse response vector
    //       Fs: Sample rate in Hz
    //       d:  Useful dynamic range in dB
    //       T:  Reverberation time
    //       p:  Parameters of the linear fit [a,b]
    //           Y = a + b*(t-to)
    //       E:  Energy reverse integration
    //       i:  Two-component vector with the indices
    //           where E is -5 dB and -35 dB below its
    //           maximum 
    //
    // ---------------------------
    // Federio Miyara
    // 2013-10-17
    // 2017-06-06
    // 2018-07-05
    // 2024-08-13

    L = length(x);
    // Reverse integration    
    E = cumsum(x($:-1:1).^2);
    E = E($:-1:1);
    
    // Detect last sample that is greater than 5 dB
    // less than the maximum (which is always the 
    // first sample)
    i(1) = max(find(E/E(1)>10^(-10/20)));
    t1 = i(1)/Fs;
     
    // Detect the last sample that is grater than 
    // 5 dB + d dB less than the maximum
    i(2) = max(find(E/E(1)>10^(-(10 + 2*d)/20)));
    t2 = i(2)/Fs;
     
    // Linear regression fitting of the energy decay
    // Gap function (error to minimize)
    function e = G(p, z)
       x = z(1,:);
       y = z(2,:);
       e = y - (p(1) + p(2)*(x - t1));
    endfunction
    
    // Empirical data
    Z = [t(i(1):i(2)); 10*log10(E(i(1):i(2)))];
    // Initial guess of the parameters 
    a = 10*log10(E(i(1)));
    b = 10*log10(E(i(2))/E(i(1)))/(t2 - t1);
    po = [a, b]; 
    // Optimum values of the parámetros
    [p, err] = datafit(G, Z, po);
    
    // Estimated reverberation time
    T = -60/p(2);
endfunction


function y = rms(x)
    // RMS value
    // 
    // Usage:
    //       y = rms(x)
    // where
    //       x:  Vector containing the discreta signal
    //       y:  RMS value of x
    
    y = sqrt(sum(x.^2)/length(x)) 
endfunction

