572 Chapter 13.Fourier and Spectral Applications 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 not to use a single set of LP coefficients for the whole data set,but rather to partition the data into segments,computing and storing different LP coefficients for each segment. If the data are really well characterized by their LP coefficients,and you can tolerate some small amount oferror,then don't bother storing all of the residuals.Just do linear prediction until you are outside of tolerances,then reinitialize(using M sequential stored residuals)and continue predicting. In some applications,most notably speech synthesis,one cares only about the spectral content of the reconstructed signal,not the relative phases. In this case,one need not store any starting values at all,but only the LP coefficients for each segment of the data.The output is reconstructed by driving these coefficients with initial conditions consisting of all zeros except for one nonzero spike.A speech synthesizer chip may have of order 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,even when the residuals zi are not small.The xi's are then interpreted as the underlying"input signal"which,when filtered through the all-poles filter defined by the LP coefficients(see $13.7),produces the observed "output University Press. signal."LPC reveals simultaneously,it is said,the nature of the filter and the particular input that is driving it.We are skeptical of these applications, the literature,however,is full of extravagant claims. Programs CITED REFERENCES AND FURTHER READING: OF SCIENTIFIC Childers,D.G.(ed.)1978,Modem 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 Spatia/Statistics and Digital Image Analysis(Washington:National Academy Press).[3] Press,W.H.,and Rybicki,G.B.1992,Astrophysica/Journal,vol.398,pp.169-176.[4] Numerical Recipes 10621 43106 13.7 Power Spectrum Estimation by the (outside Maximum Entropy (All Poles)Method Software. 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 the Nyquist interval-fe<f<fe,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 plane or z-plane,by the relation 2三e2mfA (13.7.1) where A is,as usual,the sampling interval in the time domain.Notice that the Nyquist interval on the real axis of the f-plane maps one-to-one onto the unit circle in the complex z-plane
572 Chapter 13. Fourier and Spectral Applications 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 machinereadable 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). 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 not to use a single set of LP coefficients for the whole data set, but rather to partition the data into segments, computing and storing different LP coefficients for each segment. • If the data are really well characterized by their LP coefficients, and you can tolerate some small amount of error, then don’t bother storing all of the residuals. Just do linear prediction until you are outside of tolerances, then reinitialize (using M sequential stored residuals) and continue predicting. • In some applications, most notably speech synthesis, one cares only about the spectral content of the reconstructed signal, not the relative phases. In this case, one need not store any starting values at all, but only the LP coefficients for each segment of the data. The output is reconstructed by driving these coefficients with initial conditions consisting of all zeros except for one nonzero spike. A speech synthesizer chip may have of order 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, even when the residuals xi are not small. The xi’s are then interpreted as the underlying “input signal” which, when filtered through the 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 the particular input that is driving it. We are skeptical of these applications; 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 Spatial Statistics and Digital Image Analysis (Washington: National Academy 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 the Nyquist interval −fc <f<fc, 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 plane or z-plane, by the relation z ≡ e2πif∆ (13.7.1) where ∆ is, as usual, the sampling interval in the time domain. Notice that the Nyquist interval on the real axis of the f-plane maps one-to-one onto the unit circle in the complex z-plane
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,when
13.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 machinereadable 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
574 Chapter 13.Fourier and Spectral Applications 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 a's would have allowed a more accurate linear prediction,and would have been removed.)The overall normalization of this flat spectrum is just the mean square amplitude of the r's.But this is exactly the quantity computed in equation (13.6.13)and returned by the routine memcof as xms.Thus,the coefficients ao and ak in equation (13.7.4)are related to the LP coefficients returned by memcof simply by o=xms ak=-d(k),k=1,..,M (13.7.6) There is also another way to describe the relation between the ak's and the autocorrelation components o.The Wiener-Khinchin theorem (12.0.12)says that the Fourier transform of 8 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 ao (13.7.7) 3 1+ k三1 RECIPES 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 side term by term from M to zM.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 子g Press. number of autocorrelations available.(In practice,one often chooses M much smaller than N.)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) Programs defines a certain sort of extrapolation of the autocorrelation function to lags larger than M,in fact even to lags larger than N,i.e.,larger than the run of data can actually measure.It turns SCIENTIFIC out 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 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"estimate 排 than 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 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 as 10621 N,then the method would be much slower than the Nlog N FFT methods of the previous Numerica section.In practice,however,one usually wants to limit the order (or number of poles)of the uctio Recipes 43106 MEM approximation to a few times the number of sharp spectral features that one desires it to 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,one (outside might take M=10 or 20 or 50 for N =1000 or 10000.In that case MEM estimation is not much slower than FFT estimation. North Software. We feel obliged to warn you that memcof can 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 double precision.With"peaky"data (i.e.,data with extremely sharp spectral features),the algorithm may suggest split peaks even at modest orders,and the peaks may shift with the phase of the sine wave.Also,with noisy input functions,if you choose too high an order,you will find spurious peaks galore!Some experts recommend the use of this algorithm in conjunction with more conservative methods,like periodograms,to help choose the correct model order,and to avoid getting too fooled by spurious spectral features.MEM can be finicky,but it can also do remarkable things.We recommend that you try it out,cautiously,on your own problems.We now 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 FFT case
574 Chapter 13. Fourier and Spectral Applications 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 machinereadable 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). 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 would have allowed a more accurate linear prediction, and would have been removed.) The overall normalization 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 memcof as xms. Thus, the coefficients a0 and ak in equation (13.7.4) are related to the LP coefficients returned by memcof simply by a0 = xms ak = −d(k), k = 1,...,M (13.7.6) There is also another way to describe the relation between the ak’s and the autocorrelation 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 1 + M k=1 akzk 2 ≈ M j=−M φj zj (13.7.7) 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 side term by term from z−M to zM. 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 M much smaller than N.) 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, in fact even to lags larger than N, i.e., larger than the run of data can actually measure. It turns out 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 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” estimate than 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 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 as N, then the method would be much slower than the N log N FFT methods of the previous section. In practice, however, one usually wants to limit the order (or number of poles) of the MEM approximation to a few times the number of sharp spectral features that one desires it to 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, one might 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 memcof can 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 double precision. With “peaky” data (i.e., data with extremely sharp spectral features), the algorithm may suggest split peaks even at modest orders, and the peaks may shift with the phase of the sine wave. Also, with noisy input functions, if you choose too high an order, you will find spurious peaks galore! Some experts recommend the use of this algorithm in conjunction with more conservative methods, like periodograms, to help choose the correct model order, and to avoid getting too fooled by spurious spectral features. MEM can be finicky, but it can also do remarkable things. We recommend that you try it out, cautiously, on your own problems. We now 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 FFT case.
13.8 Spectral Analysis of Unevenly Sampled Data 575 In fact,since the MEM estimate may have very sharp spectral features,one wants to be able to evaluate it on a very fine mesh near to those features,but perhaps only more coarsely farther away 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 fA(the frequency times the sampling interval).Of course,fA should lie in the Nyquist range between-1/2 and 1/2. #include 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 spectrum estimate P(f)as a function of fdt =fA. 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. nted for thetas=6.28318530717959*fdt wpr=cos(theta); Set up for recurrence relations. wpi=sin(theta); for(i=1;1<=m;i+){ Loop over the terms in the sum wr=(wtemp=wr)*wpr-wi*wpi; RECIPES wi=wi*wpr+wtemp*wpi; sumr -d[i]*wr; These accumulate the denominator of(13.7.4). sumi -d[i]*wi; return xms/(sumr*sumr+sumi*sumi); Equation (13.7.4). 子二 Press. Be sure to evaluate P(f)on a fine enough grid to find any 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 memcof and evlmem is normalized with respect to the mean square value of the input data vector.The answer is SCIENTIFIC r/ r1/2 6 P(fA)d(fA)=2 P(f△)d(f△)=mean square value of data (13.7.8) -1/2 n Sample 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 Il. Numerica 10621 Kay.S.M.,and Marple,S.L.1981,Proceedings of the /EEE,vol.69,pp.1380-1419. uctio Recipes 43108 (outside Software. 13.8 Spectral Analysis of Unevenly Sampled ying of Data Thus far,we have been dealing exclusively with evenly sampled data, hm=h(n△)n=.,-3,-2,-1,0,1,2,3,. (13.8.1) where A is the sampling interval,whose reciprocal is the sampling rate.Recall also (12.1) the significance of the Nyquist critical frequency 1 f:=2公 (13.8.2)
13.8 Spectral Analysis of Unevenly Sampled Data 575 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 machinereadable 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). In fact, since the MEM estimate may have very sharp spectral features, one wants to be able to evaluate it on a very fine mesh near to those features, but perhaps only more coarsely farther away 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 times the sampling interval). Of course, f∆ should lie in the Nyquist range between −1/2 and 1/2. #include 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 spectrum estimate 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); Set up for recurrence relations. wpi=sin(theta); for (i=1;i<=m;i++) { Loop over the terms in the sum. wr=(wtemp=wr)*wpr-wi*wpi; wi=wi*wpr+wtemp*wpi; sumr -= d[i]*wr; These accumulate the denominator of (13.7.4). sumi -= d[i]*wi; } return xms/(sumr*sumr+sumi*sumi); Equation (13.7.4). } Be sure to evaluate P(f) on a fine enough grid to find any 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 memcof and evlmem is normalized with respect to the mean square value of the input data vector. The answer is 1/2 −1/2 P(f∆)d(f∆) = 2 1/2 0 P(f∆)d(f∆) = mean square value of data (13.7.8) Sample 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 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)