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

(In practice, one often chooses M much smaller thanN .) M is called the order or number of poles of the approximation.Whatever the chosen value of M , the series expansion of the left-hand side of (13.7.7)defines a certain sort of extrapolation of the autocorrelation function to lags larger than M , infact even to lags larger than N , i.e., larger than the run of data can actually measure. It turnsout that this particular extrapolation can be shown to have, among all possible extrapolations,the maximum entropy in a definable information-theoretic sense. Hence the name maximumentropy method, or MEM. The maximum entropy property has caused MEM to acquire acertain “cult” popularity; one sometimes hears that it gives an intrinsically “better” estimatethan is given by other methods.

Don’t believe it. MEM has the very cute property ofbeing able to fit sharp spectral features, but there is nothing else magical about its powerspectrum estimates.The operations count in memcof scales as the product of N (the number of data points)and M (the desired order of the MEM approximation). If M were chosen to be as large asN , then the method would be much slower than the N log N FFT methods of the previoussection. In practice, however, one usually wants to limit the order (or number of poles) of theMEM approximation to a few times the number of sharp spectral features that one desires itto fit. With this restricted number of poles, the method will smooth the spectrum somewhat,but this is often a desirable property.

While exact values depend on the application, onemight take M = 10 or 20 or 50 for N = 1000 or 10000. In that case MEM estimation isnot much slower than FFT estimation.We feel obliged to warn you that memcof can be a bit quirky at times. If the number ofpoles or number of data points is too large, roundoff error can be a problem, even in doubleprecision. With “peaky” data (i.e., data with extremely sharp spectral features), the algorithmmay suggest split peaks even at modest orders, and the peaks may shift with the phase of thesine wave.

Also, with noisy input functions, if you choose too high an order, you will findspurious peaks galore! Some experts recommend the use of this algorithm in conjunction withmore conservative methods, like periodograms, to help choose the correct model order, and toavoid getting too fooled by spurious spectral features. MEM can be finicky, but it can also doremarkable things. We recommend that you try it out, cautiously, on your own problems. Wenow turn to the evaluation of the MEM spectral estimate from its coefficients.The MEM estimation (13.7.4) is a function of continuously varying frequency f .

Thereis no special significance to specific equally spaced frequencies as there was in the FFT case.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).they are residual discrepancies from linear prediction.

Although we will not prove it formally,it is intuitively believable that the x’s are independently random and therefore have a flat(white noise) spectrum. (Roughly speaking, any residual correlations left in the x’s wouldhave allowed a more accurate linear prediction, and would have been removed.) The overallnormalization of this flat spectrum is just the mean square amplitude of the x’s. But this isexactly the quantity computed in equation (13.6.13) and returned by the routine memcof asxms. Thus, the coefficients a0 and ak in equation (13.7.4) are related to the LP coefficientsreturned by memcof simply by57513.8 Spectral Analysis of Unevenly Sampled DataIn fact, since the MEM estimate may have very sharp spectral features, one wants to be able toevaluate it on a very fine mesh near to those features, but perhaps only more coarsely fartheraway from them.

Here is a function which, given the coefficients already computed, evaluates(13.7.4) and returns the estimated power spectrum as a function of f ∆ (the frequency timesthe sampling interval). Of course, f ∆ should lie in the Nyquist range between −1/2 and 1/2.float evlmem(float fdt, float d[], int m, float xms)Given d[1..m], m, xms as returned by memcof, this function returns the power spectrumestimate P (f ) as a function of fdt = f ∆.{int i;float sumr=1.0,sumi=0.0;double wr=1.0,wi=0.0,wpr,wpi,wtemp,theta;Trig.

recurrences in double precision.theta=6.28318530717959*fdt;wpr=cos(theta);wpi=sin(theta);for (i=1;i<=m;i++) {wr=(wtemp=wr)*wpr-wi*wpi;wi=wi*wpr+wtemp*wpi;sumr -= d[i]*wr;sumi -= d[i]*wi;}return xms/(sumr*sumr+sumi*sumi);Set up for recurrence relations.Loop over the terms in the sum.These accumulate the denominator of (13.7.4).Equation (13.7.4).}Be sure to evaluate P (f ) on a fine enough grid to find any narrow features that maybe there! Such narrow features, if present, can contain virtually all of the power in the data.You might also wish to know how the P (f ) produced by the routines memcof and evlmem isnormalized with respect to the mean square value of the input data vector.

The answer isZ 1/2Z 1/2P (f ∆)d(f ∆) = 2P (f ∆)d(f ∆) = mean square value of data(13.7.8)−1/20Sample spectra produced by the routines memcof and evlmem are shown in Figure 13.7.1.CITED REFERENCES AND FURTHER READING:Childers, D.G. (ed.) 1978, Modern Spectrum Analysis (New York: IEEE Press), Chapter II.Kay, S.M., and Marple, S.L. 1981, Proceedings of the IEEE, vol. 69, pp. 1380–1419.13.8 Spectral Analysis of Unevenly SampledDataThus far, we have been dealing exclusively with evenly sampled data,hn = h(n∆)n = . . .

, −3, −2, −1, 0, 1, 2, 3, . . .(13.8.1)where ∆ is the sampling interval, whose reciprocal is the sampling rate. Recall also (§12.1)the significance of the Nyquist critical frequencyfc ≡12∆(13.8.2)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).#include <math.h>.