function [f, Xp, phi, BW] = spectralpeaks(x, N, Mo, TH, Fs)
    // Spectral peaks detector
    // This function detects the main spectral peaks in an 
    // input signal.
    //
    // Usage:
    //        [f, Xp, phi, BW] = spectralpeaks(x, N, Mo, TH, Fs)
    // where
    //        x:   signal vector
    //        N:   size of analysis window
    //        Mo:  desired number of peaks
    //        TH:  threshold peak amplitude     
    //        Fs:  sampling rate in Hz
    //        f:   vector containing the spectral peaks in Hz
    //        Xp:  vector containing the amplitudes of peaks  
    //        phi: vector containing the phase of peaks
    //        BW:  vector containing the bandwidth of peaks
    //
    // If length of x does not match N, then it is zero-padded
    // or truncated to fit analysis window. If x is a column
    // vector output vectors are column vectors; otherwise
    // they are row vectors. Input arguments N, Mo, TH, Fs are 
    // optional. Their default values are 2048, 10, 0.01 and 
    // 44100.  Output arguments Xp, phi, BW are optional.
    //
    // Peaks are estimated by finding first the local maxima of
    // the FFT spectrum y(n), then interpolating parabolically 
    // y(n-1), y(n), y(n+1) and finally computing the frequency
    // of the vertex of the parabola. The magnitude of the peak 
    // is obtained by Parseval energy addition of side bands up
    // to the order 2. 
    // Threshold stipulates the minimum value of a peak to be
    // accepted as valid. Threshold overrides number of peaks. 
    // Use TH = 0 to control the selection of peaks by number
    // exclusively and Mo = N/2 to control it by threshold.
    // If the number of peaks exceeding TH is smaller than Mo, 
    // then output vectors are zero-padded to length Mo. To 
    // save memory, this feature works only for Mo < 100.
    // Phase is approximated by parabolic interpolation of the
    // angle of y(n-1), y(n) and y(n+1), corrected for picket
    // fence effects. Due to rapid phase change when input
    // frequency varies from a spectral line to the next one,
    // the correction is quite sensitive to minute errors in
    // the peak estimation, so the approximation is not very
    // good.
    // Bandwidth is retrieved by successively broadening a
    // band around a peak up to the point when the spectrum 
    // falls 3 dB from the peak value. Bandwidth may prove 
    // meaningless in the case of noisy or broadband signals. 
    //
    // ------------------------------
    // Author: Federico Miyara
    // Date:   2010-06-22
    //         2010-07-05
    //         2020-04-15
    

    rhs = argn(2)

    if rhs<4 
        // Sampling rate
        Fs = 44100;
    end
    if rhs<3 
        // Maximum desired number of peaks
        Mo = 10;
    end
    if rhs<2
        // Size of analysis window
        N = 2048;
    end
    
    
    // FOR TESTING -- DELETE
    if 1==2
        clear
        exec localmax.sci;
        clc
        Mo = 6
        TH = 0.1
        Fs = 44100
        N = 4096
        t = [0:N-1]/Fs;
        Fo = 800
        // Polyharmonic signal
        x = 1.0*cos(2*%pi*Fo*t) + ...
            0.5*cos(2*%pi*2*Fo*t + %pi/3) + ...
            0.2*cos(2*%pi*3*Fo*t + %pi/4) + ...
            1.0*cos(2*%pi*4*Fo*t) + ...
            1.5*cos(2*%pi*5*Fo*t + %pi/2) + ...
            1.2*cos(2*%pi*6*Fo*t + %pi/5) + ...
            0.1*cos(2*%pi*7*Fo*t) + ...
            0.05*cos(2*%pi*8*Fo*t + %pi/3);
    end
    
    
    // Ensure x is a column vector
    if size(x,1)==1
        x = x(:);
        wasrow = 1;
    else 
        wasrow = 0;
    end

    // Zero-pad or truncate signal to fit analysis window size
    x = resize_matrix(x,N,1);

    // Get the FFT spectrum with Hann window to improve precision. A
    // correction factor is included to account for the N/2 factor of 
    // the definition of DFT and for the sqrt(6)/4 factor due to the 
    // Hann window.
    // NOTE: The correction factor is obtained considering the RMS
    //       value of 0.5*(1 - cos(w0*t))*cos(w*t) where w0 is the
    //       angular frequency corresponding to the period of the
    //       window analysis. This is easily computed considering the
    //       decomposition into a carrier (w) and two lateral bands
    //       (w + w0 and w - w0). The carrier is 1/2 the amplitude 
    //       of the original signal and the lateral bands are 1/4 the
    //       amplitude, yielding sqrt(1/4 + 1/16 + 1/16) = sqrt(3)/4.
    //       This, multiplied by sqrt(2) yields the equivalent peak,
    //       i.e., sqrt(6)/4. 
    y = fft(x .* window("hn", N+1)(1:N)') * 2/N * 4/sqrt(6);

    // Remove DC and very low frequency components 
    y(1) = 0; 
    y(2) = 0;
    
    // NOTATION
    // k: spectral index; each value of k indexes a spectral FFT line
    // h: index to the spectral maxima

    // Find the local maxima
    // h contains the indexes where local maxima of the spectrum
    // occur; these indexes correspond to frequencies. ymax contains
    // the corresponding spectral values
//  [h, ymax] = localmax(abs(y(1:N/2)),2);   // Delete
    [k, ymax] = localmax(abs(y(1:N/2)),2);
    
    // Sort the local maxima in ascending order to allow 
    // discarding those smaller than the threshold. Afterwards
    // the spectral order will be recovered.
    // h contains the sorted indexes to the local maxima
    // NOTE: These indexes correspond to the order in which the 
    //       maxima must be picked to have them in ascending
    //       order
    [Ymax, h] = gsort(ymax, "g", "i");

    // Remove peaks under the threshold
    hTH = find(diff([0; Ymax>TH]));
    Ymax = Ymax(hTH:$);
    h = h(hTH:$);

    // Prevent errors due to low maxima count 
    // (if there are less maxima than Mo then one must content
    // with the existing ones) 
    M = min(Mo, length(Ymax));

    // Preserve the M largest maxima (M main peaks)
    Ymax = Ymax($-M+1:$);

    // h now contains the M indexes h that correspond to the M
    // largest maxima
    h = h($-M+1:$);

    // H contains the spectral indexes (bins) corresponding to
    // frequencies where the M largest maxima occur; they are
    // ordered in ascending order
    K = gsort(k(h), "g", "i");


    // Preallocate Ko, phi, BW (improves slightly computation
    // speed --very slightly, since the peak vector is never 
    // that huge!--)
    Hmax = zeros(1, M);       // To remove
    Ko = zeros(1, M);     // Contains fractional indexes to maxima
    Xp = zeros(1, M);     // Contains more accurate peak values
    phi = zeros(1, M);    // Contains phase value estimates
    BW = zeros(1, M);     // Contains peak bandwidth estimates

    // Main algorithm
    for n=1:M
        // Complex values of the spectrum around the n-th peak
        YY = [y(K(n)-1); y(K(n)); y(K(n)+1)];
        
        // Absolute values of the spectrum around the n-th peak
        BY = abs(YY);
        
        // Spectrum phases around the n-th peak
        Bphi = atan(imag(YY) ./ real(YY));
                
        // Fractional index of the estimated peak
        // This formula comes from solving the equation system 
        // for finding the coefficients of a parabola passing 
        // through the maximum and its neighbours and from 
        // them computing the peak
        Ko(n) = K(n) + (BY(3)-BY(1)) / (4*BY(2)-2*BY(3)-2*BY(1));
        
        // Amplitude fitting (for testing only)
        // The´parameters of the following expression are 
        // computed:
        // Y = Ypeak + CY*(k - Ko)^2  
        CY = (BY(2) - BY(1)) / (2*(K(n) - Ko(n)) - 1);
        Ypeak = BY(2) - CY*(K(n) - Ko(n))^2;
      
        // Phase fitting
        // To perform it we need phimax, Kophi and Cphi, so that
        // phi = phimax + Cphi*(k - Kophi)^2    
        Kophi = K(n) + ...
                (Bphi(3)-Bphi(1)) / (4*Bphi(2)-2*Bphi(3)-2*Bphi(1));       
        Cphi = (Bphi(2) - Bphi(1)) / (2*(K(n) - Kophi) - 1);
        phimax = Bphi(2) - Cphi*(K(n) - Kophi)^2;
        
        // The value of phi at the position of the maximum 
        // amplitude is computed using the parabolic model
        // with the parameters just derived 
        phi(n) = phimax + Cphi * (Ko(n) - Kophi)^2;

        // The following is an alternative method computing 
        // the coefficients of the parabola from the system
        // of equations for a parabola passing through three
        // points 
        // System matrix for solving the parabolic interpolation.
        // NOTE: While cubic interpolation with interp1 would be 
        //       an easier method, is about 10 times slower than
        //       the present algorithm
        // The interpolation is local and will be valid between 
        // H(n)-1 and H(n)+1
        A = [(K(n)-1)^2, K(n)-1, 1;...
        K(n)^2,    K(n),   1;...
        (K(n)+1)^2, K(n)+1, 1];
        
        // System solution for the spectrum amplitude (yields the
        // coefficients of the parabola in order of decreasing 
        // degree)
        coefY = inv(A)*BY;
     
        // Interpolated maximum of the parabola. It yields the
        // fractional index corresponding to the estimated
        // frequency
        // NOTE: As an index it is meaningless since fractional
        //       indexes have no value assigned, but it does 
        //       represent a frequency that cannot be represented
        //       by means of any of the N FFT lines
        Ko(n) = -coefY(2)/2/coefY(1);

        // System solution for the spectrum phase (yields the
        // coefficients of the interpolating parabola)
        coefphi = inv(A) * Bphi;

        // Parseval estimation of peak amplitude from neighbouring
        // spectral lines. 
        // NOTE: This is far best an estimation than the peak of 
        //       the interpolating parabola 
        Xp(n) = sqrt(sum(abs(y(K(n)-2:K(n)+2)).^2)); 

        // Estimation of the phase of partial n by parabolic
        // interpolation
        phi(n) = coefphi(1)*Ko(n).^2 + coefphi(2)*Ko(n) + coefphi(3);

        // Estimation of the bandwidth 
        // Starting values for the iterative process
        Y3dB = abs(Xp(n));
        n3dB1 = K(n);
        n3dB2 = K(n);
        // Start a process in which a band centered at peak n
        // will be gradually broadened until the spectral values
        // at its limits fall 3 dB
        while Y3dB > abs(Xp(n))/sqrt(2)
            // Widen the candidate bandwidth
            n3dB1 = n3dB1 - 1;
            n3dB2 = n3dB2 + 1;
            // Make sure the extremes of the window analysis 
            // have not been reached
            if (n3dB1>=1) * (n3dB2<=length(y))
                // Update the magnitude at the band limits
                Y3dB = max(abs(y(n3dB1)), abs(y(n3dB2)));
            else
                break
            end
        end
        BW(n) = (n3dB2 - n3dB1) * Fs/N;      
    end

    // Frequency estimation from fractional index
    f = (Ko-1)*Fs/N;

    // Phase must be corrected for picket fence effects

    // Fractional indexes corresponding to estimated frequencies
    nu = f*N/Fs;

    // Nearest integer index corresponding to peak of the
    // FFT spectrum
    n0 = round(nu);

    // Picket fence phase correction
    phi = phi - pi*(nu - n0);

    // If Mo is less than 100, zero-pad output vectors to a
    // length Mo
    if Mo < 100
        f =   [f,   zeros(1, Mo-M)];
        Xp =  [Xp,  zeros(1, Mo-M)];
        phi = [phi, zeros(1, Mo-M)];
        BW =  [BW,  zeros(1, Mo-M)];
    end

    // If input signal was a column vector then output data
    // are presented as column vectors
    if wasrow==0
        f = f(:);
        Xp = Xp(:);
        phi = phi(:);
        BW = BW(:);
    end

endfunction

// ---------------------------------------------------------
// PRIVATE FUNCTION
// ---------------------------------------------------------

function [i,y] = localmax(x)
    // This function finds local maxima within vector x
    //
    //    [i, y] = localmax(x) 
    // 
    // outputs a vector i containing the indexes where x 
    // reaches local maxima (defined as values greater than 
    // the immediately preceding and immediately following 
    // ones), and the values of x at such indexes.
    //
    // NOTE: In order to assess the first and last elements
    //       the vector x is extended with 0 values 
    //
    // ------------------------------
    //
    // Federico Miyara
    // 2007-06-28
    // 2008-08-09

    // Ensure x is a column vector
    if size(x,1)==1
        x = x(:);
        wasrow = 1;
    else 
        wasrow = 0;
    end

    // Compute a vector xdec that is positive on indexes where x
    // is decreasing. To that end, the next value is subtracted
    // from each value of x
    xdec = x - [x(2:$); 0];

    // Compute a vector xinc that is positive on indexes where x
    // is increasing. To that end, the preceding value is subtracted
    // from each value of x
    xinc = x - [0; x(1:$-1)];

    // Obtain a vector whose components are 1 if both xdec and xinc
    // are positive
    maxi = (xdec>0) .* (xinc>0);

    // Find the indexes corresponding to 1's
    i = find(maxi); 

    // If x was a row vector, so must be i
    if wasrow== 1
        i = i(:).';
    else
        i = i(:);
    end
    
    // Provide the second output argument
    y = x(i);
    // If x was a row vector, so must be y
    if wasrow==1
        y = y(:).';
    else
        y = y(:);
    end

endfunction






