13.7 Maximum Entropy(All Poles)Method 573 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 c=c(t)can be written, except for normalization convention,as 2 P( Ck (13.7.2) k=-N/2 Of course,(13.7.2)is not the true power spectrum of the underlying function c(t),but only an estimate.We can see in two related ways why the estimate is not likely to be exact.First.in the time domain,the estimate is based on only a finite range of the function c(t)which may,for all we know,have continued from t=-oo to oo.Second,in the z-plane of equation (13.7.2),the 8 finite Laurent series offers,in general,only an approximation to a general analytic function of z.In fact,a formal expression for representing "true"power spectra (up to normalization)is PO) (13.7.3) This is an infinite Laurent series which depends on an infinite number of values c.Equation (13.7.2)is just one kind of analytic approximation to the analytic function of z represented RECIPES by (13.7.3):the kind,in fact,that is implicit in the use of FFTs to estimate power spectra by periodogram 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 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 入@ Press. we could do a better job with a rational function,one with a series of type(13.7.2)in both the numerator and the denominator.Less obviously,it turns out that there are some advantages in an approximation whose free parameters all lie in the denominator,namely, 1 P(f) ao (13.7.4) SCIENTIFIC M/2 6 k=-M/2 Here the second equality brings in a new set of coefficients ak's,which can be determined from the b's using the fact that z lies on the unit circle.The be's can be thought of as being determined by the condition that power series expansion of(13.7.4)agree with the first M+1 terms of (13.7.3).In practice,as we shall see,one determines the b's or ak's by another method. The differences between the approximations (13.7.2)and(13.7.4)are not just cosmetic. Numerica 10621 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 uction Recipes 43108 real frequencies in the Nyquist interval.Such poles can provide an accurate representation for underlying power spectra that have sharp,discrete"lines"or delta-functions.By contrast, (outs (13.7.2)can have only zeros,not poles,at real frequencies in the Nyquist interval,and must thus 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, North Software. autoregressive model (AR).We need only find out how to compute the coefficients ao and the a'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 it with linear filter equations (13.5.1)and (13.5.2),and you will see that, viewed as a filter that takes input r's into output y's,linear prediction has a filter function H()= 1 (13.7.5) 1- djz-3 7=1 Thus,the power spectrum of the y's should be equal to the power spectrum of the z's multiplied by (f).Now let us think about what the spectrum of the input z's is,when13.7 Maximum Entropy (All Poles) Method 573 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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) =





N/ 2−1 k=−N/2 ckzk





2 (13.7.2) Of course, (13.7.2) is not the true power spectrum of the underlying function c(t), but only an estimate. We can see in two related ways why the estimate is not likely to be exact. First, in the time domain, the estimate is based on only a finite range of the function c(t) which may, for all we know, have continued from t = −∞ to ∞. Second, in the z-plane of equation (13.7.2), the finite Laurent series offers, in general, only an approximation to a general analytic function of z. In fact, a formal expression for representing “true” power spectra (up to normalization) is P(f) =




∞ k=−∞ ckzk




2 (13.7.3) 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 represented by (13.7.3); the kind, in fact, that is implicit in the use of FFTs to estimate power spectra by periodogram 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 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 we could do a better job with a rational function, one with a series of type (13.7.2) in both the numerator and the denominator. Less obviously, it turns out that there are some advantages in an approximation whose free parameters all lie in the denominator, namely, P(f) ≈ 1




M/ 2 k=−M/2 bkzk




2 = a0



1 + M k=1 akzk



2 (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 z lies 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 the first M + 1 terms of (13.7.3). In practice, as we shall see, one determines the bk’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 representation for 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 must thus 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 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 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 function H(f) = 1 1 − N j=1 dj z−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