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) =0
584 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). 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)
13.9 Computing Fourier Integrals Using the FFT 585 which is at any rate first-order accurate.(If we centered the h,'s and the t,'s in the intervals, we could be accurate to second order.)Now for certain values of w and M,the sum in equation (13.9.4)can be made into a discrete Fourier transform,or DFT,and evaluated by the fast Fourier transform(FFT)algorithm.In particular,we can choose M to be an integer power of 2,and define a set of special w's by n (13.9.5) where m has the values m =0,1,...,M/2-1.Then equation (13.9.4)becomes M-1 I(wm)≈△euma h,e2rim时/M=△DFT(ho.hM-i月m (13.9.6) 三 人0 8 Equation (13.9.6),while simple and clear,is emphatically not recommended for use:It is likely 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 w 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 that it is easily swamped by first-order,or even (with centered values)second-order,truncation error.Furthermore,the characteristic "small parameter"that occurs in the error term is not A/(b-a)=1/M,as it would be if the integrand were not oscillatory,but wA,which can be as large as for w's within the Nyquist interval of the DFT (cf.equation 13.9.5).The result is that equation (13.9.6)becomes systematically inaccurate as w increases. It is a sobering exercise to implement equation(13.9.6)for an integral that can be done 线 THE analytically,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 h,we ART can approximate the function h(t)everywhere in the interval a,b by interpolation on nearby hj's.The simplest case is linear interpolation,using the two nearest h's,one to the left and Programs one to the right.A higher-order interpolation,e.g.,would be cubic interpolation,using two points to the left and two to the right-except in the first and last subintervals,where we must interpolate with three hi's on one side,one on the other. The formulas for such interpolation schemes are (piecewise)polynomial in the inde- pendent variable t,but with coefficients that are of course linear in the function values hj.Although one does not usually think of it in this way,interpolation can be viewed as approximating a function by a sum of kernel functions(which depend only on the interpolation OF SCIENTIFIC COMPUTING (ISBN scheme)times sample values (which depend only on the function).Let us write 1988-19920 0（）+∑（之） (13.9.7) =0 j=endpoints Here (s)is the kernel function of an interior point:It is zero for s sufficiently negative or sufficiently positive,and becomes nonzero only when s is in the range where the Numerical Recipes 10-621 -43106 h;multiplying it is actually used in the interpolation.We always have (0)=1 and (m)=0,m=±l,±2，，since interpolation right on a sample point should give the sampled function value.For linear interpolation (s)is piecewise linear,rises from 0 to 1 (outside for s in (-1,0),and falls back to 0 for s in (0,1).For higher-order interpolation,(s)is Software. made up piecewise of segments of Lagrange interpolation polynomials.It has discontinuous 首 derivatives at integer values of s,where the pieces join,because the set of points used in the 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 the special endpoint kernels (s).Actually,for reasons that will become clearer below,we have included all the points in the first sum (with kernel )so the ,'s are actually differences between true endpoint kernels and the interior kernel It is a tedious,but straightforward, exercise to write down all the(s)'s for any particular order of interpolation,each one consisting of differences of Lagrange interpolating polynomials spliced together piecewise. Now apply the integral operator dtexp(iwt)to both sides of equation (13.9.7). interchange the sums and integral,and make the changes of variable s=(t-)/A in the
13.9 Computing Fourier Integrals Using the FFT 585 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). which 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 in equation (13.9.4) can be made into a discrete Fourier transform, or DFT, and evaluated by the fast Fourier transform (FFT) algorithm. In particular, we can choose M to be an integer power of 2, and define a set of special ω’s by ωm∆ ≡ 2πm M (13.9.5) where m has the values m = 0, 1,...,M/2 − 1. Then equation (13.9.4) becomes I(ωm) ≈ ∆eiωma M −1 j=0 hj e2πimj/M = ∆eiωma[DFT(h0 ...hM−1)]m (13.9.6) Equation (13.9.6), while simple and clear, is emphatically not recommended for use: It is likely 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 that it is easily swamped by first-order, or even (with centered values) second-order, truncation error. 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 be as large as π for ω’s within the Nyquist interval of the DFT (cf. equation 13.9.5). The result is 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 done analytically, 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 , we can approximate the function h(t) everywhere in the interval [a, b] by interpolation on nearby hj ’s. The simplest case is linear interpolation, using the two nearest hj ’s, one to the left and one to the right. A higher-order interpolation, e.g., would be cubic interpolation, using two points to the left and two to the right — except in the first and last subintervals, where we must 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 values hj . Although one does not usually think of it in this way, interpolation can be viewed as approximating a function by a sum of kernel functions (which depend only on the interpolation scheme) times sample values (which depend only on the function). Let us write h(t) ≈ M j=0 hj ψ t − tj ∆ + j=endpoints hj ϕj t − tj ∆ (13.9.7) Here ψ(s) is the kernel function of an interior point: It is zero for s sufficiently negative or sufficiently positive, and becomes nonzero only when s is in the range where the hj 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 the sampled function value. For linear interpolation ψ(s) is piecewise linear, rises from 0 to 1 for s in (−1, 0), and falls back to 0 for s in (0, 1). For higher-order interpolation, ψ(s) is made up piecewise of segments of Lagrange interpolation polynomials. It has discontinuous derivatives at integer values of s, where the pieces join, because the set of points used in the 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 the special endpoint kernels ϕj (s). Actually, for reasons that will become clearer below, we have included all the points in the first sum (with kernel ψ), so the ϕj ’s are actually differences between 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 one consisting of differences of Lagrange interpolating polynomials spliced together piecewise. Now apply the integral operator b 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 the
586 Chapter 13.Fourier and Spectral Applications first sum,s=(t-a)/A in the second sum.The result is I≈△euwa wO)∑h,eo+∑ ha(o) (13.9.8) j=0 j=endpoints Here0≡w△，and the functions W(O)anda(o)are defined by W(O)三 dses(s) (13.9.9) a(0)≡ dseis(s-j) (13.9.10) -00 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 algorithm for applying"endpoint corrections"to a sum which (as we will see)can be done using the FFT,giving a result with high-order accuracy. We will consider only interpolations that are left-right symmetric. Then symmetry implies PM-j(s)=pj(-s) aM-(0)=eaMa(0)=eu-a(0) (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 2 To do so,choose some N that is an integer power of 2 with N>M+1.(Note that M need not be a power of two,so M N-1 is allowed.)If N M+1,define h;=0,M+1<j<N-1,i.e.,"zero pad"the array of hi's so that j takes on the range THE 0<j<N-1.Then the sum can be done as a DFT for the special valuesw=wn given by 公 ART A=2=0n=0,1岁-1 (13.9.12) 9 Program For fixed M,the larger N is chosen,the finer the sampling in frequency space.The value M,on the other hand,determines the highest frequency sampled,since A decreases with increasing M (equation 13.9.3),and the largest value of wA 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 as I(wn)=△ena{W(O)DFT(ho..hN-1月n OF SCIENTIFIC COMPUTING (ISBN +ao(0)ho+a1(0)h1+a2(0)h2+a3(0)h3+.. 198918920 teiv(b-a)[ o()h+i-+-+ai(0ha-+.. 10-621 (13.9.13) For cubic (or lower)polynomial interpolation,at most the terms explicitly shown above Numerical Recipes 431085 are nonzero;the ellipses (...)can therefore be ignored,and we need explicit forms only for the functions W,ao,1,2,a3,calculated with equations (13.9.9)and (13.9.10).We have worked these out for you,in the trapezoidal (second-order)and cubic (fourth-order)cases. (outside Here are the results,along with the first few terms of their power series expansions for small 0: North Software. Trapezoidal order: ying of w0=200≈1-立2+a0- 1 1 02 20160 a0(0)=-L-cos0》+:9-sin0】 02 02 720 592+1 (6-120 5040 362880 a1=a2=ag=0
586 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). first sum, s = (t − a)/∆ in the second sum. The result is I ≈ ∆eiωa W(θ) M j=0 hj eijθ + j=endpoints hjαj (θ) (13.9.8) Here θ ≡ ω∆, and the functions W(θ) and αj (θ) are defined by W(θ) ≡ ∞ −∞ ds eiθsψ(s) (13.9.9) αj (θ) ≡ ∞ −∞ ds eiθsϕj (s − j) (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 algorithm for applying “endpoint corrections” to a sum which (as we will see) can be done using the FFT, giving a result with high-order accuracy. We will consider only interpolations that are left-right symmetric. Then symmetry implies ϕ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 that M need not be a power of two, so M = N − 1 is allowed.) If N>M + 1, define hj ≡ 0, M + 1 < j ≤ N − 1, i.e., “zero pad” the array of hj ’s so that j takes on the range 0 ≤ j ≤ N − 1. Then the sum can be done as a DFT for the special values ω = ωn given by ωn∆ ≡ 2πn N ≡ θ n = 0, 1,..., N 2 − 1 (13.9.12) For fixed M, the larger N is chosen, the finer the sampling in frequency space. The value M, on the other hand, determines the highest frequency sampled, since ∆ decreases with 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 as I(ωn)=∆eiωna W(θ)[DFT(h0 ...hN−1)]n + α0(θ)h0 + α1(θ)h1 + α2(θ)h2 + α3(θ)h3 + ... + 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 above are nonzero; the ellipses (...) can therefore be ignored, and we need explicit forms only for the functions W, α0, α1, α2, α3, calculated with equations (13.9.9) and (13.9.10). We have worked 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(θ) = 2(1 − cos θ) θ2 ≈ 1 − 1 12 θ2 + 1 360 θ4 − 1 20160 θ6 α0(θ) = − (1 − cos θ) θ2 + i (θ − sin θ) θ2 ≈ − 1 2 + 1 24 θ2 − 1 720 θ4 + 1 40320 θ6 + iθ 1 6 − 1 120 θ2 + 1 5040 θ4 − 1 362880 θ6 α1 = α2 = α3 = 0
13.9 Computing Fourier Integrals Using the FFT 587 Cubic order: w间)=()B-4ms0+cs291-品0+品的 /6+92 23 a00=-42+502)+6+62)8cos9-cos20+,-120+603)+6+92)sin29 604 604 86 3+ 22a6+0（+05-飞 -2835 04+ 467775 a1(0=143-g2)-76+0)cos0+:300-56+92)sim0 604 609 25920006+0(7、1 2+ 11 13 83 3456 （72-1681 72576 5987520 鱼 20=-48-02)+26+9)0s0+:-120+26+g2)sim0 1.00 304 304 7 1 284+立一8°+0-90十206 83 11 13 64800 90720 +7484400 a3(0)=28-69)-6+02)os0+00-6+02)sm0 604 604 1 7 .11g2+01020-2592009+8 1 24-180 360840 82+ 11 13 6 America Press. 362880 29937600 ART The program dftcor,below,implements the endpoint corrections for the cubic case. Given input valuesofw,A,a,b,and an array with the eight values ho,...,hs,hM-3,...,hM, it returns the real and imaginary parts of the endpoint corrections in equation (13.9.13),and the factor W().The code is turgid,but only because the formulas above are complicated.The formulas have cancellations to high powers of 0.It is therefore necessary to compute the right- hand sides in double precision,even when the corrections are desired only to single precision. to dir It is also necessary to use the series expansion for small values of 6.The optimal cross-over value of 6 depends on your machine's wordlength,but you can always find it experimentally SCIENTIFIC COMPUTING (ISBN as the largest value where the two methods give identical results to machine precision. 19841820 #include <math.h> void dftcor(float w,float delta,float a,float b,float endpts[], float *corre,float *corim,float *corfac) For an integral approximated by a discrete Fourier transform,this routine computes the cor- idge.org uurggoglrion Numerical Recipes 10.621 43106 rection factor that multiplies the DFT and the endpoint correction to be added.Input is the angular frequency w,stepsize delta,lower and upper limits of the integral a and b,while the array endpts contains the first 4 and last 4 function values.The correction factor W()is (outside returned as corfac,while the real and imaginary parts of the endpoint correction are returned North Software. as corre and corim. void nrerror(char error_text []) float a0i,aor,ali,alr,a2i,a2r,a3i,a3r,arg,c,cl,cr,s,sl,sr,t; float t2,t4,t6; double cth,ctth,spth2,sth,sth4i,stth,th,th2,th4,tmth2,tth4i; th=y*delta: if (a >b lI th 0.0e0 Il th 3.1416e0)nrerror("bad arguments to dftcor"); if (fabs(th)<5.0e-2){ Use series. t=th; t2=t*t; t4=t2*t2; t6=t4*t2;
13.9 Computing Fourier Integrals Using the FFT 587 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). Cubic order: W(θ) = 6 + θ2 3θ4 (3 − 4 cos θ + cos 2θ) ≈ 1 − 11 720 θ4 + 23 15120 θ6 α0(θ) = (−42 + 5θ2) + (6 + θ2)(8 cos θ − cos 2θ) 6θ4 + i (−12θ + 6θ3) + (6 + θ2) sin 2θ 6θ4 ≈ − 2 3 + 1 45 θ2 + 103 15120 θ4 − 169 226800 θ6 + iθ 2 45 + 2 105 θ2 − 8 2835 θ4 + 86 467775 θ6 α1(θ) = 14(3 − θ2) − 7(6 + θ2) cos θ 6θ4 + i 30θ − 5(6 + θ2) sin θ 6θ4 ≈ 7 24 − 7 180 θ2 + 5 3456 θ4 − 7 259200 θ6 + iθ 7 72 − 1 168 θ2 + 11 72576 θ4 − 13 5987520 θ6 α2(θ) = −4(3 − θ2) + 2(6 + θ2) cos θ 3θ4 + i −12θ + 2(6 + θ2) sin θ 3θ4 ≈ − 1 6 + 1 45 θ2 − 5 6048 θ4 + 1 64800 θ6 + iθ − 7 90 + 1 210 θ2 − 11 90720 θ4 + 13 7484400 θ6 α3(θ) = 2(3 − θ2) − (6 + θ2) cos θ 6θ4 + i 6θ − (6 + θ2) sin θ 6θ4 ≈ 1 24 − 1 180 θ2 + 5 24192 θ4 − 1 259200 θ6 + iθ 7 360 − 1 840 θ2 + 11 362880 θ4 − 13 29937600 θ6 The program dftcor, below, implements the endpoint corrections for the cubic case. Given input values of ω, ∆, a, b, and an array with the eight values h0 ,...,h3, hM−3,...,hM , it returns the real and imaginary parts of the endpoint corrections in equation (13.9.13), and the factor W(θ). The code is turgid, but only because the formulas above are complicated. The formulas have cancellations to high powers of θ. It is therefore necessary to compute the righthand sides in double precision, even when the corrections are desired only to single precision. It is also necessary to use the series expansion for small values of θ. The optimal cross-over value of θ depends on your machine’s wordlength, but you can always find it experimentally as the largest value where the two methods give identical results to machine precision. #include <math.h> void dftcor(float w, float delta, float a, float b, float endpts[], float *corre, float *corim, float *corfac) For an integral approximated by a discrete Fourier transform, this routine computes the correction factor that multiplies the DFT and the endpoint correction to be added. Input is the angular frequency w, stepsize delta, lower and upper limits of the integral a and b, while the array endpts contains the first 4 and last 4 function values. The correction factor W(θ) is returned as corfac, while the real and imaginary parts of the endpoint correction are returned as corre and corim. { void nrerror(char error_text[]); float a0i,a0r,a1i,a1r,a2i,a2r,a3i,a3r,arg,c,cl,cr,s,sl,sr,t; float t2,t4,t6; double cth,ctth,spth2,sth,sth4i,stth,th,th2,th4,tmth2,tth4i; th=w*delta; if (a >= b || th < 0.0e0 || th > 3.1416e0) nrerror("bad arguments to dftcor"); if (fabs(th) < 5.0e-2) { Use series. t=th; t2=t*t; t4=t2*t2; t6=t4*t2;
588 Chapter 13. Fourier and Spectral Applications *c0rfac=1.0-(11.0/720.0)*t4+(23.0/15120.0)*t6: a0r=(-2.0/3.0)+t2/45.0+(103.0/15120.0)*t4-(169.0/226800.0)*t6; a1r=(7.0/24.0)-(7.0/180.0)*t2+(5.0/3456.0)*t4-(7.0/259200.0)*t6; a2r=(-1.0/6.0)+t2/45.0-(5.0/6048.0)*t4+t6/64800.0; a3r=(1.0/24.0)-t2/180.0+(5.0/24192.0)*t4-t6/259200.0: a01=t*(2.0/45.0+(2.0/105.0)*t2-(8.0/2835.0)*t4+(86.0/467775.0)*t6); a11=t*(7.0/72.0-t2/168.0+(11.0/72576.0)*t4-(13.0/5987520.0)*t6); a21=t*(-7.0/90.0+t2/210.0-(11.0/90720.0)*t4+(13.0/7484400.0)*t6) 31=t*(7.0/360.0-t2/840.0+(11.0/362880.0)*t4-(13.0/29937600.0)*t6); else Use trigonometric formulas in double precision. cth=cos(th); sth=sin(th); http://www.nr read able files Permission is ctth=cth*cth-sth*sth; stth=2.0e0*sth*cth; th2-th*th; th4=th2*th2: tmth2=3.0e0-th2; spth2=6.0e0+th2: (including this one) 19881992 sth41=1.0/(6.0e0*th4) tth4i=2.0e0*sth4i; 11-800-872 *corfac=tth4i*spth2*(3.0e0-4.0e0*cth+ctth); aOr=sth4i*(-42.0e0+5.0e0*th2+spth2*(8.0e0*cth-ctth)); a0i=sth4i*(th*(-12.0e0+6.0e0*th2)+spth2*stth); alr=sth4i*(14.0e0*tmth2-7.0e0*spth2*cth); ali=sth4i*(30.0e0*th-5.0e0*spth2*sth); a2r=tth4i*(-4.0e0*tmth2+2.0e0*spth2*cth): a21=tth4i*(（-12.0e0*th+2.0e0*spth2*sth); (North America to any server computer,is users to make one paper /Cambridge University Press.Programs from NUMERICAL RECIPES IN THE a3r=sth4i*(2.0e0*tmth2-spth2*cth); ART a3i=sth4i*(6.0e0*th-spth2*sth); cl=aOr*endpts[1]+alr*endpts [2]+a2r*endpts [3]+a3r*endpts [4]; sl=aoi*endpts[1]+ali*endpts [2]+a2i*endpts [3]+a3i*endpts [4]; copy for their cr=aOr*endpts [8]+alr*endpts [7]+a2r*endpts [6]+a3r*endpts [5]; sr =-aoi*endpts[8]-ali*endpts[7]-a2i*endpts[6]-a3i*endpts[5]; st st arg=w*(b-a); c=cos(arg); s=sin(arg); *corre=cl+c*cr-s*sr; *corim=sl+s*cr+c*sr; rectcustser OF SCIENTIFIC COMPUTING(ISBN 0-521 Since the use of dftcor can be confusing,we also give an illustrative program dftint which uses dftcor to compute equation (13.9.1)for general a,b,w,and h(t).Several points within this program bear mentioning:The parameters M and NDFT correspond to M and N in the above discussion.On successive calls,we recompute the Fourier transform only if 1988-1992 by Numerical Recipes -431085 a or b or h(t)has changed. Since dftint is designed to work for any value of w satisfying wA<not just the (outside special values returned by the DFT (equation 13.9.12),we do polynomial interpolation of degree MPOL on the DFT spectrum.You should be warned that a large factor of oversampling Software. (NM)is required for this interpolation to be accurate.After interpolation,we add the North endpoint corrections from dftcor,which can be evaluated for any w. Ame While dftcor is good at what it does,dftint is illustrative only.It is not a general purpose program,because it does not adapt its parameters M,NDFT,MPOL,or its interpolation scheme,to any particular function h(t).You will have to experiment with your own application #include <math.h> #include "nrutil.h" #define M 64 #define NDFT 1024 #define MPOL 6 #def1neTW0PI(2.0*3.14159265)
588 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). *corfac=1.0-(11.0/720.0)*t4+(23.0/15120.0)*t6; a0r=(-2.0/3.0)+t2/45.0+(103.0/15120.0)*t4-(169.0/226800.0)*t6; a1r=(7.0/24.0)-(7.0/180.0)*t2+(5.0/3456.0)*t4-(7.0/259200.0)*t6; a2r=(-1.0/6.0)+t2/45.0-(5.0/6048.0)*t4+t6/64800.0; a3r=(1.0/24.0)-t2/180.0+(5.0/24192.0)*t4-t6/259200.0; a0i=t*(2.0/45.0+(2.0/105.0)*t2-(8.0/2835.0)*t4+(86.0/467775.0)*t6); a1i=t*(7.0/72.0-t2/168.0+(11.0/72576.0)*t4-(13.0/5987520.0)*t6); a2i=t*(-7.0/90.0+t2/210.0-(11.0/90720.0)*t4+(13.0/7484400.0)*t6); a3i=t*(7.0/360.0-t2/840.0+(11.0/362880.0)*t4-(13.0/29937600.0)*t6); } else { Use trigonometric formulas in double precision. cth=cos(th); sth=sin(th); ctth=cth*cth-sth*sth; stth=2.0e0*sth*cth; th2=th*th; th4=th2*th2; tmth2=3.0e0-th2; spth2=6.0e0+th2; sth4i=1.0/(6.0e0*th4); tth4i=2.0e0*sth4i; *corfac=tth4i*spth2*(3.0e0-4.0e0*cth+ctth); a0r=sth4i*(-42.0e0+5.0e0*th2+spth2*(8.0e0*cth-ctth)); a0i=sth4i*(th*(-12.0e0+6.0e0*th2)+spth2*stth); a1r=sth4i*(14.0e0*tmth2-7.0e0*spth2*cth); a1i=sth4i*(30.0e0*th-5.0e0*spth2*sth); a2r=tth4i*(-4.0e0*tmth2+2.0e0*spth2*cth); a2i=tth4i*(-12.0e0*th+2.0e0*spth2*sth); a3r=sth4i*(2.0e0*tmth2-spth2*cth); a3i=sth4i*(6.0e0*th-spth2*sth); } cl=a0r*endpts[1]+a1r*endpts[2]+a2r*endpts[3]+a3r*endpts[4]; sl=a0i*endpts[1]+a1i*endpts[2]+a2i*endpts[3]+a3i*endpts[4]; cr=a0r*endpts[8]+a1r*endpts[7]+a2r*endpts[6]+a3r*endpts[5]; sr = -a0i*endpts[8]-a1i*endpts[7]-a2i*endpts[6]-a3i*endpts[5]; arg=w*(b-a); c=cos(arg); s=sin(arg); *corre=cl+c*cr-s*sr; *corim=sl+s*cr+c*sr; } Since the use of dftcor can be confusing, we also give an illustrative program dftint which uses dftcor to compute equation (13.9.1) for general a, b, ω, and h(t). Several points within this program bear mentioning: The parameters M and NDFT correspond to M and N in the above discussion. On successive calls, we recompute the Fourier transform only if a or b or h(t) has changed. Since dftint is designed to work for any value of ω satisfying ω∆ < π, not just the special values returned by the DFT (equation 13.9.12), we do polynomial interpolation of degree MPOL on the DFT spectrum. You should be warned that a large factor of oversampling (N M) is required for this interpolation to be accurate. After interpolation, we add the endpoint corrections from dftcor, which can be evaluated for any ω. While dftcor is good at what it does, dftint is illustrative only. It is not a general purpose program, because it does not adapt its parameters M, NDFT, MPOL, or its interpolation scheme, to any particular function h(t). You will have to experiment with your own application. #include <math.h> #include "nrutil.h" #define M 64 #define NDFT 1024 #define MPOL 6 #define TWOPI (2.0*3.14159265)