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

See text.sum=discrp;for (k=1;k<=m;k++) sum += d[k]*reg[k];for (k=m;k>=2;k--) reg[k]=reg[k-1];[If you want to implement circularfuture[j]=reg[1]=sum;arrays, you can avoid this shift}ing of coefficients.]free_vector(reg,1,m);}Removing the Bias in Linear PredictionYou might expect that the sum of the dj ’s in equation (13.6.11) (or, moregenerally, in equation 13.6.2) should be 1, so that (e.g.) adding a constant to all thedata points yi yields a prediction that is increased by the same constant. However,the dj ’s do not sum to 1 but, in general, to a value slightly less than one. This factreveals a subtle point, that the estimator of classical linear prediction is not unbiased,even though it does minimize the mean square discrepancy.

At any place where themeasured autocorrelation does not imply a better estimate, the equations of linearprediction tend to predict a value that tends towards zero.Sometimes, that is just what you want. If the process that generates the yi ’sin fact has zero mean, then zero is the best guess absent other information. Atother times, however, this behavior is unwarranted. If you have data that showonly small variations around a positive value, you don’t want linear predictionsthat droop towards zero.Often it is a workable approximation to subtract the mean off your data set,perform the linear prediction, and then add the mean back.

This procedure containsthe germ of the correct solution; but the simple arithmetic mean is not quite thecorrect constant to subtract. In fact, an unbiased estimator is obtained by subtractingfrom every data point an autocorrelation-weighted mean defined by [3,4]y≡Xβ[φµν +−1ηµν ]αβyβX−1[φµν + ηµν ]αβ(13.6.18)αβWith this subtraction, the sum of the LP coefficients should be unity, up to roundoffand differences in how the φk ’s are estimated.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).void predic(float data[], int ndata, float d[], int m, float future[],int nfut)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 nfut data points, which it returns in the arrayfuture[1..nfut]. Note that the routine references only the last m values of data, as initialvalues for the prediction.{int k,j;float sum,discrp,*reg;13.6 Linear Prediction and Linear Predictive Coding571Linear Predictive Coding (LPC)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).A different, though related, method to which the formalism above can beapplied is the “compression” of a sampled signal so that it can be stored morecompactly. The original form should be exactly recoverable from the compressedversion. Obviously, compression can be accomplished only if there is redundancyin the signal.

Equation (13.6.11) describes one kind of redundancy: It says thatthe 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, or LPC.The basic idea of LPC (in its simplest form) is to record as a compressed file (i)the number of LP coefficients M , (ii) their M values, e.g., as obtained by memcof,(iii) the first M data points, and then (iv) for each subsequent data point only itsresidual discrepancy xi (equation 13.6.1). When you are creating the compressedfile, you find the residual by applying (13.6.1) to the previous M points, subtractingthe sum from the actual value of the current point.

When you are reconstructing theoriginal file, you add the residual back in, at the point indicated in the routine predic.It may not be obvious why there is any compression at all in this scheme. Afterall, we are storing one value of residual per data point! Why not just store the originaldata point? The answer depends on the relative sizes of the numbers involved.

Theresidual is obtained by subtracting two very nearly equal numbers (the data and thelinear prediction). Therefore, the discrepancy typically has only a very small numberof nonzero bits. These can be stored in a compressed file. How do you do it in ahigh-level language? Here is one way: Scale your data to have integer values, saybetween +1000000 and −1000000 (supposing that you need six significant figures).Modify equation (13.6.1) by enclosing the sum term in an “integer part of” operator.The discrepancy will now, by definition, be an integer. Experiment with differentvalues of M , to find LP coefficients that make the range of the discrepancy as smallas you can.

If you can get to within a range of ±127 (and in our experience this is notat all difficult) then you can write it to a file as a single byte. This is a compressionfactor of 4, compared to 4-byte integer or floating formats.Notice that the LP coefficients are computed using the quantized data, and thatthe 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 theinitial quantization of the data, you will not introduce even a single bit of roundofferror into the compression-reconstruction process: While the evaluation of the sumin (13.6.11) may have roundoff errors, the residual that you store is the value which,when added back to the sum, gives exactly the original (quantized) data value.

Noticealso that you do not need to massage the LP coefficients for stability; by adding theresidual back in to each point, you never depart from the original data, so instabilitiescannot grow. There is therefore no need for fixrts, above.Look at §20.4 to learn about Huffman coding, which will further compress theresiduals by taking advantage of the fact that smaller values of discrepancy will occurmore 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 occasional one isoutside, then reserve the value 127 to mean “out of range,” and then record on the file(immediately following the 127) a full-word value of the out-of-range discrepancy.§20.4 explains how to do much better.572Chapter 13.Fourier and Spectral ApplicationsCITED REFERENCES AND FURTHER READING:Childers, D.G. (ed.) 1978, Modern Spectrum Analysis (New York: IEEE Press), especially thepaper 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 Spatial Statistics and Digital Image Analysis (Washington: National AcademyPress). [3]Press, W.H., and Rybicki, G.B. 1992, Astrophysical Journal, vol. 398, pp. 169–176. [4]13.7 Power Spectrum Estimation by theMaximum Entropy (All Poles) MethodThe FFT is not the only way to estimate the power spectrum of a process, nor is itnecessarily 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 −fc < f < fc , but also the entire complex frequency plane. From thatvantage point, let us transform the complex f -plane to a new plane, called the z-transformplane or z-plane, by the relationz ≡ e2πif ∆(13.7.1)where ∆ is, as usual, the sampling interval in the time domain. Notice that the Nyquist intervalon the real axis of the f -plane maps one-to-one onto the unit circle in the complex z-plane.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).There are many variant procedures that all fall under the rubric of LPC.• If the spectral character of the data is time-variable, then it is best notto use a single set of LP coefficients for the whole data set, but ratherto 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 youcan tolerate some small amount of error, then don’t bother storing all of theresiduals.