13.6 Linear Prediction and Linear Predictive Coding 569 *xms*=(1.0-S0R(d[k]); for(1=1:1<=(k-1):1++) d[i]=wkm[i]-d[k]*wkm[k-i]; The algorithm is recursive,building up the answer for larger and larger values of m until the desired value is reached.At this point in the algorithm,one could return the vector d and scalar xms for a set of LP coefficients with k (rather than m) terms. 1f(k=m)[ free_vector(wkm,1,m); free_vector(wk2,1,n); free_vector(wk1,1,n); return: for (i=1;i<=k;i++)wkm[i]=d[i]; for(j=1;j<=(n-k-1):j+){ wk1[i]-wkm[k]*wk2[i]: granted for 19881992 k2[j]=wk2[j+1]-wkm[k]*wk1[j+1]; 1600 (including this one) nrerror("never get here in memcof."); Cambridge NUMERICAL RECIPES Here are procedures for rendering the LP coefficients stable(if you choose to (Nort server to make 令 do so).and for extrapolating a data set by linear prediction.using the original or massaged LP coefficients.The routine zroots (89.5)is used to find all complex merica computer THE roots of a polynomial ART Programs #include <math.h> #include "complex.h" #define NMAX 100 Largest expected value of m. OF SCIENTIFIC #define ZERO Complex(0.0,0.0) #define ONE Complex(1.0,0.0) to dir void fixrts(float d],int m) Given the LP coefficients d[1..m],this routine finds all roots of the characteristic polynomial (13.6.14),reflects any roots that are outside the unit circle back inside,and then returns a modified set of coefficients d[1..m] 1920 COMPUTING(ISBN void zroots(fcomplex a],int m,fcomplex roots[],int polish); int i,j,polish; Numerical 1021 fcomplex a[NMAX],roots [NMAX]; a[m]=ONE; 43108 for(j=m-1;j>=0;j--） Set up complex coefficients for polynomial root Recipes a[j]=Complex(-d[m-j],0.0); finder. polish=1; (outside zroots(a,m,roots,polish); Find all the roots Software. for (j=1;j<=mij++) Look for a.. if (Cabs(roots[j])>1.0) root outside the unit circle, roots [j]=Cdiv(ONE,Conjg(roots [j]));and reflect it back inside a [O]=Csub(ZERO,roots[1]); Now reconstruct the polynomial coefficients, a[1]=0WE; for (j=2;j<=m;j++) by looping over the roots a[j]=ONE; for(1=j;1>=2:i=-) and synthetically multiplying. a[i-1]=Csub(a[i-2],Cmul(roots [j],a[i-1])); a[O]=Csub(ZERO,Cmul(roots[j],a [O])); for(j=0;j<=m-1;j+) The polynomial coefficients are guaranteed to be d[m-j]=-a[j].r; real,so we need only return the real part as new LP coefficients13.6 Linear Prediction and Linear Predictive Coding 569 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). *xms *= (1.0-SQR(d[k])); for (i=1;i<=(k-1);i++) d[i]=wkm[i]-d[k]*wkm[k-i]; The algorithm is recursive, building up the answer for larger and larger values of m until the desired value is reached. At this point in the algorithm, one could return the vector d and scalar xms for a set of LP coefficients with k (rather than m) terms. if (k == m) { free_vector(wkm,1,m); free_vector(wk2,1,n); free_vector(wk1,1,n); return; } for (i=1;i<=k;i++) wkm[i]=d[i]; for (j=1;j<=(n-k-1);j++) { wk1[j] -= wkm[k]*wk2[j]; wk2[j]=wk2[j+1]-wkm[k]*wk1[j+1]; } } nrerror("never get here in memcof."); } Here are procedures for rendering the LP coefficients stable (if you choose to do so), and for extrapolating a data set by linear prediction, using the original or massaged LP coefficients. The routine zroots (§9.5) is used to find all complex roots of a polynomial. #include <math.h> #include "complex.h" #define NMAX 100 Largest expected value of m. #define ZERO Complex(0.0,0.0) #define ONE Complex(1.0,0.0) void fixrts(float d[], int m) Given the LP coefficients d[1..m], this routine finds all roots of the characteristic polynomial (13.6.14), reflects any roots that are outside the unit circle back inside, and then returns a modified set of coefficients d[1..m]. { void zroots(fcomplex a[], int m, fcomplex roots[], int polish); int i,j,polish; fcomplex a[NMAX],roots[NMAX]; a[m]=ONE; for (j=m-1;j>=0;j--) Set up complex coefficients for polynomial root a[j]=Complex(-d[m-j],0.0); finder. polish=1; zroots(a,m,roots,polish); Find all the roots. for (j=1;j<=m;j++) Look for a... if (Cabs(roots[j]) > 1.0) root outside the unit circle, roots[j]=Cdiv(ONE,Conjg(roots[j])); and reflect it back inside. a[0]=Csub(ZERO,roots[1]); Now reconstruct the polynomial coefficients, a[1]=ONE; for (j=2;j<=m;j++) { by looping over the roots a[j]=ONE; for (i=j;i>=2;i--) and synthetically multiplying. a[i-1]=Csub(a[i-2],Cmul(roots[j],a[i-1])); a[0]=Csub(ZERO,Cmul(roots[j],a[0])); } for (j=0;j<=m-1;j++) The polynomial coefficients are guaranteed to be real, so we need only return the real part as new LP coefficients. d[m-j] = -a[j].r; }