function [Lx2,Lx,Lmin,Lmax,fo] = octavemean(x,D,Fs,N,L,cal,dB,w,graph)
    // Octave-band analysis of calibrated audio signal
    // This procedure computes the equivalent mean octave
    // band spectrum of an audio signal x. In order to 
    // get meaningful results it is necessary to provide
    // calibration information. Argument dB refers to the
    // sound pressure level of the calibration signal
    // (usually 94.0 dB), while cal refers to the level
    // corresponding to such signal referenced to maximum
    // signal (usually -12 dB to -24 dB depending on the
    // desired headroom, general calibration of audio
    // chain and intended maximum level). Usually x will
    // have been read from an uncompressed audio file in
    // WAV format using the function wavread, so maximum
    // signal will be +/- 1.
    //
    // Usage:
    //     Lx2 = octavemean(x,D,Fs,N,L,cal,dB,w,graph)
    //     [Lx2,Lx,Lmin,Lmax,fo] = octavemean(x,D,Fs,cal,dB,w,graph)
    //
    // where
    //     x:     Audio signal column or row vector or
    //            array of columns 
    //     D:     Denominator of fractional octave band
    //            i.e., 1 for octave, 3 for 1/3 octave, 
    //            etc. Custom bands provided as a vector
    //            of consecutive frequency limits are also 
    //            allowed.     
    //     Fs:    Sampling rate in Hz
    //     N:     Length of the analysis window
    //     L:     A positive number indicating the offset
    //            in samples between successive analysis
    //            windows or a list whose first element
    //            is a string indicating the method of
    //            selection of the windows to average
    //            ("resample", "uniform", or "poisson",
    //            see Notes) and whose second element
    //            indicates the number of windows to
    //            select for averaging.
    //     cal:   Level of calibration tone referenced to
    //            maximum possible sinusoidal signal (+/- 1)
    //            in dBFS (AES17) 
    //     dB:    Sound pressure level of calibration
    //            tone in dB
    //     w:     Windowing function to apply to the signal. 
    //            It may be a custom window given as a vector
    //            with N components or a string naming a
    //            supported window:
    //               "box":   Boxcar window (no window)
    //               "hann":  Periodic Hann window
    //               "black": Periodic Blackman window 
    //     graph: Determines graphic output
    //            If it is a string: 
    //               'fig': new figure
    //               'cuf': current figure 
    //               'nof': no figure
    //            If it is a list, the first element
    //            is a string as in the preceding case 
    //            and the next elements are the indices
    //            of the spectrum to plot
    //     Lx2:   Column vector or array of column vectors 
    //            with octave band energy-averaged spectrum
    //            of the input signal in dB
    //     Lx:    Column vector or array of column vectors 
    //            with octave band linearly averaged spectrum
    //            of the input signal in dB
    //     Lmax:  Column vector or array of column vectors
    //            with octave band maximum spectrum of the
    //            input signal in dB
    //     Lmin:  Column vector or array of column vectors 
    //            with octave band minimum spectrum of the
    //            input signal in dB
    //     fo:    Column vector containing the standard
    //            octave or fractional band frequencies or 
    //            the geometric mean of consecutive 
    //            frequencies in case of custom bands
    //
    // NOTES:
    // Arguments D, Fs, N, L, cal, dB, w and graph are
    // optional;  their default values are, respectively,
    // 1, 44100, 8192, 8192, 0, 94.0, "hann" and 'nof'.
    // When L is a list, the methods of window selection
    // are:
    // "resample": the windows are selected randomly from
    //             the whole length of the signal (a
    //             statistical method called bootstrap 
    //             resampling). Depending on the number 
    //             of windows to be selected, the windows
    //             might overlap, repeating parts of the
    //             signal, or might leave parts of the 
    //             signal unanalyzed.
    // "uniform":  The windows are selected so that they
    //             stick as tightly as possible to uniform
    //             spacing.
    // "poisson":  The analysis windows are located randomly 
    //             with a Poisson's distribution so that the
    //             first and last window cover the beginning 
    //             and end of the signal.
    //               
    // If the graph argument is a string different from "nof"
    // the output is provided only for the first signal, i.e.,
    // x(:,1, ..., 1). If it is a list such as ["cuf",3,2]
    // then the spectrum to be plotted is that of x(:,3,2). 
    // If the minimum spectrum is 0, then it will turn -inf
    // when converted to dB and consequently it will not show
    // on the graph.
    //
    // -------------------------------------------------------------
    // Author: Federico Miyara
    // Date:   2005-12-05
    //         2009-08-25
    //         2010-06-07
    //         2015-12-30
    //         2020-01-21
    //         2020-02-19
    //         2023-12-04 


    // INPUT DATA --------------------------------------------------
    
    // Number of input (right-hand side) arguments
    rhs = argn(2);
 
    // Default values of arguments if missing
    if rhs<2
        // Octave band (the most common one)
        D = 1;
    end
    if rhs<3
        // Usual CD-quality sampling rate  
        Fs = 44100;
    end
    if rhs<4 
        // Select a window length to ensure accuracy
        // at low frequency. The smallest band must 
        // contain at least 3 spectral lines
        if isvector(D)
            N = 2^nextpow2(3*Fs/min(diff(D)));
        else
            N = 2^nextpow2(8192*D);     
        end
    end
    if rhs<5
        L = N; 
    end
    if rhs<6
        // Assumption that cal corresponds to reference 
        // wav level (+/-1)
        cal = 0;
    end
    if rhs<7
        // Standard sound pressure level of calibration tone
        dB = 94.0;
    end
    if rhs<8
        // Use Hann window
        w = "hann";
    end
    if rhs<9
        // No graphic output if not expressly asked for
        graph = 'nof';
    end


    if 1==2
        // --------------------------------
        // TEST-CASE DATA - DELETE
        // --------------------------------
        exec pinknoise.sci;
        Fs = 44100;
        D = 1;
        N = 8192;
        L = 4096;
        len = 100000;
        x = zeros(len, 2, 4);
        siz = size(x);
        for h=1:siz(3)
            // h-th sound 
            for j=1:siz(2)-1 
                // j-th channel
                x(1:len,j,h) = pinknoise((len-1)/Fs, Fs, N);
            end
            x(1:len,siz(2),h) = sin(2*%pi*1000*[0:len-1]/Fs);
        end
        cal = 0;
        dB = 94;
        w = "black";
        graph = "fig";
        //--------------------------------
    end

    // Load required functions
    prot = funcprot();
    funcprot(0);
    exec octaves.sci;
    funcprot(prot)

    // Size of input array
    siz = size(x);
    
    // If the argument is a row vector (or else it is an N-D array
    // made of row vectors, meaning that the first --or column--
    // dimension is 1), the first two dimensions are permuted so
    // that it turns to be a column vector (or an N-D array made of
    // column vectors). Its size is then recalculated. The result 
    // is a series of layers, each one containing column vectors or
    // matrices formed by column vectors. An FFT will be performed
    // on these vectors.    
    if siz(1)==1
        x = permute(x, [2, 1, [3:length(siz)]]);        
    end
    
    // New signal size
    siz = size(x);

    // Length of each column containing signal
    len = siz(1);
       
    // L is a list whose first element indicates the method
    // for choosing the windows to average and whose second
    // element indicates the number of windows to average
    if type(L)==15
        // Number of FFT windows
        m = L(2);
        // Number of complete windows
        q = m;
        // Length of incomplete window
        H = 0;
        // Sampling type
        select L(1)
        case "resample"
            // Using bootstrap resampling
            // Windows are taken randomly with uniform
            // distribution along the signal duration.
            // They may randomly overlap or leave parts 
            // of the signal unanalyzed 
            M = floor(grand(1, m, "unf", 1, len-N+2));
            M = gsort(M,"g","i");
        case "uniform"
            // Using quasi-uniform sampling, covering
            // the full signal length with possibly
            // overlapping windows.
            // Note: the offset between windows is 
            // constant +/-1  
            M = [1, round((len-N+1)/(m-1)*(1:m-1))];
        case "poison"
            // The first M version contains a sequence of
            // m-1 instants with a Poisson distribution
            M =  cumsum(-log(1 - rand(1, m-1)));
            // The Posisson sequence is normalized so that
            // the last instant is the starting sample of
            // the last window, covering up to the last
            // available sample
            M = round(M/max(M)*(len-N+1));
            // 1 is added at the beginning
            M =  [1, M];
        end
        // Windows to be analyzed are concatenated 
        // to form a new signal to be reshaped to
        // have m columns
        xx = [];
        for k = 1:m
            xx = [xx; x(M(k):M(k)+N-1, :)];  
        end
        x = xx;               
        clear("xx")
        // Reshape back the array to recover the original
        // structure in the higher dimensions. The result
        // is an array of columns possibly longer than the
        // original ones where the overlapped windows have
        // been unoverlapped    
        x = matrix(x, [m*N, siz(2:$)]);
        // Reshape the array splitting the columns into
        // m analysis windows of length N each.
        x = matrix(x, [N, m, siz(2:$)]);
    end
    
    // L is a positive integer
    if type(L)==1 
        // Number of FFT windows, including possibly 
        // a last incomplete window,
        m = ceil((len - N)/L) + 1;    
        // Number of complete windows. 
        q = floor((len - N)/L) + 1;
        if q<m
            // Zero-padding of the last, incomplete window
            // of each column so that columns have m 
            // complete windows 
            x = resize_matrix(x, [q*L+N, siz(2:$)]);
            // Length of incomplete window
            H = len - q*L;
        else
            H = 0;
        end
        // New signal size
        siz1 = size(x);
        // Length of each column containing signal
        len1 = siz1(1);
        // Detect if offset L is different from N
        // If this is the case, rebuild x by concatenating
        // successively the originally overlapping FFT 
        // windows so that they do not overlap any longer.
        if  L<>N
            // Vector containing starting samples of windows
            M = (0:m-1)*L + 1;
            // Reshape array so that the remaining dimensions
            // colapse into a single one (columns)
            x = matrix(x, [len1, prod(siz1(2:$))]);
            // For each column select successive (possibly
            // overlapping) windows and concatenate them  
            // into a single, longer column containing m 
            // non-overlapping windows
            xx = [];
            for k = 1:m
                xx = [xx; x(M(k):M(k)+N-1, :)];    
            end
            x = xx;
            clear("xx")
            // Reshape back the array to recover the original
            // structure in the higher dimensions. The result
            // is an array of columns (longer than the original
            // ones if L < N and shorter if L > N), where the
            // overlapped windows have been unoverlapped (or
            // the non-overlapped ones have been juxtaposed)    
            x = matrix(x, [m*N, siz1(2:$)]);
        end
        // New signal size
        siz2 = size(x);    
        // New number of FFT windows
        m = siz2(1)/N;
        // Reshape the array to split each column into
        // m columns of N samples so that the spectrum
        // can be computed and averaged columnwise. The
        // remaining dimensions are pushed to the right. 
        x = matrix(x, [N, m, siz2(2:$)]);          
    end

    // Prepare for the computation of the FFT spectrum
    // First we create an array without the last incomplete
    // windows (q instead of m), if they indeed exist.
    x1 = resize_matrix(x, [N, q, size(x)(3:$)]); 
    
    // Then, if there is indeed a last incomplete window  
    // (q < m), we create an array containing only that 
    // window for all dimensions
    if q<m
        // Rearrange x colapsing all dimensions from the 
        // third on in a single dimension and extract the
        // last column (column m)
        x2 = x(:,m,:);
        // Return to the original structure for the 
        // dimensions 3 to last
        x2 = matrix(x2, [N,1,size(x)(3:$)]);
    else
        // If there is no incomplete window
        x2 = [];
    end

    // Windowing function
    // In case of unrecognized window, Blackman is selected
    // by default 
    if w=="black"
        w = blackman(N)(:);
    elseif w=="hann"
        w = hann(N)(:);
    elseif w=="box" | w==1
        w = ones(N,1);
    elseif type(w)==1 & length(w)==N;
        w = w(:);
    else
        w = blackman(N)(:);
    end
  
    // Obtain the spectrum of the complete windows
    // Use the regular default blackman window
    [X1, fo] = octaves(x1, D, Fs, N, w);
        
    // If necessary, a second version of the window is
    // computed, covering the length H of the last, 
    // incomplete window, and the whole procedure is
    // repeated for this window. 
    if H<>0
        if w=="black"
            w$ = blackman(H)(:);
        elseif w=="hann"
            w$ = hann(H)(:);
        elseif w=="box" | w==1
            w$ = ones(H,1);
        elseif type(w)==1 & length(w)==N;
            // This window has the same shape as w but
            // it has H components
            w$ = interp1(1:N, w, (1:H)*N/H)(:);
        else
            w$ = blackman(H)(:);        
        end
        // Zero-pad last window function to get uniform
        // length
        w$ = resize_matrix(w$,[N,1]);
        // Obtain the spectrum of the incomplete window
        // Use a Blackman window extending to the real
        // duration of the signal (less than a frame)

        X2 = octaves(x2, D, Fs, N, w$);
    else
        // If there is no incomplete window then X2 
        // has the required structure and size but 
        // filled with zeros
        X2 = zeros(size(X1)(1), prod(size(x)(3:$)));
        X2 = matrix(X2,[size(X1)(1),1,size(x)(3:$)]);
    end

    // Computation of average, energy average, maximum and minimum
    // of band spectra
    // The average is computed along the different columns of each
    // row. To average along a row means that the values of the 
    // octave band corresponding to that row for different frames 
    // are averaged. This process yields a single column vector 
    // containing the bands for each signal.  
    
    // Energy mean of bands
    // Spectra of complete frames are squared and weighted
    // by their length N 
    y2N = X1.^2*N;
    // Spectra of last incomplete frames are squared and
    // weighted by their active length H 
    y2H = X2.^2*H;
    // Complete and incomplete frames are concatenated, 
    // summed along columns (each row corresponds to a 
    // band, all the columns in the same row are summed)
    // and divided by the total active length. Finally
    // the square root is computed
    y2 = sqrt(sum([y2N, y2H], 2) / (N*q + H));
    
    // Linear mean of bands
    // Spectra of complete frames are weighted by their 
    // length N 
    yN = X1*N;
    // Spectra of last incomplete frames are 
    // weighted by their active length H 
    yH = X2*H;  
    // Complete and incomplete frames are concatenated, 
    // summed along columns and divided by the total 
    // active length. 
    y = sum([yN, yH], 2) / (N*q + H);
      
    // Maximum and minimum
    // Spectra of complete frames and incomplete frames
    // are concatenated and the maximum and minimum for
    // each row are obtained along columns.
    ymax = max([X1, X2], "c");
    if q<m
        ymin = min([X1, X2], "c");  
    else
        // This is to avoid getting a spurious 0 due
        // to having set X2 to a null value
        ymin = min(X1, "c");
    end
    
    // Remove singletone column dimension  
    y2 = matrix(y2, [size(y2)(1),size(y2)(3:$)]);
    y = matrix(y, [size(y)(1),size(y)(3:$)]);  
    ymax = matrix(ymax, [size(ymax)(1),size(ymax)(3:$)]);  
    ymin = matrix(ymin, [size(ymin)(1),size(ymin)(3:$)]);  
    
    
    // Calibration correction taking into account that
    // cal is the level that corresponds to a calibration
    // signal referenced to maximum (+/- 1 = 0 dB)
    calib = dB - cal; 
    
    // Convert to levels and apply calibration
    // Note: Decimal logarithm currently does
    // not support hyperarrays
    Lx  =  20*log10(y) + calib;
    Lx2 =  20*log10(y2) + calib;
    Lmax = 20*log10(ymax) + calib;
    Lmin = 20*log10(ymin) + calib;

    if type(graph)==15
        // graph is included as an input argument and 
        // is a list whose first element is a string
        // and whose remaining elements are indexes for
        // the upper dimensions of x
        // Replicate list
        G = graph();
        // Delete the first argument and keep only 
        // the index elements 
        G(1) = null();
        // Replace graph by its string first element 
        graph = graph(1);
    else
        // Otherwise create a list with the first index
        // of each dimension 
        G = list();
        for n=1:length(size(Lx))-1
            G(n)=1;
        end; 
    end
    
    // Open new figure
    if graph=="fig"
        scf();
        clf();
    end

    // Plot on current figure
    if graph=="fig" | graph=="cuf"
        // Convert indices of signal to plot to string
        s = string(G(1));
        for n=2:length(G)
            s = s + ", " + string(G(n)); 
        end        
        // Give a name to the type of band
        if isvector(D)
            band = "custom band ";
        else
            select D
            case 1
                band = "octave band ";
            case 2
                band = "half octave band "
            case 3
                band = "one third octave band "
            case 4
                band = "one quarter octave band "
            case 6 
                band = "one sixth octave band "
            case 12
                band = "one twelfth octave band "
            else
                band = "band "
            end
        end
        // Title
        tit = "Mean " + band + "spectra of signal " + "(" + s + ")";
        
        // Prepare data for graph of the octave band
        // minimum, linearly averaged, energy averaged 
        // and maximum spectra. Only the first signal
        // will be displayed 
        LLmax = Lmax(:, G(:));
        LLx = Lx(:, G(:));
        LLx2 = Lx2(:, G(:)); 
        LLmin = Lmin(:, G(:)); 
        
        // Plot
        plot("ln", fo, [LLmax, LLx, LLx2, LLmin],"o-")
        legend("L_max", "L_mean", "L2_mean", "L_min");
        title(tit, "font_size", 2);
        xlabel("f, Hz", "font_size", 2);
        ylabel("L, dB", "font_size", 2);
        gca().children(1).background = -2;
        gca().grid=[1 1 1];
    end

endfunction

function w = hann(N)
    // Periodic Hann window
    // Usage:
    //       w = hann(N)
    // Where 
    //       N: Length of Hann window
    //       w: Vector containing the periodic Hann window
    // ---------------------------------- 
    w = 0.5 - 0.5*cos(2*%pi*(0:N-1)/N);
endfunction

function w = blackman(N)
    // Periodic Blackman window
    // Usage:
    //       w = blackman(N)
    // Where 
    //       N: Length of Blackman window
    //       w: Vector containing the periodic Blackman window
    // ---------------------------------- 
    w = 0.42 - 0.5*cos(2*%pi*(0:N-1)/N) + 0.08*cos(4*%pi*(0:N-1)/N);
endfunction
