function [y2, y, ymin, ymax, fo] = fftmean(x, Fs, N, w, L)
    // FFT spectrum analysis     
    // This procedure computes the mean FFT spectrum of a signal. 
    // Usually x will have been read from an uncompressed audio
    // file in WAV format using the function wavread, so maximum 
    // signal will be +/- 1. The spectrum is presented as a real
    // unilateral spectrum where each spectral line is the RMS
    // value of the corresponding sine wave. The first line is 
    // the mean value (DC component) of the signal.
    //
    // Usage:
    //      y2 = fftmean(x, Fs, N, w, L)
    //      [y2, y, ymin, ymax, fo] c= fftmean(x, Fs, N, w, L)
    //
    // where
    //      x:    Signal column or row vector or array of columns
    //      Fs:   Sampling rate in Hz
    //      N:    FFT window length
    //      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 (equivalent to no window)
    //            "hann":  Periodic Hann window
    //            "black": Periodic Blackman window
    //      L:    Offset in samples between successive FFT windows
    //            or a list whose first element is a string that
    //            indicates the type of window sampling: "resample"
    //            for bootstrap resampling, "uniform" for uniform
    //            distribution of analysis windows, and "poisson" 
    //            for ensuring that the choice of windows is 
    //            controlled by a Poisson's process; and whose
    //            second element indicates the number of window
    //            spectra to be averaged.
    //      y2:   Column vector or array of column vectors with 
    //            energy-averaged FFT spectrum of the input signal
    //      y:    Column vector or array of column vectors with
    //            linearly averaged FFT spectrum of the input signal
    //      ymax: Column vector or array of column vectors with
    //            the maximum FFT detected for each spectral line
    //      ymin: Column vector or array of column vectors with
    //            the minimum FFT detected for each spectral line
    //      fo:   Column vector containing the frequencies of the
    //            spectral lines
    //
    // Arguments Fs, N, w an L are optional; their default values 
    // are, respectively, 44100, 4096, "box" and N.
    //
    // Application notes
    // 1) The result is expressed as the RMS values of signal 
    //    components at frequencies 0, Fs/N, 2*Fs/N, ..., 
    //    (N/2-1)*Fs/N
    // 2) Argument L allows a variety of methods to pick the
    //    spectral windows to be averaged. It allows an offset
    //    in samples, bootstrap resampling (which means that 
    //    a requested number m of analysis windows are randomly 
    //    chosen from the signal) or uniform sampling, which
    //    means that m windows are chosen uniformly distributed 
    //    along the signal    
    // 3) The function accepts any multidimensional array as
    //    long as the signals are distributed as columns. The
    //    average process is carried out for each signal, i.e.,
    //    along columns. If the signals are rows, instead, but
    //    the first (column) dimension is 1, the function 
    //    permutes the first two dimensions.
    // 4) In the case of pure tones, the spectral peak is not 
    //    exactly the RMS due to spectral leakage. A better
    //    estimate can be obtained with an energy superposition
    //    of the peak and a few of its side lines, such as:
    //    [m,i] = max(y2(:,1,2,1))
    //    yrms = sqrt(sum(y2(i-1:i+1,1,2,1).^2))
    // 5) For pure tones whose period is a submultiple of N,
    //    if the signal length is a multiple of N the spectral
    //    leakage is vanishingly small but if there is an 
    //    incomplete final window then the spectral leakage
    //    rises considerably    
    // 6) It is possible to estimate the power spectral density 
    //    as y2.^2*N/Fs
    // 7) If the sound pressure level Lcal of a sine wave of 
    //    amplitude 1 is known then the spectral level can be 
    //    obtained as 
    //    lp = Lcal + 10*log10(2*N/Fs)
    //    where Fs is in Hz and lp isin dB referred to 1 Pa^2/Hz.
    //
    // -------------------------------------------------------------
    // Author:   Federico Miyara
    // Date:     2008-04-24
    // Revision: 2008-08-02
    //           2010-04-14
    //           2011-03-20
    //           2011-10-15
    //           2012-12-24
    //           2016-01-03
    //           2020-01-22
    //           2021-01-02
    //           2021-01-05
    //           2021-01-24
    //           2021-02-09
    //           2024-09-15

    // INPUT DATA --------------------------------------------------
    
    // Number of input (right-hand) arguments
    rhs = argn(2);
    
    // Default values of missing arguments 
    if rhs<2
        // Usual sampling rate in CD quality
        Fs = 44100;
    end
    if rhs<3
        N = 4096;
    end    
    if rhs<4
        w = "box";
    end
    if rhs<5
        L = N;
    end


    // ------------------------------------------------
    // TEST-CASE DATA - DELETE
    // ------------------------------------------------
    if 1==2
        Fs = 44100;
        w = "black";
        caso = 1
        select caso
        case 1
            exec pinknoise.sci;
            N = 4096;
            L = 2048;
            //L = list("resample", 40)
            //len = 100000;        
            len = 98304;
            Npink = 8192;
            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, Npink);
                end
                x(1:len, siz(2), h) = sin(2*%pi*1000*[0:len-1]/Fs);
            end
        case 2
            N = 8
            L = 8
            x = 1:240;
            x = matrix(x, 20, 4, 3);
        end

    end
    //-------------------------------------------------


    // 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;
        // 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 "poisson"
            // 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;
        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((k-1)*L+1:(k-1)*L+N, :)];
                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
    


    // 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 frame
        x2 = [];
    end

    // Windowing function
    // In case of unrecognized window, Blackman is selected
    // by default 
    if w=="black"
        exec blackman.sci;
        w = blackman(N)(:);
    elseif w=="hann"
        exec hann.sci;
        w = hann(N)(:);
    elseif w=="box" | w==1
        w = ones(N,1);
    elseif type(w)==1 & length(w)==N;
        w = w(:);
    else
        exec blackman.sci;
        w = blackman(N)(:);
    end    
    // Calculate RMS value of windows in order to
    // compensate the energy reduction caused by
    // windowing
    wrms = sqrt(mean(w.^2));
    // Replicate window to all remaining dimensions
    wa = repmat(w, [1, q, size(x1)(3:$)]);

    // Obtain the spectrum of the complete frames
    // applying the specified windowing function
    X1 = abs(fft(x1.*wa, -1, 1)) / wrms;
    // Keep only the first N/2 spectral lines, which
    // correspond to the single-sided spectrum
    X1 = resize_matrix(X1, [N/2, size(X1)(2:$)]);
    // Size of X1
    X1siz = size(X1);
    // Apply factor 1/N to the first component of
    // the columns representing complete FFT windows 
    // (across all the remaining dimensions) to get
    // the mean value
    X1(1,:) = X1(1,:)/N;
    // Apply factor 2/N/sqrt(2) to the remaining
    // components to get the RMS value of the sine
    // represented by such components
    X1(2:N/2,:) = X1(2:N/2,:)*2/N/sqrt(2); 
            
    // 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 isdef("H")
        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;
            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]);
        // Calculate RMS value of the last windows
        w$rms = sqrt(mean(w$.^2));
        // Replicate last window to all remaining dimensions
        w$a = repmat(w$, [1, 1, size(x2)(3:$)]);
            
        // Obtain the spectrum of the incomplete window
        X2 = abs(fft(x2.*w$a, -1, 1)) / w$rms;
        // Keep only the first N/2 spectral lines
        X2 = resize_matrix(X2, [N/2, size(X2)(2:$)]);
        // Size of X2
        X2siz = size(X2);    
        // Apply factor 1/N to the first component
        X2(1,:) = X2(1,:)/N;
        // Apply factor 2/N/sqrt(2) to the remaining
        // components
        X2(2:N/2,:) = X2(2:N/2,:)*2/N/sqrt(2);
    else
        // To avoid errors and simplify the code, dummy 
        // (harmless) values are provided for the
        // spectrum of the nonexistent incomplete window 
        // and its length
        X2 = [];
        H = 0;         
    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 FFF component corresponding to that 
    // row for different FFT windows are averaged. This 
    // process yields a single column vector containing
    // the FFT values for each signal.  
    
    // Energy mean of FFT spectral lines
    // 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 across columns (each row corresponds to a 
    // spectral line, 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 across 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");
    ymin = min([X1, X2], "c");
    
    // Dimensional adjustmnt to remove singletone rows
    y2 = matrix(y2, [N/2, size(y2)(3:$)]);
    y = matrix(y, [N/2, size(y)(3:$)]);
    ymax = matrix(ymax, [N/2, size(ymax)(3:$)]);
    ymin = matrix(ymin, [N/2, size(ymin)(3:$)]);


    // FFT frequency vector
    fo = [0:N/2-1]'*Fs/N;


    // Remove conditional to add calibration capability. In such
    // case caldB and calwav must be added as input arguments and
    // Lx, Lx2, Lmax, Lmin as output arguments
    if 1==2 
        // Calibration correction taking into account that
        // calwav is the level that corresponds to a calibration
        // signal referenced to maximum (+/- 1 = 0 dB)
        cal = caldB - calwav; 

        // Convert to levels and apply calibration
        // Note: Decimal logarithm currently does
        // not support hyperarrays
        Lx  =  20*log(y)/log(10) + cal;
        Lx2 =  20*log(y2)/log(10) + cal;
        Lmax = 20*log(ymax)/log(10) + cal;
        Lmin = 20*log(ymin)/log(10) + cal;
    end
    
    // Remove conditional to add graphic capability. In such
    // case graph must be added as input variable. When set
    // to "fig" a new figure is created and all of the output
    // variables are plotted. When set to "cuf" the current
    // figure is used to plot the variables
    if 1==2
    
        // Open new figure
        if graph=='fig'
            figure
            gcf().background=8
        end

        // Plot on current figure
        if graph=='fig' | graph=='cuf'
            // Prepares data for graph of the octave band
            // minimum, linearly averaged, energy averaged 
            // and maximum spectra. Only the first signal
            // will be displayed 
            LLmax = Lmax(:,1,1,1);
            LLx = Lx(:,1,1,1);
            LLx2 = Lx2(:,1,1,1); 
            LLmin = Lmin(:,1,1,1); 

            // Plot
            plot('ln', fo, [LLmax, LLx, LLx2, LLmin])
            legend('Loct_max', 'Loct_mean', 'Loct_mean2', 'Loct_min');
            gca().background = 8;
            gca().grid=[1 1 1];
        end
    end

endfunction

