function [y, fo] = octaves(x, D, Fs, N, w)
    // Octave or fractional octave band spectrum 
    // This function computes the octave or fractional octave 
    // band spectrum of vector x corresponding to a signal
    // sampled at a sample rate Fs. It yields a column vector  
    // y whose components are the effective values in the 
    // octave or fractional octave bands centred at standard
    // frequencies between 20 Hz and 20000 Hz. For instance, 
    // for octave band spectrum these frequencies are 31,5 Hz, 
    // 63 Hz,  125 Hz,  250 Hz,  500 Hz,  1000 Hz,  2000 Hz, 
    // 4000 Hz, 8000 Hz and 16000 Hz. It also allows a custom
    // band set. A windowing function is applied to reduce
    // spectral leakage.
    // 
    // Usage:
    //        [y, fo] = octaves(x, D, Fs, N, w)
    //
    // where
    //        x:  signal vector or array
    //        D:  Denominator of fractional octave band. 
    //            For instance, D = 1 for octave band 
    //            analysis and D = 3 for one third octave
    //            band analysis. D may also be a vector of 
    //            frequencies so that any two consecutive 
    //            frequencies defines a custom band. In any 
    //            case the analysis limits to the highest
    //            band with content below Fs/2.  
    //        Fs: sampling rate in Hz
    //        N:  FFT size
    //        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
    //        y:  Column vector or array containing the RMS 
    //            spectrum of the signal
    //        fo: Column vector containing band center
    //            frequencies. In the case of a custom band
    //            set they are computed as the geometric 
    //            mean of each band's bounds.     
    //
    // For single signals, this function accepts column or 
    // row vectors, but if x is a matrix (or N-D array) 
    // containing several signals, the analysis is carried 
    // out by columns, yielding a matrix (or N-D array), y,
    // whose columns have as many components as frequency 
    // bands according to D (for instance, 10 for octave 
    // band analysis) and are distributed according to the
    // original matrix or N-D array. Each column is the  
    // octave band spectrum of the respective columns of x 
    // (or row if columns are singleton). This function 
    // uses by default a Blackman window to minimize
    // spectral leakage, but it also allows to specify 
    // other windows, including custom windows. The energy 
    // reduction of the signal caused by whatever window
    // is used is fixed by dividing by the window RMS value. 
    // Fs, D, N and w are optional. Their default values 
    // are 44100, 1 (octave), 8192 and "black". If 
    // Fs < 44100 Hz then the number of nonzero bands may
    // be reduced to avoid exceeding the Nyquist frequency 
    // Fs/2. N should be a power of 2 for dramatic speed 
    // increase. However, if it is too small then precision
    // at the lower bands is impaired.
    //
    // NOTES: 
    // 1) If signal vector is shorter than N samples it is
    //    zero-padded. If it is longer, it is truncated 
    //    so that only the first N samples are used. For
    //    longer signals, use octavemean().
    // 2) For signals of length len shorter than N samples
    //    the octave bands are referred to the real signal
    //    length. To refer it to N instead, just multiply
    //    by sqrt(len/N).
    // 3) While octaves() yields the result as an RMS 
    //    spectrum, octavemean() provides a calibrated
    //    spectrum in dB. o convert the octaves() spectrum
    //    to the octavemean() format, use
    //    Loct = 20*log10(y) + dB - cal
    //    where dB is the sound pressure level of the 
    //    acoustic calibration signal and cal is the level
    //    of the corresponding digital signal, referenced 
    //    to an undistorted sine wave of maximum amplitude
    //    in dBFS.    
    //
    // ---------------------------------------------------------
    // Author: Federico Miyara
    // Date:   2005-10-09
    //         2009-09-06
    //         2010-06-08
    //         2011-01-31
    //         2011-02-06
    //         2020-01-18
    //         2024-10-20
    
    // Input argument check
    rhs = argn(2);

    if rhs<5
        w = "black"; 
    end
    if rhs==3
        // Reasonable length for analysis window.
        // FFT size is selected to increase resolution
        // when the spectral bands become narrower.
        N = 2^nextpow2(8192*D);        
    end
    if rhs==2
        // Octave band analysis 
        D = 1;
        N = 8192;
    end
    if rhs==1
        // Most common sampling rate
        Fs = 44100;
        D = 1;
        N = 8192;
    end


    if 1==2
        // --------------------------------
        // TEST-CASE DATA - DELETE
        // --------------------------------
        exec pinknoise.sci;
        w = "black";               // Window type
        N = 8192;                  // Window size  
        len = 4000;                // Signal length 
        Fs = 44100;                // Sample rate
        D = 1;                     // Band type 
        x = zeros(1,len,3);
        siz = size(x)
        for i=1:siz(1)
            for h=1:siz(3)
                x(i,1:len,h) = pinknoise((len-1)/Fs, Fs, N);
            end
        end
        //--------------------------------
    end


    // Size of argument x (may be multidimensional)
    // The first two dimensions are most important since
    // the spectrum will be calculated along the first 
    // non-singleton of them.
    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);
    
    // Actual length of signals
    len = siz(1);

    // Band retrieval
    if isvector(D)
        // D is a custom set of frequencies of band
        // boundaries. 
        // Center frequencies are the geometric mean
        // of each pair of consecutive frequencies
        fo = sqrt(D(1:$-1).*D(2:$));
        // Index to the maximum frequency less than
        // Fs/2. This is the lower limit of the maximum
        // band
        Nsup = max(find(D<Fs/2));
        bands = Nsup - 1;
        // Lower index for band h
        k1 = ceil(D(1:bands)*N/Fs);
        // Upper index for band h 
        k2 = floor(D(2:bands+1)*N/Fs);
        k2($) = min(N/2-1,k2($));
    else
        // D is indeed the denominator of an octave
        // fraction
        // Fractional octave band center frequencies
        fo = fstandard(D)';
        // Index to frequency 1000 Hz
        No = find(fo==1000);
        // Number of bands compatible with Nyquist frequency
        bands = floor(No + 0.5 + D*log10(Fs/2/1000)/log10(2));
        bands = min([bands, length(fo)]);
        // Lower index for band h
        k1 = ceil(1000 * 2^(((1:bands) - No - 0.5)/D)*N/Fs);  
        // Upper index for band h 
        k2 = floor(1000 * 2^(((1:bands) - No + 0.5)/D)*N/Fs);
        k2($) = min(N/2-1,k2($));
    end

    // Window retrieval
    if type(w)==1 & length(w)==N 
        // A valid custom window, just ensure it 
        // is a column vector and go ahead
        w = w(:);
    elseif type(w)==1 & length(w)==len
        // A custom window whose length is equal
        // to the signal length, which is less
        // than the FFT window. Zero-pad to 
        // complete FFT length
        w = resize_matrix(w(:), [N,1]);
    elseif w=="box" | w==1 
        // Rectangular window (no window, "do no harm")
        w = resize_matrix(repmat(1, [len,1]), [N,1]);
    elseif w=="hann"
        // Periodic Hann window as column vector
        // extended to the actual signal length 
        // and zero-padded to FFT window size 
        w = resize_matrix(hann(len)', [N,1]);
    elseif w=="black"
        // Periodic Blackman window as column vector
        // extended to the actual signal length 
        // and zero-padded to FFT window size 
        w = resize_matrix(blackman(len)', [N,1]);
    else
        // In case of any unrecognized window type
        // select the default Blackman window
        // zero-padded to FFT window size 
        w = resize_matrix(blackman(len)', [N,1]);
    end
   
    // Zero-pad columns if they have less than N
    // components; or trim out components in excess
    // if columns have more than N components
    x = resize_matrix(x, [N, siz(2:$)]);
    
    // Replicate window along columns (first dimension)
    ww = repmat(w, [1, siz(2:$)]);
   
    // Apply window
    xw = x.*ww;
    
    // Apply FFT along columns and adjust components
    // for real spectrum
    e = 2/N * abs(fft(xw, -1, 1));
    e(1, :) = e(1, :)/2;
    
    // Retain the N/2 spectral lines from each column
    e = resize_matrix(e, [N/2, siz(2:$)]);

    // Computation of the octave band spectrum of each
    // column from the FFT line spectrum already 
    // obtained. Then each line is classified according
    // to the band to which it belongs and finally their
    // mean square values within each band are added. 
    //
    // Notes: 1) Exact values of centre frequencies of
    // the octave bands, i.e., 1000*2^(h-6), are used
    // instead of the standard ones. Due to the high
    // slope at the bands' limits, this may cause that
    // a signal very close to the cutoff of a band be
    // considered to be at a different band. This effect
    // may be compensated by the spectral widening due
    // to the Blackman window
    // 2) The number of bands compatible with the sample
    // rate and audibility (f < 20000 Hz) is computed
    // 3) Audible bands beyond the Nyquist limit are set
    // to 0 
    
    // Allocate space for bands 
    sizoct = [bands, siz(2:$)];
    y = repmat(0, sizoct);
    
    // Reshape and colapse all dimensions other than
    //  columnsinto a single dimension
    ee = matrix(e, [N/2, prod(siz(2:$))]);

    // Spectral lines within band h are added on an
    // energy basis
    for h = 1:bands
        y(h,:) = sqrt(sum(e(k1(h):k2(h), :).^2, 1));
    end
    
    // Bands beyond Nyquist rate are set to 0
    for h = (bands+1):size(y)(1)  
        y(h,:) = 0;
    end
    
    // Recover original non-column dimensions 
    y = matrix(y, [bands, siz(2:$)]);
    
    // Correction for window attenuation
    y = y / sqrt(mean(w.^2));
    
    // Conversion to RMS (prior to conversion, each
    // spectral line is the amplitude of a sine wave)
    y = y / sqrt(2); 

endfunction

// -------------------------------------------------------------
// LOCAL FUNCTION
// -------------------------------------------------------------

function y = fstandard(n)
    // This function yields a column vector with 
    // center frequencies of standard fractional
    // octave bands according to ISO 3.
    //
    // Usage:
    //       y = fstandard(n)
    //
    // If n = 1, 2, 3, 4, 6, 12 a vector with the
    // center frequencies of bands of 1/N of an
    // octave between 20 Hz and 20000 Hz. When n 
    // is 3, 6 or 12 the result ar the Renard 
    // series R10, R20 and R40. If n is different
    // from the standard values, then the exact
    // values between 20 Hz and 20000 Hz are 
    // chosen in such a way that 1000 Hz belong
    // to the series.
    //
    // NOTICE: From the 2024 version arguments 2 
    //         and 3 have been reversed, so that 
    //         octaves() and octavemean() have 
    //         these arguments in the same order.
    //         All functions and scripts of this
    //         package using octaves() have been
    //         adapted to these change so that   
    //         they keep working.
    //  
    // ------------------------
    // Federico Miyara
    // 2007-06-01
    // 2010-12-11
    // 2011-02-06
    // 2017-08-08
    // 2024-10-21

    // The vector is filled out with standard
    // tables of cente frequencies of fractional
    // octave bands. If non-standard, the 
    // frequencies are calculated
    switch n
    case 1
        y = [31.5 63 125 250 500 1000 2000 4000 8000 16000];
    case 2
        y = [22.4 31.5 45 63 90 125 180 250 355 500 710 ...
        1000 1400 2000 2800 4000 5600 8000 11200 16000];
    case 3
        y = [20 25 31.5 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];
    case 4
        y = [22.4 26.5 31.5 37.5 45 53 63 75 90 106 125 ...
        150 180 212 250 300 355  425 500 600 710 850 1000 ...
        1180 1400 1700 2000 2360 2800 3350 4000 4750 5600 ...
        6700 8000 9500 11200 13200 16000 19000];
    case 6
        y = [20 22.4 25 28 31.5 35.5 40 45 50 56 63 71 ...
        80 90 100 112 125 140 160 180 200 224 250 280 ... 
        315 355 400 450 500 560 630 710 800 900 1000 ... 
        1120 1250 1400 1600 1800 2000 2240 2500 2800 3150 ...
        3550 4000 4500 5000 5600 6300 7100 8000 9000 ... 
        10000 11200 12500 14000 16000 18000 20000]; 
    case 12
        y = [20 21.2 22.4 23.6 25 26.5 28 30 31.5 33.5 ...
        35.5 37.5 40 42.5 45 47.5 50 53 56 60 63 67 ...
        71 75 80 85 90 95 100 106 112 118 125 132 140 ...
        150 160 170 180 190 200 212 224 236 250 265 280 ... 
        300 315 335 355 375 400 425 450 475 500 530 560 ...
        600 630 670 710 750 800 850 900 950 1000 1060 ... 
        1120 1180 1250 1320 1400 1500 1600 1700 1800 1900 ...
        2000 2120 2240 2360 2500 2650 2800 3000 3150 3350 ...
        3550 3750 4000 4250 4500 4750 5000 5300 5600 6000 ... 
        6300 6700 7100 7500 8000 8500 9000 9500 10000 ...
        10600 11200 11800 12500 13200 14000 15000 16000 ...
        17000 18000 19000 20000];
    otherwise
        k1 = ceil(n*log2(1/50));
        k2 = floor(n*log2(20));
        y = 1000 * 2.^([k1:k2]/n);
    end
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

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
