584 Chapter 13.Fourier and Spectral Applications ix=(int)x; if (x ==(float)ix)yy[ix]+=y; else 11o=LMIN(LMAX((1ong)(x-0.5*m+1.0),1),n-m+1); 1h1=11o+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)); http://www.nr. readable files 83g 11-800 (including this one) granted for from NUMERICAL RECIPES IN 19881992 CITED REFERENCES AND FURTHER READING: Lomb,N.R.1976,Astrophysics and Space Science,vol.39,pp.447-462.[1] /Cambridge Barning,F.J.M.1963,Bulletin of the Astronomical Institutes of the Netherlands,vol.17,pp.22- 28.[2] server Vanicek,P.1971,Astrophysics and Space Science,vol.12,pp.10-33.[3] (Nort Scargle,J.D.1982,Astrophysical Journal,vol.263,pp.835-853.[4] computer, Ameri to make one paper UnN电.t THE Horne,J.H.,and Baliunas,S.L.1986,Astrophysical Journal,vol.302,pp.757-763.[5] Press,W.H.and Rybicki,G.B.1989,Astrophysica/Journal,vol.338,pp.277-280.[6] ART ictly proh Programs 13.9 Computing Fourier Integrals Using the FFT 6 Not uncommonly,one wants to calculate accurate numerical values for integrals of the form OF SCIENTIFIC COMPUTING (ISBN 198918920 I= eh(t)dt, (13.9.1) or the equivalent real and imaginary parts Numerica 10621 ^6 Ic= cos(wt)h(t)dt I= sin(wt)h(t)dt, (13.9.2) 43106 and one wants to evaluate this integral for many different values of w.In cases of interest,h(t) is often a smooth function,but it is not necessarily periodic in a,b,nor does it necessarily (outside Recipes go to zero at a or b.While it seems intuitively obvious that the force majeure of the FFT ought to be applicable to this problem,doing so turns out to be a surprisingly subtle matter, North as we will now see. Let us first approach the problem naively,to see where the difficulty lies.Divide the interval a,b into M subintervals,where M is a large integer,and define A=-a ，与三a+jA,h三h6),j=0,M (13.9.3) Notice that ho h(a)and hh(b),and that there are M+1 values hj.We can approximate the integral I by a sum, M-1 I≈△) hiexp(iwtj) (13.9.4) =0584 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 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). 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)); } } } 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´ıcek, 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 FFT Not uncommonly, one wants to calculate accurate numerical values for integrals of the form I =
b a e iωth(t)dt , (13.9.1) or the equivalent real and imaginary parts Ic =
b a cos(ωt)h(t)dt Is =
b a sin(ωt)h(t)dt , (13.9.2) and 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 necessarily go to zero at a or b. While it seems intuitively obvious that the force majeure of the FFT ought 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 the interval [a, b] into M subintervals, where M is a large integer, and define ∆ ≡ b − a M , tj ≡ a + j∆ , hj ≡ h(tj ), j = 0,...,M (13.9.3) Notice that h0 = h(a) and hM = h(b), and that there are M + 1 values hj . We can approximate the integral I by a sum, I ≈ ∆ M −1 j=0 hj exp(iωtj ) (13.9.4)