[Scilab-users] linear prediction model

classic Classic list List threaded Threaded
1 message Options
fmiyara fmiyara
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] linear prediction model


Dear All,

Is there a script or function somewhere which calculates the Linear prediction model (LPC) of a signal frame?

I tried to implement a function (below) without using the algorithm of Levinson & Durbin but it doesn't seem to work properly. seemingly the problem is due to inaccuracy of the matrix equation solution.

Regards,

Federico Miyara





function a=lpc(x, p)
    // Calculation of the LPC coeficients
    // Usage:
    //       a = lpc(x, p)
    // where
    //       x: Signal vector
    //       p: Order of the linear prediction model
    //       a: Coefficients of the linear prediction model
    //
    // The coefficients are computed so to get minimize the 
    // prediction error by the least-square method
    //
    // ------------------------------
    // Author: Federico Miyara
    // Date:   20202-09-27
    
    if 1==2
        // Test data --delete--
        // LPC order
        p = 10
        // Signal length
        N = 512
        // Create whiter noise
        exec whitenoise.sci;
        y = whitenoise(2*N-1,1);
        // 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
        x = filter(BB, AA, y);
        x = x(N+1:2*N);
        plot(x)
    end
    
    
    for h = 1:p
        // Row h
        for i = 1:p
            // Column i
            A(h,i) = sum(x(p+1-i:$-i).*x(p+1-h:$-h)); 
        end
        B(h) = sum(x(p+1:$).*x(p+1-h:$-h));
    end
    
    a = ((inv(A)*B))';
    
    
    if 1==2
        H = 1/(1 - sum(a.*%z^(-(1:p))))
        // Numerator and denominator 
        B1 = H.num
        A1 = H.den
        // Filter white noise
        x1 = filter(B1, A1, y);
        x1 = x1(N+1:2*N);
        plot(x1)
    end
              
endfunction

Libre de virus. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users