c13-6 (Numerical Recipes in C), страница 2

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 dataset will cause some of these roots to be inside the unit circle, others outside. Insome applications, where the resulting instabilities are slowly growing and the linearprediction is not pushed too far, it is best to use the “unmassaged” LP coefficientsthat come directly out of (13.6.12).

For example, one might be extrapolating to fill ashort gap in a data set; then one might extrapolate both forwards across the gap andbackwards from the data beyond the gap. If the two extrapolations agree tolerablywell, then instability is not a problem.When instability is a problem, you have to “massage” the LP coefficients.

Youdo this by (i) solving (numerically) equation (13.6.14) for its N complex roots; (ii)moving the roots to where you think they ought to be inside or on the unit circle; (iii)reconstituting the now-modified LP coefficients. You may think that step (ii) soundsa 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 incommensurateperiods), then you will want simply to move each root zi onto the unit circle,zi → zi / |zi |(13.6.16)In other circumstances it may seem appropriate to reflect a bad root across theunit circlezi → 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 onSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).the future, imagining the unknown “future” discrepancies xi to be zero.

In thisapplication, (13.6.11) is a kind of extrapolation formula. In many situations, thisextrapolation turns out to be vastly more powerful than any kind of simple polynomialextrapolation. (By the way, you should not confuse the terms “linear prediction” and“linear extrapolation”; the general functional form used by linear prediction is muchmore complex than a straight line, or even a low-order polynomial!)However, to achieve its full usefulness, linear prediction must be constrained inone additional respect: One must take additional measures to guarantee its stability.Equation (13.6.11) is a special case of the general linear filter (13.5.1).

The conditionthat (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 polynomial568Chapter 13.Fourier and Spectral Applications#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 coefficients 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;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).identifying its time sense so that signals that should be damped as time proceeds endup growing in amplitude. The choice between (13.6.16) and (13.6.17) sometimesmight 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.

Youshould choose M to be as small as works for you, that is, you should choose it byexperimenting with your data. Try M = 5, 10, 20, 40. If you need larger M ’s thanthis, be aware that the procedure of “massaging” all those complex roots is quitesensitive to roundoff error.

Use double precision.Linear prediction is especially successful at extrapolating signals that are smoothand oscillatory, though not necessarily periodic. In such cases, linear prediction oftenextrapolates accurately through many cycles of the signal. By contrast, polynomialextrapolation 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 theheight of ocean tides, for which the fundamental 12-hour period is modulated inphase and amplitude over the course of the month and year, and for which localhydrodynamic effects may make even one cycle of the curve look rather differentin shape from a sine wave.We already remarked that equation (13.6.10) is not necessarily the best wayto estimate the covariances φk from the data set. In fact, results obtained fromlinear 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 procedurefor increasing the order M by one unit at a time, at each stage re-estimating thecoefficients 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 methodis implemented in the following routine [1,2] for estimating the LP coefficients djof a data set.13.6 Linear Prediction and Linear Predictive Coding569}nrerror("never get here in memcof.");}Here are procedures for rendering the LP coefficients stable (if you choose todo 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 complexroots of a polynomial.#include <math.h>#include "complex.h"#define NMAX 100#define ZERO Complex(0.0,0.0)#define ONE Complex(1.0,0.0)Largest expected value of m.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 amodified 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 roota[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 rootsa[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 bed[m-j] = -a[j].r;real, so we need only return the real part asnew LP coefficients.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (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 muntil the desired value is reached. At this point in the algorithm, one could returnthe 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];}570Chapter 13.Fourier and Spectral Applications#include "nrutil.h"reg=vector(1,m);for (j=1;j<=m;j++) reg[j]=data[ndata+1-j];for (j=1;j<=nfut;j++) {discrp=0.0;This is where you would put in a known discrepancy if you were reconstructing afunction by linear predictive coding rather than extrapolating a function by linear prediction.