f13-6

13.6LinearPredictionandLinearPredictiveCoding 557Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).13.6 Linear Prediction and Linear Predictive Coding Webeginwithaverygeneralformulationthatwillallowustomakeconnections to various special cases. Let {y/prime α}be a set of measured values for some underlying set of true values of a quantity y, denoted {yα}, related to these true values by the addition of random noise, y/prime α=yα+nα (13.6.1 ) (compareequation 13.3.2, with a somewhat different notation). Our use of a Greek subscript to index the members of the set is meant to indicate that the data points are not necessarily equally spaced along a line, or even ordered: they might be“random”pointsinthree-dimensionalspace,forexample. Now,supposewewantto constructthe “best”estimate ofthe truevalueof someparticularpoint y ⋆as alinear combination of the known, noisy, values. Writing y⋆=/summationdisplay αd⋆αy/prime α+x⋆ (13.6.2 ) we want to find coefficients d⋆αthat minimize, in some way, the discrepancy x⋆. Thecoefficients d⋆αhavea“star”subscripttoindicatethattheydependonthechoice of point y⋆. Later, we might want to let y⋆be one of the existing yα’s. In that case, our problem becomes one of optimal filtering or estimation, closely related to thediscussion in §13.3. On the other hand, we might want y ⋆to be a completely new point. In that case, our problem will be one of linear prediction . A natural way to minimize the discrepancy x⋆is in the statistical mean square sense. Ifanglebracketsdenotestatisticalaverages,thenweseek d⋆α’sthatminimize /angbracketleftbig x2 ⋆/angbracketrightbig =/angbracketleftBigg/bracketleftbigg/summationdisplay αd⋆α(yα+nα)−y⋆/bracketrightbigg2/angbracketrightBigg =/summationdisplay αβ(/angbracketleftyαyβ/angbracketright+/angbracketleftnαnβ/angbracketright)d⋆αd⋆β−2/summationdisplay α/angbracketlefty⋆yα/angbracketrightd⋆α+/angbracketleftbig y2 ⋆/angbracketrightbig(13.6.3 ) Here we have used the fact that noise is uncorrelatedwith signal, e.g., /angbracketleftnαyβ/angbracketright=0. The quantities /angbracketleftyαyβ/angbracketrightand/angbracketlefty⋆yα/angbracketrightdescribe the autocorrelation structure of the underlying data. We have already seen an analogous expression, (13.2.2), for thecase of equally spaced data points on a line; we will meet correlation several times againinitsstatisticalsenseinChapters14and15. Thequantities /angbracketleftn αnβ/angbracketrightdescribethe autocorrelationpropertiesof the noise. Often, forpoint-to-pointuncorrelatednoise, we have /angbracketleftnαnβ/angbracketright=/angbracketleftbig n2 α/angbracketrightbig δαβ. It is convenient to think of the various correlation quantities as comprising matrices and vectors, φαβ≡/angbracketleftyαyβ/angbracketright φ⋆α≡/angbracketlefty⋆yα/angbracketright ηαβ≡/angbracketleftnαnβ/angbracketrightor/angbracketleftbig n2 α/angbracketrightbig δαβ (13.6.4 ) Setting the derivative of equation (13.6.3) with respect to the d⋆α’s equal to zero, one readily obtains the set of linear equations, /summationdisplay β[φαβ+ηαβ]d⋆β=φ⋆α (13.6.5 ) 558 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If we write the solution as a matrix inverse, then the estimation equation (13.6.2) becomes, omitting the minimized discrepancy x⋆, y⋆≈/summationdisplay αβφ⋆α[φµν+ηµν]−1 αβy/prime β (13.6.6 ) Fromequations(13.6.3)and(13.6.5)onecanalsocalculatetheexpectedmeansquare value of the discrepancy at its minimum, denoted/angbracketleftbig x2 ⋆/angbracketrightbig 0, /angbracketleftbig x2 ⋆/angbracketrightbig 0=/angbracketleftbig y2 ⋆/angbracketrightbig −/summationdisplay βd⋆βφ⋆β=/angbracketleftbig y2 ⋆/angbracketrightbig −/summationdisplay αβφ⋆α[φµν+ηµν]−1 αβφ⋆β (13.6.7 ) A final general result tells how much the mean square discrepancy/angbracketleftbig x2 ⋆/angbracketrightbig is increasedifweusetheestimationequation(13.6.2)notwiththebestvalues d⋆β,b ut with some other values /hatwided⋆β. The above equations then imply /angbracketleftbig x2 ⋆/angbracketrightbig =/angbracketleftbig x2 ⋆/angbracketrightbig 0+/summationdisplay αβ(/hatwided⋆α−d⋆α)[φαβ+ηαβ](/hatwided⋆β−d⋆β)( 13.6.8 ) Since the second term is a pure quadratic form, we see that the increase in the discrepancyis only second order in any error made in estimating the d⋆β’s. Connection toOptimalFiltering If we change “star” to a Greek index, say γ, then the above formulas describe optimal filtering, generalizing the discussion of §13.3. One sees, for example, that if the noise amplitudes nαgo to zero, so likewise do the noise autocorrelations ηαβ, and, canceling a matrix times its inverse, equation (13.6.6) simply becomes yγ=y/prime γ. Another special case occurs if the matrices φαβandηαβare diagonal. In that case, equation (13.6.6) becomes yγ=φγγ φγγ+ηγγy/prime γ (13.6.9 ) whichisreadilyrecognizableasequation(13.3.6)with S2→φγγ,N2→ηγγ. What is going on is this: For the case of equally spaced data points, and in the Fourierdomain, autocorrelations become simply squares of Fourier amplitudes (Wiener- Khinchin theorem, equation 12.0.12), and the optimal filter can be constructed algebraically, as equation (13.6.9), without inverting any matrix. Moregenerally,inthe timedomain,oranyotherdomain,anoptimalfilter (one that minimizes the square of the discrepancy from the underlying true value in thepresenceofmeasurementnoise)canbeconstructedbyestimatingtheautocorrelation matrices φ αβandηαβ, and applying equation (13.6.6) with ⋆→γ. (Equation 13.6.8 is in fact the basis for the §13.3’s statement that even crude optimal filtering can be quite effective.) 13.6LinearPredictionandLinearPredictiveCoding 559Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Linear Prediction Classical linear prediction specializes to the case where the data points yβ are equally spaced along a line, yi,i=1,2,...,N, and we want to use M consecutivevaluesof yito predictan M+1st. Stationarityis assumed. Thatis, the autocorrelation /angbracketleftyjyk/angbracketrightis assumed to dependonly on the difference |j−k|, and not onjorkindividually, so that the autocorrelation φhas only a single index, φj≡/angbracketleftyiyi+j/angbracketright≈1 N−jN−j/summationdisplay i=1yiyi+j (13.6.10 ) Here, the approximate equality shows one way to use the actual data set values to estimatetheautocorrelationcomponents. (Infact,thereisabetterwaytomakethese estimates; see below.) In thesituationdescribed,theestimationequation(13.6.2)is yn=M/summationdisplay j=1djyn−j+xn (13.6.11 ) (compareequation13.5.1)andequation(13.6.5)becomesthesetof Mequationsfor theMunknown dj’s, now called the linear prediction (LP) coefficients , M/summationdisplay j=1φ|j−k|dj=φk (k=1,...,M )( 13.6.12 ) Notice that while noise is not explicitly included in the equations, it is properly accounted for, ifit is point-to-point uncorrelated: φ0, as estimated by equation (13.6.10)using measured values y/prime i,actuallyestimatesthediagonalpartof φαα+ηαα, above. The mean squarediscrepancy/angbracketleftbig x2 n/angbracketrightbig is estimated by equation(13.6.7)as /angbracketleftbig x2 n/angbracketrightbig =φ0−φ1d1−φ2d2−···− φMdM (13.6.13 ) To use linear prediction, we first compute the dj’s, using equations (13.6.10) and (13.6.12). We then calculate equation (13.6.13) or, more concretely, apply (13.6.11) to the known record to get an idea of how large are the discrepancies xi. If the discrepanciesare small, then we can continueapplying(13.6.11)rightoninto the future, imagining the unknown “future” discrepancies xito be zero. In this application, (13.6.11) is a kind of extrapolation formula. In many situations, thisextrapolationturnsouttobevastlymorepowerfulthananykindofsimplepolynomial extrapolation. (Bytheway,youshouldnotconfusetheterms“linearprediction”and “linearextrapolation”;thegeneralfunctionalformusedbylinearpredictionis much more complex than a straight line, or even a low-order polynomial!) However,toachieveitsfullusefulness,linearpredictionmustbeconstrainedin oneadditionalrespect: One must take additionalmeasures to guaranteeits stability. Equation(13.6.11)isaspecialcaseofthegenerallinearfilter(13.5.1). Thecondition that (13.6.11) be stable as a linear predictor is precisely that given in equations(13.5.5) and (13.5.6), namely that the characteristic polynomial z N−N/summationdisplay j=1djzN−j=0 ( 13.6.14 ) 560 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).have all Nof its roots inside the unit circle, |z|≤ 1( 13.6.15 ) There is no guaranteethat the coefficients producedbyequation(13.6.12)will have this property. If the data contain many oscillations without any particular trendtowards increasing or decreasing amplitude, then the complex roots of (13.6.14) will generally all be rather close to the unit circle. The finite length of the data set will cause some of these roots to be inside the unit circle, others outside. Insomeapplications,wheretheresultinginstabilitiesareslowlygrowingandthelinear prediction is not pushed too far, it is best to use the “unmassaged” LP coefficients thatcomedirectlyoutof(13.6.12). Forexample,onemightbeextrapolatingtofill ashort gap in a data set; then one might extrapolateboth forwardsacross the gap and backwards from the data beyond the gap. If the two extrapolations agree tolerably well, then instability is not a problem. When instability isa problem,youhave to “massage” the LP coefficients. You do this by (i) solving (numerically) equation (13.6.14)for its Ncomplex roots; (ii) movingtherootstowhereyouthinktheyoughttobeinsideorontheunitcircle;(iii) reconstitutingthenow-modifiedLPcoefficients. Youmaythinkthatstep(ii)sounds a little vague. It is. There is no “best” procedure. If you think that your signalis truly a sum of undamped sine and cosine waves (perhaps with incommensurate periods), then you will want simply to moveeach root z ionto the unit circle, zi→zi/|zi| (13.6.16 ) In other circumstances it may seem appropriate to reflect a bad root across the unit circle zi→ 1/zi*( 13.6.17 ) This alternative has the property that it preserves the amplitude of the output of (13.6.11) when it is driven by a sinusoidal set of xi’s. It assumes that (13.6.12) has correctly identified the spectral width of a resonance, but only slipped up onidentifyingitstimesensesothatsignalsthatshouldbedampedastimeproceedsend 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 Mto be as small as works for you, that is, you should choose it by experimentingwith 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. Linearpredictionisespeciallysuccessfulatextrapolatingsignalsthataresmooth andoscillatory,thoughnotnecessarilyperiodic. Insuchcases,linearpredictionoften 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 prototypicalexample 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 13.6LinearPredictionandLinearPredictiveCoding 561Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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 φkfrom the data set. In fact, results obtained from linear prediction are remarkably sensitive to exactly how the φk’s are estimated. One particularlygood methodis due to Burg [1], and involvesa recursive procedure for increasing the order Mby one unit at a time, at each stage re-estimating the coefficients dj,j=1,...,Mso as to minimize the residual in equation (13.6.13). AlthoughfurtherdiscussionoftheBurgmethodisbeyondourscopehere,themethod is implemented in the following routine [1,2]for estimating the LP coefficients dj of a data set. SUBROUTINE memcof(data,n,m,xms,d) INTEGER m,n,MMAX,NMAXREAL xms,d(m),data(n)PARAMETER (MMAX=60,NMAX=2000) Given a real vector of data(1:n) ,a n dg i v e n m, this routine returns mlinear prediction coefficients as d(1:m) , and returns the mean square discrepancy as xms. INTEGER i,j,k REAL denom,p,pneum,wk1(NMAX),wk2(NMAX),wkm(MMAX) if (m.gt.MMAX.or.n.gt.NMAX) pause ’workspace too small in memcof’p=0.do 11j=1,n p=p+data(j)**2 enddo 11 xms=p/nwk1(1)=data(1) wk2(n-1)=data(n) do 12j=2,n-1 wk1(j)=data(j) wk2(j-1)=data(j) enddo 12 do17k=1,m pneum=0. denom=0. do13j=1,n-k pneum=pneum+wk1(j)*wk2(j) denom=denom+wk1(j)**2+wk2(j)**2 enddo 13 d(k)=2.*pneum/denomxms=xms*(1.-d(k)**2) do 14i=1,k-1 d(i)=wkm(i)-d(k)*wkm(k-i) enddo 14 The algorithm is recursive, building up the answer for larger and larger values of muntil the desired value is reached. At this point in the algorithm, one could return the vectordand scalar xmsfor a set of LP coefficients with k(rather than m)terms. if(k.eq.m)return do 15i=1,k wkm(i)=d(i) enddo 15 do16j=1,n-k-1 wk1(j)=wk1(j)-wkm(k)*wk2(j) wk2(j)=wk2(j+1)-wkm(k)*wk1(j+1) enddo 16 enddo 17 pause ’never get here in memcof’END 562 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).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 ormassaged LP coefficients. The routine zroots(§9.5) is used to find all complex roots of a polynomial. SUBROUTINE fixrts(d,m) INTEGER m,MMAX REAL d(m)PARAMETER (MMAX=100) Largest expected value of m. C USES zroots 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 returnsa modified set of coefficients d(1:m) . INTEGER i,j LOGICAL polishCOMPLEX a(MMAX),roots(MMAX)a(m+1)=cmplx(1.,0.) do 11j=m,1,-1 Set up complex coefficients for polynomial root finder. a(j)=cmplx(-d(m+1-j),0.) enddo 11 polish=.true.call zroots(a,m,roots,polish) Find all the roots. do 12j=1,m Look for a... if(abs(roots(j)).gt.1.)then root outside the unit circle, roots(j)=1./conjg(roots(j)) and reflect it back inside. endif enddo 12 a(1)=-roots(1) Now reconstruct the polynomial coefficients, a(2)=cmplx(1.,0.)do 14j=2,m by looping over the roots a(j+1)=cmplx(1.,0.) do13i=j,2,-1 and synthetically multiplying. a(i)=a(i-1)-roots(j)*a(i) enddo 13 a(1)=-roots(j)*a(1) enddo 14 do15j=1,m The polynomial coefficients are guaranteed to be real, d(m+1-j)=-real(a(j)) so we need only return the real part as new LP coefficients. enddo 15 return END SUBROUTINE predic(data,ndata,d,m,future,nfut) INTEGER ndata,nfut,m,MMAXREAL d(m),data(ndata),future(nfut) PARAMETER (MMAX=100) Given data(1:ndata) , and given the data’s LP coefficients d(1:m) , this routine applies equation (13.6.11)to predict the next nfutdata points, which it returns in the array future(1:nfut) . Note that the routine references only the last mvalues of data,a s initial values for the prediction. Parameter: MMAXis the largest expected value of m. INTEGER j,k REAL discrp,sum,reg(MMAX) do11j=1,m reg(j)=data(ndata+1-j) enddo 11 do14j=1,nfut discrp=0. This is where you would put in a known discrepancy if you were reconstructing a function by linear predictive coding rather than extrapolating a function by linear prediction. See text.sum=discrp do12k=1,m 13.6LinearPredictionandLinearPredictiveCoding 563Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).sum=sum+d(k)*reg(k) enddo 12 do13k=m,2,-1 [If you want to implement circular arrays, you can avoid this shifting of coefficients!] reg(k)=reg(k-1) enddo 13 reg(1)=sumfuture(j)=sum enddo 14 return END Removingthe Bias inLinearPrediction You might expect that the sum of the dj’s in equation (13.6.11) (or, more generally, in equation 13.6.2)should be 1, so that (e.g.) adding a constant to all the data points yiyields a prediction that is increased by the same constant. However, thedj’s do not sum to 1 but, in general, to a value slightly less than one. This fact revealsasubtlepoint,thattheestimatorofclassicallinearpredictionisnot unbiased, even thoughit does minimize the mean square discrepancy. At any place where themeasured autocorrelation does not imply a better estimate, the equations of linear prediction tend to predict a value that tends towards zero. Sometimes, that is just what you want. If the process that generates the y i’s in fact has zero mean, then zero is the best guess absent other information. At other times, however, this behavior is unwarranted. If you have data that show only small variations around a positive value, you don’t want linear predictions that droop towards zero. Often it is a workable approximation to subtract the mean off your data set, performthe linear prediction,and then add the mean back. This procedurecontains the germ of the correct solution; but the simple arithmetic mean is not quite the correctconstanttosubtract. Infact,anunbiasedestimatorisobtainedbysubtractingfrom every data point an autocorrelation-weightedmean defined by [3,4] y≡/summationdisplay β[φµν+ηµν]−1 αβyβ/slashbigg/summationdisplay αβ[φµν+ηµν]−1 αβ(13.6.18 ) With this subtraction,the sum of the LP coefficients shouldbe unity, up to roundoff and differences in how the φk’s are estimated. Linear Predictive Coding(LPC) A different, though related, method to which the formalism above can be applied is the “compression” of a sampled signal so that it can be stored more compactly. The original form should be exactlyrecoverable from the compressed version. Obviously, compression can be accomplished only if there is redundancy in the signal. Equation (13.6.11) describes one kind of redundancy: It says that the signal, except for a small discrepancy, is predictable from its previous valuesand from a small number of LP coefficients. Compression of a signal by the use of (13.6.11) is thus called linear predictive coding ,o rLPC. ThebasicideaofLPC(initssimplestform)istorecordasacompressedfile(i) the number of LP coefficients M, (ii) their Mvalues, e.g., as obtained by memcof, 564 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).(iii) the first Mdata points, and then (iv) for each subsequent data point only its residual discrepancy xi(equation 13.6.1). When you are creating the compressed file, youfindthe residualbyapplying(13.6.1)to the previous Mpoints,subtracting the sum fromthe actual valueof the currentpoint. Whenyou are reconstructingthe originalfile,youaddtheresidualbackin,atthepointindicatedintheroutine predic. Itmaynotbeobviouswhythereisanycompressionatallinthisscheme. After all,wearestoringonevalueofresidualperdatapoint! Whynotjuststoretheoriginal data point? The answer depends on the relative sizes of the numbers involved. The residual is obtained by subtracting two very nearly equal numbers (the data and the linearprediction). Therefore,thediscrepancytypicallyhasonlyaverysmallnumberof nonzero bits. These can be stored in a compressed file. How do you do it in a high-level language? Here is one way: Scale your data to have integer values, say between +1000000 and−1000000(supposingthatyouneedsixsignificantfigures). Modifyequation(13.6.1)byenclosingthesumterminan“integerpartof”operator. The discrepancy will now, by definition, be an integer. Experiment with different values of M, to find LP coefficientsthat make the rangeof the discrepancyas small asyoucan. Ifyoucangettowithinarangeof ±127(andinourexperiencethisisnot at all difficult) then you can write it to a file as a single byte. This is a compression factor of 4, compared to 4-byte integer or floating formats. Notice that the LP coefficientsare computedusing the quantized data, and that the discrepancy is also quantized, i.e., quantization is done both outside and insidethe LPC loop. If you are careful in following this prescription, then, apart from the initial quantization of the data, you will not introduce even a single bit of roundoff error into the compression-reconstructionprocess: While the evaluation of the sumin(13.6.11)mayhaveroundofferrors,theresidualthat youstoreis thevaluewhich, whenaddedbacktothesum,gives exactlytheoriginal(quantized)datavalue. Notice also that you do not need to massage the LP coefficients for stability; by addingthe residualbackintoeachpoint,youneverdepartfromtheoriginaldata,soinstabilities cannot grow. There is therefore no need for fixrts, above. Look at §20.4 to learn about Huffmancoding , which will further compress the residualsbytakingadvantageofthefactthatsmallervaluesofdiscrepancywilloccur more often than larger values. A very primitive version of Huffman coding wouldbe this: If most of the discrepancies are in the range ±127,but an occasionalone is outside,thenreservethevalue127tomean“outofrange,”andthenrecordonthefile (immediately following the 127) a full-word value of the out-of-range discrepancy. §20.4 explains how to do much better. There are many variant proceduresthat all fall under the rubric of LPC. •If the spectral character of the data is time-variable, then it is best not to use a single set of LP coefficients for the whole data set, but rather to partition the data into segments, computing and storing different LPcoefficients for each segment. •If the data are really well characterized by their LP coefficients, and you cantoleratesomesmallamountoferror,thendon’tbotherstoringallofthe residuals. Justdolinearpredictionuntilyouareoutsideoftolerances,then reinitialize(using Msequentialstoredresiduals)andcontinuepredicting. •Insomeapplications,mostnotablyspeechsynthesis,onecaresonlyabout the spectral content of the reconstructed signal, not the relative phases. In this case, one need not store any starting values at all, but only the 13.7MaximumEntropy(AllPoles)Method 565Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).LP coefficients for each segment of the data. The output is reconstructed by driving these coefficients with initial conditions consisting of all zerosexcept for one nonzero spike. A speech synthesizer chip may have of order10LPcoefficients,whichchangeperhaps20to50timespersecond. •Somepeoplebelievethatit isinterestingtoanalyzeasignalbyLPC, even when the residuals x iarenotsmall. The xi’s are then interpreted as the underlying“input signal” which, when filtered throughthe all-poles filter defined by the LP coefficients (see §13.7), produces the observed“output signal.” LPC reveals simultaneously,it is said, the nature of the filter and theparticularinputthatisdrivingit. Weareskepticaloftheseapplications;the literature, however, is full of extravagant claims. CITED REFERENCES AND FURTHER READING: Childers, D.G. (ed.) 1978, Modern Spectrum Analysis (New York: IEEE Press), especially the paper by J. Makhoul (reprinted from Proceedings of the IEEE , vol. 63, p. 561, 1975). Burg, J.P. 1968, reprinted in Childers, 1978. [1] Anderson, N. 1974, reprinted in Childers, 1978. [2] Cressie,N.1991,in SpatialStatisticsandDigitalImageAnalysis (Washington:NationalAcademy Press). [3] Press, W.H., and Rybicki, G.B. 1992, Astrophysical Journal , vol. 398, pp. 169–176. [4] 13.7 Power Spectrum Estimation by the Maximum Entropy (All Poles) Method The FFT is not the only way to estimate the power spectrum of a process, nor is it necessarily the best way for all purposes. To see how one might devise another method,let us enlarge our view for a moment, so that it includes not only real frequencies in theNyquist interval −f c<f <f c, but also the entire complex frequency plane. From that vantage point, let us transform the complex f-plane to a new plane, called the z-transform planeorz-plane, by the relation z≡e2πif∆(13.7.1 ) where ∆is,asusual,thesamplingintervalinthetimedomain. NoticethattheNyquistinterval on the real axis of the f-plane maps one-to-one onto the unit circle in the complex z-plane. If we now compare (13.7.1) to equations (13.4.4) and (13.4.6), we see that the FFT power spectrum estimate (13.4.5) for any real sampled function ck≡c(tk)can be written, except for normalization convention, as P(f)=/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleN/2−1/summationdisplay k=−N/2ckzk/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle2 (13.7.2 ) Ofcourse, (13.7.2)isnotthe truepower spectrumoftheunderlying function c(t),butonlyan estimate. Wecanseeintworelatedwayswhytheestimateisnotlikelytobeexact. First,inthetimedomain,theestimateisbasedononlyafiniterangeofthefunction c(t)whichmay,forall weknow,havecontinuedfrom t=−∞to∞. Second,inthe z-planeofequation(13.7.2),the finiteLaurentseriesoffers,ingeneral,only anapproximation toageneralanalytic functionofz. In fact,a formal expression for representing “true” power spectra (up to normalization) is P(f)=/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle ∞/summationdisplay k=−∞ckzk/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle2 (13.7.3 )