[Scilab-users] linear prediction model

[Federico Miyara]

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



-- 
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20201203/f7a2e0cd/attachment.htm>