568 Chapter 13. Fourier and Spectral Applications identifying its time sense so that signals that should be damped as time proceeds end up growing in amplitude.The choice between(13.6.16)and(13.6.17)sometimes might as well be based on voodoo.We prefer (13.6.17). Also magical is the choice of M,the number of LP coefficients to use.You should choose M to be as small as works for you,that is,you should choose it by experimenting with your data.Try M=5.10,20,40.If you need larger M's than this,be aware that the procedure of"massaging"all those complex roots is quite sensitive to roundoff error.Use double precision. Linear prediction is especially successful at extrapolating signals that are smooth and oscillatory,though not necessarily periodic.In such cases,linear prediction often extrapolates accurately through many cycles of the signal.By contrast,polynomial extrapolation in general becomes seriously inaccurate after at most a cycle or two. 鱼 19881992 A prototypical example of a signal that can successfully be linearly predicted is the height of ocean tides,for which the fundamental 12-hour period is modulated in -200 phase and amplitude over the course of the month and year,and for which local hydrodynamic effects may make even one cycle of the curve look rather different from NUMERICAL RECIPESI in shape from a sine wave. We already remarked that equation (13.6.10)is not necessarily the best way (Nort server to estimate the covariances o from the data set.In fact,results obtained from linear prediction are remarkably sensitive to exactly how the o's are estimated. America computer, make one paper UnN电.t THE One particularly good method is due to Burg [1],and involves a recursive procedure ART for increasing the order M by one unit at a time,at each stage re-estimating the coefficients d,j=1,...,M so as to minimize the residual in equation(13.6.13). Programs Although further discussion of the Burg method is beyond our scope here,the method is implemented in the following routine [1.2]for estimating the LP coefficients dj of a data set 会 #include <math.h> #include "nrutil.h' void memcof(float data[],int n,int m,float *xms,float d[) Given a real vector of data[1..n],and given m,this routine returns m linear prediction coef- OF SCIENTIFIC COMPUTING(ISBN 0-521- ficients as d[1..m],and returns the mean square discrepancy as xms. int k,j,i; float p=0.0,*wk1,*wk2,*wkm; v@cambridge.org er eba 1988-1992 by Numerical Recipes -431085 yki-vector(1,n); wk2-vector(1,n); wkm=vector(1,m); for (j=1;j<=n;j++)p +SQR(data[j]); Software. 来xmg=o/n: wk1[1]=data[1]; wk2[n-1]=data[n]; (outside North America) for(j=2;j<=n-1;j++)( wk1[i]=data[i]; wk2[j-1]=data[j]; for (k=1;k<=m;k++){ float num=0.0,denom=0.0; for(j=1;j<=(n-k);j+)1 num +wk1[j]*wk2[i] denom +SQR(wk1[j])+SQR(wk2[]); d[k]=2.0*num/denom;568 Chapter 13. Fourier and Spectral Applications 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). identifying its time sense so that signals that should be damped as time proceeds end up growing in amplitude. The choice between (13.6.16) and (13.6.17) sometimes might as well be based on voodoo. We prefer (13.6.17). Also magical is the choice of M, the number of LP coefficients to use. You should choose M to be as small as works for you, that is, you should choose it by experimenting with your data. Try M = 5, 10, 20, 40. If you need larger M’s than this, be aware that the procedure of “massaging” all those complex roots is quite sensitive to roundoff error. Use double precision. Linear prediction is especially successful at extrapolating signals that are smooth and oscillatory, though not necessarily periodic. In such cases, linear prediction often extrapolates accurately through many cycles of the signal. By contrast, polynomial extrapolation in general becomes seriously inaccurate after at most a cycle or two. A prototypical example of a signal that can successfully be linearly predicted is the height of ocean tides, for which the fundamental 12-hour period is modulated in phase and amplitude over the course of the month and year, and for which local hydrodynamic effects may make even one cycle of the curve look rather different in shape from a sine wave. We already remarked that equation (13.6.10) is not necessarily the best way to estimate the covariances φk from the data set. In fact, results obtained from linear prediction are remarkably sensitive to exactly how the φ k’s are estimated. One particularly good method is due to Burg [1], and involves a recursive procedure for increasing the order M by one unit at a time, at each stage re-estimating the coefficients dj , j = 1,...,M so as to minimize the residual in equation (13.6.13). Although further discussion of the Burg method is beyond our scope here, the method is implemented in the following routine [1,2] for estimating the LP coefficients dj of a data set. #include <math.h> #include "nrutil.h" void memcof(float data[], int n, int m, float *xms, float d[]) Given a real vector of data[1..n], and given m, this routine returns m linear prediction coef- ficients as d[1..m], and returns the mean square discrepancy as xms. { int k,j,i; float p=0.0,*wk1,*wk2,*wkm; wk1=vector(1,n); wk2=vector(1,n); wkm=vector(1,m); for (j=1;j<=n;j++) p += SQR(data[j]); *xms=p/n; wk1[1]=data[1]; wk2[n-1]=data[n]; for (j=2;j<=n-1;j++) { wk1[j]=data[j]; wk2[j-1]=data[j]; } for (k=1;k<=m;k++) { float num=0.0,denom=0.0; for (j=1;j<=(n-k);j++) { num += wk1[j]*wk2[j]; denom += SQR(wk1[j])+SQR(wk2[j]); } d[k]=2.0*num/denom;