c13-9 (Numerical Recipes in C)

584Chapter 13.Fourier and Spectral Applications}CITED REFERENCES AND FURTHER READING:Lomb, N.R. 1976, Astrophysics and Space Science, vol. 39, pp. 447–462. [1]Barning, F.J.M. 1963, Bulletin of the Astronomical Institutes of the Netherlands, vol. 17, pp. 22–28. [2]Vanı́ček, P. 1971, Astrophysics and Space Science, vol. 12, pp. 10–33. [3]Scargle, J.D.

1982, Astrophysical Journal, vol. 263, pp. 835–853. [4]Horne, J.H., and Baliunas, S.L. 1986, Astrophysical Journal, vol. 302, pp. 757–763. [5]Press, W.H. and Rybicki, G.B. 1989, Astrophysical Journal, vol. 338, pp. 277–280. [6]13.9 Computing Fourier Integrals Using the FFTNot uncommonly, one wants to calculate accurate numerical values for integrals ofthe formZ bI=eiωt h(t)dt ,(13.9.1)aor the equivalent real and imaginary partsZ bIc =cos(ωt)h(t)dtaZIs =bsin(ωt)h(t)dt ,(13.9.2)aand one wants to evaluate this integral for many different values of ω.

In cases of interest, h(t)is often a smooth function, but it is not necessarily periodic in [a, b], nor does it necessarilygo to zero at a or b. While it seems intuitively obvious that the force majeure of the FFTought to be applicable to this problem, doing so turns out to be a surprisingly subtle matter,as we will now see.Let us first approach the problem naively, to see where the difficulty lies. Divide theinterval [a, b] into M subintervals, where M is a large integer, and defineb−a(13.9.3), tj ≡ a + j∆ , hj ≡ h(tj ) , j = 0, .

. . , MMNotice that h0 = h(a) and hM = h(b), and that there are M + 1 values hj . We canapproximate the integral I by a sum,∆≡I ≈∆M−1Xj=0hj exp(iωtj )(13.9.4)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).ix=(int)x;if (x == (float)ix) yy[ix] += y;else {ilo=LMIN(LMAX((long)(x-0.5*m+1.0),1),n-m+1);ihi=ilo+m-1;nden=nfac[m];fac=x-ilo;for (j=ilo+1;j<=ihi;j++) fac *= (x-j);yy[ihi] += y*fac/(nden*(x-ihi));for (j=ihi-1;j>=ilo;j--) {nden=(nden/(j+1-ilo))*(j-ihi);yy[j] += y*fac/(nden*(x-j));}}13.9 Computing Fourier Integrals Using the FFT585which is at any rate first-order accurate.

(If we centered the hj ’s and the tj ’s in the intervals,we could be accurate to second order.) Now for certain values of ω and M , the sum inequation (13.9.4) can be made into a discrete Fourier transform, or DFT, and evaluated bythe fast Fourier transform (FFT) algorithm. In particular, we can choose M to be an integerpower of 2, and define a set of special ω’s byI(ωm ) ≈ ∆eiωm aM−1Xhj e2πimj/M = ∆eiωm a [DFT(h0 . . . hM −1 )]m(13.9.6)j=0Equation (13.9.6), while simple and clear, is emphatically not recommended for use: It islikely to give wrong answers!The problem lies in the oscillatory nature of the integral (13.9.1). If h(t) is at all smooth,and if ω is large enough to imply several cycles in the interval [a, b] — in fact, ωm in equation(13.9.5) gives exactly m cycles — then the value of I is typically very small, so small thatit is easily swamped by first-order, or even (with centered values) second-order, truncationerror.

Furthermore, the characteristic “small parameter” that occurs in the error term is not∆/(b − a) = 1/M , as it would be if the integrand were not oscillatory, but ω∆, which can beas large as π for ω’s within the Nyquist interval of the DFT (cf. equation 13.9.5). The resultis that equation (13.9.6) becomes systematically inaccurate as ω increases.It is a sobering exercise to implement equation (13.9.6) for an integral that can be doneanalytically, and to see just how bad it is. We recommend that you try it.Let us therefore turn to a more sophisticated treatment.

Given the sampled points hj , wecan approximate the function h(t) everywhere in the interval [a, b] by interpolation on nearbyhj ’s. The simplest case is linear interpolation, using the two nearest hj ’s, one to the left andone to the right. A higher-order interpolation, e.g., would be cubic interpolation, using twopoints to the left and two to the right — except in the first and last subintervals, where wemust interpolate with three hj ’s on one side, one on the other.The formulas for such interpolation schemes are (piecewise) polynomial in the independent variable t, but with coefficients that are of course linear in the function valueshj .

Although one does not usually think of it in this way, interpolation can be viewed asapproximating a function by a sum of kernel functions (which depend only on the interpolationscheme) times sample values (which depend only on the function). Let us writeMXXt − tjt − tjh(t) ≈hj ψhj ϕj+(13.9.7)∆∆j=0j=endpointsHere ψ(s) is the kernel function of an interior point: It is zero for s sufficiently negativeor sufficiently positive, and becomes nonzero only when s is in the range where thehj multiplying it is actually used in the interpolation.

We always have ψ(0) = 1 andψ(m) = 0, m = ±1, ±2, . . . , since interpolation right on a sample point should give thesampled function value. For linear interpolation ψ(s) is piecewise linear, rises from 0 to 1for s in (−1, 0), and falls back to 0 for s in (0, 1). For higher-order interpolation, ψ(s) ismade up piecewise of segments of Lagrange interpolation polynomials. It has discontinuousderivatives at integer values of s, where the pieces join, because the set of points used inthe interpolation changes discretely.As already remarked, the subintervals closest to a and b require different (noncentered)interpolation formulas. This is reflected in equation (13.9.7) by the second sum, with thespecial endpoint kernels ϕj (s).

Actually, for reasons that will become clearer below, we haveincluded all the points in the first sum (with kernel ψ), so the ϕj ’s are actually differencesbetween true endpoint kernels and the interior kernel ψ. It is a tedious, but straightforward,exercise to write down all the ϕj (s)’s for any particular order of interpolation, each oneconsisting of differences of Lagrange interpolating polynomials spliced together piecewise.RbNow apply the integral operator a dt exp(iωt) to both sides of equation (13.9.7),interchange the sums and integral, and make the changes of variable s = (t − tj )/∆ in theSample 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).2πm(13.9.5)Mwhere m has the values m = 0, 1, .

. . , M/2 − 1. Then equation (13.9.4) becomesωm ∆ ≡586Chapter 13.Fourier and Spectral Applicationsfirst sum, s = (t − a)/∆ in the second sum. The result is"MXXhj eijθ +I ≈ ∆eiωa W (θ)j=0#hj αj (θ)(13.9.8)j=endpoints(13.9.9)(13.9.10)−∞The key point is that equations (13.9.9) and (13.9.10) can be evaluated, analytically,once and for all, for any given interpolation scheme. Then equation (13.9.8) is an algorithmfor applying “endpoint corrections” to a sum which (as we will see) can be done using theFFT, giving a result with high-order accuracy.We will consider only interpolations that are left-right symmetric.

Then symmetryimpliesϕM −j (s) = ϕj (−s)αM −j (θ) = eiθM α*j (θ) = eiω(b−a) α*j (θ)(13.9.11)where * denotes complex conjugation. Also, ψ(s) = ψ(−s) implies that W (θ) is real.Turn now to the first sum in equation (13.9.8), which we want to do by FFT methods.To do so, choose some N that is an integer power of 2 with N ≥ M + 1. (Note thatM need not be a power of two, so M = N − 1 is allowed.) If N > M + 1, definehj ≡ 0, M + 1 < j ≤ N − 1, i.e., “zero pad” the array of hj ’s so that j takes on the range0 ≤ j ≤ N − 1. Then the sum can be done as a DFT for the special values ω = ωn given by2πnNωn ∆ ≡≡θn = 0, 1, . . . ,−1(13.9.12)N2For fixed M , the larger N is chosen, the finer the sampling in frequency space.

Thevalue M , on the other hand, determines the highest frequency sampled, since ∆ decreaseswith increasing M (equation 13.9.3), and the largest value of ω∆ is always just under π(equation 13.9.12). In general it is advantageous to oversample by at least a factor of 4, i.e.,N > 4M (see below). We can now rewrite equation (13.9.8) in its final form asiωn aI(ωn ) = ∆eW (θ)[DFT(h0 .

. . hN −1 )]n+ α0 (θ)h0 + α1 (θ)h1 + α2 (θ)h2 + α3(θ)h3 + . . .hi+ eiω(b−a) α*0 (θ)hM + α*1 (θ)hM −1 + α*2 (θ)hM −2 + α*3 (θ)hM −3 + . . .(13.9.13)For cubic (or lower) polynomial interpolation, at most the terms explicitly shown aboveare nonzero; the ellipses (. . .) can therefore be ignored, and we need explicit forms only forthe functions W, α0 , α1 , α2 , α3 , calculated with equations (13.9.9) and (13.9.10). We haveworked these out for you, in the trapezoidal (second-order) and cubic (fourth-order) cases.Here are the results, along with the first few terms of their power series expansions for small θ:Trapezoidal order:W (θ) =1 21 412(1 − cos θ)θ +θ −θ6≈1−θ21236020160(1 − cos θ)(θ − sin θ)+iθ2θ21 21 211 411 411θ −θ +θ6 + iθ−θ +θ −θ6≈− +2247204032061205040362880α0 (θ) = −α1 = α2 = α3 = 0Sample 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.