f13-7

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 ) 566 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).This is an infinite Laurent series which depends on an infinite number of values ck. Equation (13.7.2) is just one kind of analytic approximation to the analytic function of zrepresented by (13.7.3); the kind, in fact, that is implicit in the use of FFTs to estimate power spectra byperiodogram methods. It goes under several names, including direct method, all-zero model, andmoving ave rage(MA) model . The term “all-zero” in particular refers to the fact that the model spectrum can have zeros in the z-plane, but not poles. If we look at the problem of approximating (13.7.3) more generally it seems clear that wecould do abetterjob witha rationalfunction, one withaseriesof type (13.7.2) inboth thenumerator and thedenominator. Lessobviously, itturnsout thattherearesome advantages inan approximation whose free parameters all lie in the denominator , namely, P(f)≈1 /vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleM/2/summationtext k=−M/2bkzk/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle2=a0/vextendsingle/vextendsingle/vextendsingle/vextendsingle1+ M/summationtext k=1akzk/vextendsingle/vextendsingle/vextendsingle/vextendsingle2(13.7.4 ) Here the second equality brings in a new set of coefficients ak’s, which can be determined from the bk’s using the fact that zlies on the unit circle. The bk’s can be thought of as being determined by the condition that power series expansion of (13.7.4) agree with thefirstM+1terms of (13.7.3). In practice, as we shall see, one determines the b k’s or ak’s by another method. The differences between the approximations (13.7.2) and (13.7.4) are not just cosmetic. They are approximations with very different character. Most notable is the fact that (13.7.4) can have poles, corresponding to infinite power spectral density, on the unit z-circle, i.e., at real frequencies in the Nyquist interval. Such poles can provide an accurate representationfor underlying power spectra that have sharp, discrete “lines” or delta-functions. By contrast,(13.7.2) can have only zeros, not poles, at real frequencies in the Nyquist interval, and mustthus attempt to fit sharp spectral features with, essentially, a polynomial. The approximation(13.7.4) goes under several names: all-poles model, maximum entropy method (MEM), autoregressive model (AR) .We need only find out how to compute the coefficients a 0and the ak’s from a data set, so that we can actually use (13.7.4) to obtain spectral estimates. A pleasant surprise is that we already know how! Look at equation (13.6.11) for linear prediction. Compare itwithlinearfilterequations (13.5.1)and (13.5.2),and you willsee that,viewed as a filter that takes input x’sinto output y’s,linear prediction has a filter function H(f)=1 1−N/summationtext j=1djz−j(13.7.5 ) Thus, the power spectrum of the y’s should be equal to the power spectrum of the x’s multiplied by |H(f)|2. Now let us think about what the spectrum of the input x’s is, when theyareresidualdiscrepancies fromlinearprediction. Althoughwewillnotproveitformally,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 would have 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 is exactly the quantity computed in equation (13.6.13) and returned by the routine memcofas xms. Thus, the coefficients a 0andakin equation (13.7.4) are related to the LP coefficients returned by memcofsimply by a0=xms ak=−d(k),k =1,...,M (13.7.6 ) Thereisalsoanotherwaytodescribetherelationbetweenthe ak’sandtheautocorrelation components φk. The Wiener-Khinchin theorem (12.0.12) says that the Fourier transform of the autocorrelation is equal to the power spectrum. In z-transform language, this Fourier transform is just a Laurent series in z. The equation that is to be satisfied by the coefficients in equation (13.7.4) is thus a0/vextendsingle/vextendsingle/vextendsingle/vextendsingle1+ M/summationtext k=1akzk/vextendsingle/vextendsingle/vextendsingle/vextendsingle2≈M/summationdisplay j=−Mφjzj(13.7.7 ) 13.7Maximum Entropy(AllPoles)Method 567Sample 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).The approximately equal sign in (13.7.7) has a somewhat special interpretation. It means that the series expansion of the left-hand side is supposed to agree with the right-hand sideterm by term from z −MtozM. Outside this range of terms, the right-hand side is obviously zero, while the left-hand side will still have nonzero terms. Notice that M, the number of coefficients in the approximation on the left-hand side, can be any integer up to N, the total number of autocorrelations available. (In practice, one often chooses Mmuch smaller than N.)Mis called the orderornumber of poles of the approximation. Whatever the chosen value of M, the series expansion of the left-hand side of (13.7.7) defines acertain sortof extrapolation ofthe autocorrelation function to lagslargerthan M,in fact even to lags larger than N, i.e.,larger than the run of data can actually measure . Itturns out that this particular extrapolation can be shown to have, among all possible extrapolations,the maximum entropyin a definable information-theoretic sense. Hence the name maximum entropy method , or MEM. The maximum entropy property has caused MEM to acquire a certain “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 of being able to fit sharp spectral features, but there is nothing else magical about its power spectrum estimates. The operations count in memcofscales as the product of N(the number of data points) andM(the desired order of the MEM approximation). If Mwere chosen to be as large as N, then the method would be much slower than the NlogNFFT methods of the previous section. In practice, however, one usually wants to limitthe 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 is not much slower than FFT estimation. We feel obliged to warn you that memcofcan be a bit quirky at times. If the number of poles 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 thesinewave.Also, with noisy input functions, if you choose too high an order, you will find spuriouspeaksgalore! Someexpertsrecommendtheuseofthisalgorithminconjunctionwithmoreconservative methods, likeperiodograms, tohelpchoose thecorrectmodelorder, andtoavoid getting too fooled by spurious spectral features. MEM can be finicky, but itcan also doremarkable things. Werecommend 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. There is no special significance to specific equally spaced frequencies as there was in the FFTcase. In fact, since the MEM estimate may have very sharp spectral features, one wantsto be able to evaluate it on a very fine mesh near to those features, but perhaps only morecoarsely farther away from them. Here is a subroutine which, given the coefficients alreadycomputed, evaluates (13.7.4) and returns the estimated power spectrum as a function of f∆ (the frequency times the sampling interval). Of course, f∆should lie in the Nyquist range between −1/2and1/2. FUNCTION evlmem(fdt,d,m,xms) INTEGER mREAL evlmem,fdt,xms,d(m) Given d, m, xms as returned by memcof , this function returns the power spectrum estimate P(f)as a function of fdt= f∆. INTEGER i REAL sumi,sumr DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp Trigonometric recurrences in double precision. theta=6.28318530717959d0*fdt wpr=cos(theta) Set up for recurrence relations. wpi=sin(theta) wr=1.d0wi=0.d0 sumr=1. These will accumulate the denominator of (13.7.4). 568 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).power spectral densitty 0.11101001000 .1 .15 .2 .25 .3 frequency f Figure 13.7.1. Sample output of maximum entropy spectral estimation. The input signal consists of 512 samples of the sum of two sinusoids of very nearly the same frequency, plus white noise with about equal power. Shown is an expanded portion of the full Nyquist frequency interval (which would extendfrom zero to 0.5). The dashed spectral estimate uses 20 poles; the dotted, 40; the solid, 150. With thelarger number of poles, the method can resolve the distinct sinusoids; but the flat noise background is beginning to show spurious peaks. (Note logarithmic scale.) sumi=0. do 11i=1,m Loop over the terms in the sum. wtemp=wrwr=wr*wpr-wi*wpiwi=wi*wpr+wtemp*wpi sumr=sumr-d(i)*sngl(wr) sumi=sumi-d(i)*sngl(wi) enddo 11 evlmem=xms/(sumr**2+sumi**2) Equation (13.7.4). returnEND Be sure to evaluate P(f)on afine enough grid to findany narrow features that may be 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 memcofandevlmemis normalized with respect to the mean square value of the input data vector. The answer is /integraldisplay 1/2 −1/2P(f∆)d(f∆) = 2/integraldisplay1/2 0P(f∆)d(f∆) =mean square value ofdata (13.7.8 ) Samplespectraproducedbytheroutines memcofandevlmemareshowninFigure13.7.1. 13.8SpectralAnalysisofUnevenlySampledData 569Sample 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).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 Sampled Data Thus 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 frequency fc≡1 2∆(13.8.2 ) as codified by the sampling theorem: A sampled data set like equation (13.8.1) contains complete information about all spectral components in a signal h(t)up to the Nyquist frequency, and scrambled or aliasedinformation about any signal components at frequencies largerthan theNyquist frequency. The sampling theorem thus de fines both the attractiveness, and the limitation, of any analysis of an evenly spaced data set. Therearesituations,however,whereevenly spaceddatacannotbeobtained. Acommon caseiswhereinstrumentaldrop-outs occur,sothatdataisobtained onlyon a(notconsecutiveinteger) subset of equation (13.8.1), the so-called missing data problem. Another case, common in observational sciences like astronomy, is that the observer cannot completelycontrol the time of the observations, but must simply accept a certain dictated set of t i’s. Therearesomeobvious waystogetfromunevenly spaced ti’stoevenly spaced ones,as inequation(13.8.1). Interpolationisoneway: laydownagridofevenlyspacedtimesonyourdataandinterpolatevaluesontothatgrid;thenuseFFTmethods. Inthemissingdataproblem,you only havetointerpolate on missingdata points. Ifa lotofconsecutive points aremissing,youmightaswelljustsetthemtozero,orperhaps “clamp”thevalueatthelastmeasuredpoint. However, the experience of practitioners of such interpolation techniques is not reassuring . Generally speaking, such techniques perform poorly. Long gaps in the data, for example,oftenproduceaspuriousbulgeofpoweratlowfrequencies( wavelengthscomparabletogaps). A completely different method of spectral analysis for unevenly sampled data, one that mitigates these dif ficulties and has some other very desirable properties, was developed by Lomb [1], based in part on earlier work by Barning [2]and Van´ıˇcek[3], and additionally elaborated by Scargle [4]. The Lomb method (as we will call it) evaluates data, and sines and cosines, only at times tithat are actually measured. Suppose that there are Ndata points hi≡h(ti),i=1,...,N. Thenfirstfind the mean and variance of the data by the usual formulas, h≡1 NN/summationdisplay 1hi σ2≡1 N−1N/summationdisplay 1(hi−h)2(13.8.3 ) Now, the Lomb normalized periodogram (spectral power as a function of angular frequency ω≡2πf > 0)i sd efined by PN(ω)≡1 2σ2  /bracketleftBig/summationtext j(hj−h)c o s ω(tj−τ)/bracketrightBig2 /summationtext jcos2ω(tj−τ)+/bracketleftBig/summationtext j(hj−h)s i nω(tj−τ)/bracketrightBig2 /summationtext jsin2ω(tj−τ)   (13.8.4 )