ix=(int)x; if (x == (float)ix) yy[ix] += y; else { ilo=LMIN(LMAX((long)(*m+),1),n-m+1); ihi=ilo+m-1; nden=nfac[m]; fac=x-ilo; for (j=ilo+1;j=ilo;j--) | 584 Chapter 13. Fourier and Spectral Applications ix int x if x float ix yy ix y else ilo LMIN LMAX long m 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 . 1976 Astrophysics and Space Science vol. 39 pp. 447-462. 1 Barning . 1963 Bulletin of the Astronomical Institutes of the Netherlands vol. 17 pp. 22- 28. 2 Vanicek P. 1971 Astrophysics and Space Science vol. 12 pp. 10-33. 3 Scargle . 1982 Astrophysical Journal vol. 263 pp. 835-853. 4 Horne . and Baliunas . 1986 Astrophysical Journal vol. 302 pp. 757-763. 5 Press . and Rybicki . 1989 Astrophysical Journal vol. 338 pp. 277-280. 6 Computing FourierIntegrals Using the FFT Not uncommonly one wants to calculate accurate numerical values for integrals of the form I f a e th t dt or the equivalent real and imaginary parts Ic cos t h t dt a Is sin t h t dt a 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 A tj a jA hj h tj j 0 . . M Notice that ho h a and hM h b and that there are M 1 values hj. We can approximate the integral I by a sum Sample page from NUMERICAL RECIPES IN C THE ART OF SCIENTIFIC COMPUTING ISBN 0-521-43108-5 M-1 I A hj exp i tj j 0 Computing Fourier Integrals Using the FFT 585 which is at any rate first-order accurate. If we centered the hj s and the tj s in the intervals we could be .
