c13-4 (Numerical Recipes in C)

54913.4 Power Spectrum Estimation Using the FFT N 2 (extrapolated) S 2 (deduced)fFigure 13.3.1. Optimal (Wiener) filtering. The power spectrum of signal plus noise shows a signal peakadded to a noise tail. The tail is extrapolated back into the signal region as a “noise model.” Subtractinggives the “signal model.” The models need not be accurate for the method to be useful. A simplealgebraic combination of the models gives the optimal filter (see text).new signal which you could improve even further with the same filtering technique.Don’t waste your time on this line of thought. The scheme converges to a signal ofS(f) = 0.

Converging iterative methods do exist; this just isn’t one of them.You can use the routine four1 (§12.2) or realft (§12.3) to FFT your datawhen you are constructing an optimal filter. To apply the filter to your data, youcan use the methods described in §13.1. The specific routine convlv is not neededfor optimal filtering, since your filter is constructed in the frequency domain tobegin with. If you are also deconvolving your data with a known response function,however, you can modify convlv to multiply by your optimal filter just before ittakes the inverse Fourier transform.CITED REFERENCES AND FURTHER READING:Rabiner, L.R., and Gold, B.

1975, Theory and Application of Digital Signal Processing (EnglewoodCliffs, NJ: Prentice-Hall).Nussbaumer, H.J. 1982, Fast Fourier Transform and Convolution Algorithms (New York: SpringerVerlag).Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (NewYork: Academic Press).13.4 Power Spectrum Estimation Using the FFTIn the previous section we “informally” estimated the power spectral density of afunction c(t) by taking the modulus-squared of the discrete Fourier transform of someSample 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).log scale C 2 (measured)550Chapter 13.Fourier and Spectral ApplicationsN−1X2|cj | ≡ “sum squared amplitude”(13.4.1)j=01TZTZT2|c(t)| dt ≈02|c(t)| dt ≈ ∆0N−11 X2|cj | ≡ “mean squared amplitude”N j=0N−1X2|cj | ≡ “time-integral squared amplitude”(13.4.2)(13.4.3)j=0PSD estimators, as we shall see, have an even greater variety. In this section,we consider a class of them that give estimates at discrete values of frequency fi ,where i will range over integer values.

In the next section, we will learn abouta different class of estimators that produce estimates that are continuous functionsof frequency f. Even if it is agreed always to relate the PSD normalization to aparticular description of the function normalization (e.g., 13.4.2), there are at leastthe following possibilities: The PSD is• defined for discrete positive, zero, and negative frequencies, and its sumover these is the function mean squared amplitude• defined for zero and discrete positive frequencies only, and its sum overthese is the function mean squared amplitude• defined in the Nyquist interval from −fc to fc , and its integral over thisrange is the function mean squared amplitude• defined from 0 to fc , and its integral over this range is the function meansquared amplitudeIt never makes sense to integrate the PSD of a sampled function outside of theNyquist interval −fc and fc since, according to the sampling theorem, power therewill have been aliased into the Nyquist interval.It is hopeless to define enough notation to distinguish all possible combinationsof normalizations.

In what follows, we use the notation P (f) to mean any of theabove PSDs, stating in each instance how the particular P (f) is normalized. Bewarethe inconsistent notation in the literature.The method of power spectrum estimation used in the previous section is asimple version of an estimator called, historically, the periodogram. If we take anN -point sample of the function c(t) at equal intervals and use the FFT to computeSample 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).finite, sampled stretch of it.

In this section we’ll do roughly the same thing, but withconsiderably greater attention to details. Our attention will uncover some surprises.The first detail is power spectrum (also called a power spectral density orPSD) normalization. In general there is some relation of proportionality between ameasure of the squared amplitude of the function and a measure of the amplitudeof the PSD. Unfortunately there are several different conventions for describingthe normalization in each domain, and many opportunities for getting wrong therelationship between the two domains. Suppose that our function c(t) is sampled atN points to produce values c0 . . . cN−1 , and that these points span a range of timeT , that is T = (N − 1)∆, where ∆ is the sampling interval. Then here are severaldifferent descriptions of the total power:55113.4 Power Spectrum Estimation Using the FFTits discrete Fourier transformCk =N−1Xcj e2πijk/Nk = 0, .

. . , N − 1(13.4.4)j=012P (0) = P (f0 ) = 2 |C0 |Ni1 h22P (fk ) = 2 |Ck | + |CN−k |N21 P (fc ) = P (fN/2 ) = 2 CN/2 Nk = 1, 2, . . . ,N−12(13.4.5)where fk is defined only for the zero and positive frequenciesfk ≡kk= 2fcN∆Nk = 0, 1, . . . ,N2(13.4.6)By Parseval’s theorem, equation (12.1.10), we see immediately that equation (13.4.5)is normalized so that the sum of the N/2 + 1 values of P is equal to the meansquared amplitude of the function cj .We must now ask this question. In what sense is the periodogram estimate(13.4.5) a “true” estimator of the power spectrum of the underlying function c(t)?You can find the answer treated in considerable detail in the literature cited (see,e.g., [1] for an introduction). Here is a summary.First, is the expectation value of the periodogram estimate equal to the powerspectrum, i.e., is the estimator correct on average? Well, yes and no.

We wouldn’treally expect one of the P (fk )’s to equal the continuous P (f) at exactly fk , since fkis supposed to be representative of a whole frequency “bin” extending from halfwayfrom the preceding discrete frequency to halfway to the next one. We should beexpecting the P (fk ) to be some kind of average of P (f) over a narrow windowfunction centered on its fk . For the periodogram estimate (13.4.6) that windowfunction, as a function of s the frequency offset in bins, is2sin(πs)1W (s) = 2Nsin(πs/N )(13.4.7)Notice that W (s) has oscillatory lobes but, apart from these, falls off only about asW (s) ≈ (πs)−2 . This is not a very rapid fall-off, and it results in significant leakage(that is the technical term) from one frequency to another in the periodogram estimate.Notice also that W (s) happens to be zero for s equal to a nonzero integer.

This meansthat if the function c(t) is a pure sine wave of frequency exactly equal to one of thefk ’s, then there will be no leakage to adjacent fk ’s. But this is not the characteristiccase! If the frequency is, say, one-third of the way between two adjacent fk ’s, thenthe leakage will extend well beyond those two adjacent bins. The solution to theproblem of leakage is called data windowing, and we will discuss it below.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).then the periodogram estimate of the power spectrum is defined at N/2 + 1frequencies as552Chapter 13.Fourier and Spectral ApplicationsSample 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).Turn now to another question about the periodogram estimate. What is thevariance of that estimate as N goes to infinity? In other words, as we take moresampled points from the original function (either sampling a longer stretch of data atthe same sampling rate, or else by resampling the same stretch of data with a fastersampling rate), then how much more accurate do the estimates Pk become? Theunpleasant answer is that the periodogram estimates do not become more accurateat all! In fact, the variance of the periodogram estimate at a frequency fk is alwaysequal to the square of its expectation value at that frequency.