function [a, e] = lpc(y, p)
    // Calculation of the LPC coeficients
    //
    // Usage:
    //       [a, e] = lpc(y, p)
    // where
    //       y: Signal vector
    //       p: Order of the linear prediction model
    //       a: Coefficients of the linear prediction model
    //       e: Total quadratic error  
    //
    // The coefficients are computed so to get minimize the 
    // prediction error by the least-square method. Once the 
    // coefficients have been computed, the autorregressive
    // transfer function
    //
    // H = 1/(1 - sum(a.*%z^(-(1:p))))
    // 
    // can be used to filter either a pulse train or white 
    // noise x to get a signal spectrally similar to the one
    // used to compute the LPC coefficients. For instance,
    //
    // y1 = filter(H.num, H.den, x);
    //
    // NOTE: The total quadratic error refers to the error
    //       between the signal and the predictor of the
    //       signal, not to be confused with the signal
    //       syhnthesized using an autoregressive filter
    //       based on the linear prediction model and a 
    //       suitabln input. This error is much bigger,
    //       especially when the signal is white noise
    //
    // References
    // Atal, B.S.; Hanauer, Susanne (1971) "Speech analysis
    //      and synthesis bylinear prediction of speech 
    //      wave"". JASA Vol 50 No 2 Part 2 pp.637-655
    // Proakis, John G.; Manolakis, Dimitris G. (1998). 
    //      "Tratamiento digital de señales. Principios, 
    //      algoritmos y aplicaciones”. Prentice Hall, 
    //      Madrid, España
    //
    // ------------------------------
    // Author: Federico Miyara
    // Date:   2020-09-27
    //         2022-08-13
    
    if 1==2
        // Test data --delete--
        // LPC order
        p = 10
        // Signal length
        N = 8192
        // Create white noise
        x = rand(N, 1, "normal");
        // Create a filter
        // Create some poles within the unit circle for stability
        pol = 0.7 + %i*linspace(-0.3, 0.3, 10)
        // Denominator polynomial in z
        AAA = real(prod(%z - pol))
        // Theoretical LPC coefficients
        aa = -coeff(AAA)($-1:-1:1)'
        // Autorregressive (AR) transfer function
        HH = 1/(1 + sum(aa'.*%z^(-(1:p))))
        // Numerator and denominator 
        BB = HH.num
        AA = HH.den
        // Filter white noise
        y = filter(BB, AA, x);
        scf();
        plot(y)
    end
    
    // Covariance matrix
    for h = 1:p
        // Row h
        for i = 1:p
            // Column i
            PHI(h,i) = sum(y(p+1-i:$-i).*y(p+1-h:$-h)); 
        end
        PHIo(h) = sum(y(p+1:$).*y(p+1-h:$-h));
    end
    
    // Atal & Hanauer's method for solving the 
    // system. A triangular matrix V such that 
    // V*V' = PHI is found using Cholesky 
    // factorization. This is more efficient 
    // and accurate than a = -(inv(PHI)*PHIo); 
    V = chol(PHI)';
    Vinv = inv(V);
    X = -Vinv*PHIo;
    a = Vinv'*X;
    
    // Total quadratic error (Atal & Hanauer)
    e = sum(y(p+1:$).^2) + a'*PHIo

    // Autoregressive transfer function
    H = 1/(1 + sum(a' .* %z^(-(1:p))));
    
    // Check whether there are unstable poles and 
    // apply a correction similar to Atal & Hanauer's
    // (Optimality will be lost). Instead of just 
    // replacing the poles by ones rejected on the 
    // unit circumference, they are reflected inside
    // the unit circle, trying to achieve similar 
    // local amplitude effect    
    pols = roots(H.den);
    unst = abs(pols)>=1;
    for i=1:length(pols)
       if unst(i)==%t
          pols(i) = pols(i)/(abs(pols(i))^2);
       end
    end
    
    if or(unst)==%t
       // Denominator with corrected coefficients
       den = real(prod(%z - pols));
       // Corrected LPC
       a = coeff(den)($-1:-1:1)';
    end
    
    
    if 1==2
        // Test data --delete--
        // LPC order
        H = 1/(1 + sum(a.*%z^(-(1:p))))
        // Numerator and denominator 
        B1 = H.num
        A1 = H.den
        // Filter white noise
        y1 = filter(B1, A1, x);
        y1 = x1(N+1:2*N);
        plot(y1)
    end
              
endfunction
