c13-7 (Numerical Recipes in C)

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. Just do linear prediction until you are outside of tolerances, thenreinitialize (using M sequential stored residuals) and continue predicting.• In some applications, most notably speech synthesis, one cares only aboutthe spectral content of the reconstructed signal, not the relative phases.In this case, one need not store any starting values at all, but only theLP coefficients for each segment of the data.

The output is reconstructedby driving these coefficients with initial conditions consisting of all zerosexcept for one nonzero spike. A speech synthesizer chip may have oforder 10 LP coefficients, which change perhaps 20 to 50 times per second.• Some people believe that it is interesting to analyze a signal by LPC, evenwhen the residuals xi are not small. The xi ’s are then interpreted as theunderlying “input signal” which, when filtered through the all-poles filterdefined by the LP coefficients (see §13.7), produces the observed “outputsignal.” LPC reveals simultaneously, it is said, the nature of the filter andthe particular input that is driving it. We are skeptical of these applications;the literature, however, is full of extravagant claims.13.7 Maximum Entropy (All Poles) Method573If we now compare (13.7.1) to equations (13.4.4) and (13.4.6), we see that the FFTpower spectrum estimate (13.4.5) for any real sampled function ck ≡ c(tk ) can be written,except for normalization convention, as2 N/2−1 XP (f ) = ck z k (13.7.2)k=−N/2k=−∞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 z representedby (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,and moving average (MA) model. The term “all-zero” in particular refers to the fact that themodel 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 thatwe could do a better job with a rational function, one with a series of type (13.7.2) in both thenumerator and the denominator. Less obviously, it turns out that there are some advantages inan approximation whose free parameters all lie in the denominator, namely,1a0P (f ) ≈ 2 = 2M M/2PPkak z bk z k 1 +k=1k=−M/2(13.7.4)Here the second equality brings in a new set of coefficients ak ’s, which can be determinedfrom the bk ’s using the fact that z lies on the unit circle.

The bk ’s can be thought of asbeing determined by the condition that power series expansion of (13.7.4) agree with thefirst M + 1 terms of (13.7.3). In practice, as we shall see, one determines the bk ’s orak ’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., atreal 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 a0 and theak ’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 linearprediction. Compare it with linear filter equations (13.5.1) and (13.5.2), and you will see that,viewed as a filter that takes input x’s into output y’s, linear prediction has a filter functionH(f ) =1−1NP(13.7.5)dj z jj=1Thus, the power spectrum of the y’s should be equal to the power spectrum of the x’smultiplied by |H(f )|2.

Now let us think about what the spectrum of the input x’s is, whenSample 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).Of course, (13.7.2) is not the true power spectrum of the underlying function c(t), but only anestimate. We can see in two related ways why the estimate is not likely to be exact.

First, in thetime domain, the estimate is based on only a finite range of the function c(t) which may, for allwe know, have continued from t = −∞ to ∞. Second, in the z-plane of equation (13.7.2), thefinite Laurent series offers, in general, only an approximation to a general analytic function ofz. In fact, a formal expression for representing “true” power spectra (up to normalization) is2∞ XkP (f ) = ck z (13.7.3)574Chapter 13.Fourier and Spectral Applicationsa0 = xmsak = −d(k),k = 1, .

. . , M(13.7.6)There is also another way to describe the relation between the ak ’s and the autocorrelationcomponents φk . The Wiener-Khinchin theorem (12.0.12) says that the Fourier transform ofthe autocorrelation is equal to the power spectrum. In z-transform language, this Fouriertransform is just a Laurent series in z. The equation that is to be satisfied by the coefficientsin equation (13.7.4) is thusMXa0φj z j2 ≈MPj=−Mk1 +ak z (13.7.7)k=1The approximately equal sign in (13.7.7) has a somewhat special interpretation. It meansthat the series expansion of the left-hand side is supposed to agree with the right-hand sideterm by term from z −M to z M . Outside this range of terms, the right-hand side is obviouslyzero, while the left-hand side will still have nonzero terms. Notice that M , the number ofcoefficients in the approximation on the left-hand side, can be any integer up to N , the totalnumber of autocorrelations available.