f13-9

13.9ComputingFourierIntegralsUsingtheFFT 577Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).yy(j)=yy(j)+y*fac/(nden*(x-j)) enddo 12 endifreturnEND 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 ComputingFourierIntegralsUsingtheFFT Not uncommonly, one wants to calculate accurate numerical values for integrals of the form I=/integraldisplayb aeiωth(t)dt , (13.9.1 ) or the equivalent real and imaginary parts Ic=/integraldisplayb acos(ωt)h(t)dt I s=/integraldisplayb asin(ωt)h(t)dt , (13.9.2 ) andonewantstoevaluatethisintegralformanydifferentvaluesof ω. Incasesofinterest, h(t) is often a smooth function, but it is not necessarily periodic in [a, b], nor does it necessarily go to zero at aorb. 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]intoMsubintervals, where Mis a large integer, and define ∆≡b−a M,t j≡a+j∆,h j≡h(tj),j =0,...,M (13.9.3 ) Notice that h0=h(a)andhM=h(b), and that there are M+1values hj. We can approximate the integral Iby a sum, I≈∆M−1/summationdisplay j=0hjexp(iωtj)( 13.9.4 ) 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 ωandM, the sum in equation (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 Mto be an integer power of 2, and define a set of special ω’s by ω m∆≡2πm M(13.9.5 ) 578 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).where mhas the values m=0,1,...,M/ 2−1. Then equation (13.9.4) becomes I(ωm)≈∆eiω maM−1/summationdisplay j=0hje2πimj/M=∆eiω ma[DFT(h0...h M−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! Theproblemliesintheoscillatorynatureoftheintegral(13.9.1). If h(t)isatallsmooth, andif ωislargeenough toimplyseveralcyclesintheinterval [a, b]—infact, ωminequation (13.9.5) gives exactly mcycles — then the value of Iis typically very small, so small that it 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,asitwouldbeiftheintegrand werenotoscillatory,but ω∆,whichcanbe 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. Letusthereforeturntoamoresophisticatedtreatment. Giventhesampledpoints h j,we can approximate the function h(t)everywhere intheinterval [a, b]by interpolation onnearby 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 twopoints to the left and two to the right — except in the first and last subintervals, where wemust interpolate with three h j’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 approximatingafunctionbyasumofkernelfunctions(whichdependonlyontheinterpolationscheme) times sample values (which depend only on the function). Let us write h(t)≈M/summationdisplay j=0hjψ/parenleftbiggt−tj ∆/parenrightbigg +/summationdisplay j=endpointshjϕj/parenleftbiggt−tj ∆/parenrightbigg (13.9.7 ) Here ψ(s)is the kernel function of an interior point: It is zero for ssufficiently negative or sufficiently positive, and becomes nonzero only when sis in the range where the hjmultiplying 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 forsin(−1,0), and falls back to 0 for sin(0,1). For higher-order interpolation, ψ(s)is made 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 in the interpolation changes discretely. As already remarked, the subintervals closest to aandbrequire different (noncentered) interpolation formulas. This is reflected in equation (13.9.7) by the second sum, with thespecialendpoint kernels ϕ j(s). Actually, forreasons thatwillbecome clearerbelow, wehave included allthe points in the firstsum (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/integraltextb adtexp(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 first sum, s=(t−a)/∆in the second sum. The result is I≈∆eiωa/bracketleftBigg W(θ)M/summationdisplay j=0hjeijθ+/summationdisplay j=endpointshjαj(θ)/bracketrightBigg (13.9.8 ) Here θ≡ω∆, and the functions W(θ)andαj(θ)are defined by W(θ)≡/integraldisplay∞ −∞ds eiθsψ(s)( 13.9.9 ) 13.9ComputingFourierIntegralsUsingtheFFT 579Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).αj(θ)≡/integraldisplay∞ −∞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 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 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 Nthat is an integer power of 2 with N≥M+1. (Note that Mneed not be a power of two, so M=N−1is allowed.) If N>M +1, define hj≡0,M +1<j≤N−1, i.e., “zero pad” the array of hj’s so that jtakes on the range 0≤j≤N−1. Then the sum can be done as a DFTfor the special values ω=ωngiven by ωn∆≡2πn N≡θn =0,1,...,N 2−1( 13.9.12 ) For fixed M, the larger Nis chosen, the finer the sampling in frequency space. The value M, on the other hand, determines the highestfrequency 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 leasta 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/braceleftbigg W(θ)[DFT(h0...h N−1)]n +α0(θ)h0+α1(θ)h1+α2(θ)h2+α3(θ)h3+... +eiω(b−a)/bracketleftBig α* 0(θ)hM+α* 1(θ)hM−1+α* 2(θ)hM−2+α* 3(θ)hM−3+.../bracketrightBig/bracerightbigg (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.Herearetheresults,alongwiththefirstfewtermsoftheirpowerseriesexpansionsforsmall θ: 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θ/parenleftbigg1 6−1 120θ2+1 5040θ4−1 362880θ6/parenrightbigg α1=α2=α3=0 580 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Cubic order: W(θ)=/parenleftbigg6+θ2 3θ4/parenrightbigg (3−4c o s θ+c o s2 θ)≈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)s i n2 θ 6θ4 ≈−2 3+1 45θ2+103 15120θ4−169 226800θ6+iθ/parenleftbigg2 45+2 105θ2−8 2835θ4+86 467775θ6/parenrightbigg α1(θ)=14(3−θ2)−7(6 + θ2)c o sθ 6θ4+i30θ−5(6 + θ2)s i nθ 6θ4 ≈7 24−7 180θ2+5 3456θ4−7 259200θ6+iθ/parenleftbigg7 72−1 168θ2+11 72576θ4−13 5987520θ6/parenrightbigg α2(θ)=−4(3−θ2)+2 ( 6+ θ2)c o sθ 3θ4+i−12θ+2 ( 6+ θ2)s i nθ 3θ4 ≈−1 6+1 45θ2−5 6048θ4+1 64800θ6+iθ/parenleftbigg −7 90+1 210θ2−11 90720θ4+13 7484400θ6/parenrightbigg α3(θ)=2(3−θ2)−(6 + θ2)c o sθ 6θ4+i6θ−(6 + θ2)s i nθ 6θ4 ≈1 24−1 180θ2+5 24192θ4−1 259200θ6+iθ/parenleftbigg7 360−1 840θ2+11 362880θ4−13 29937600θ6/parenrightbigg The program dftcor, below, implements the endpoint corrections for the cubic case. Giveninputvaluesof ω,∆,a ,b ,andanarraywiththeeightvalues h0,...,h 3,hM−3,...,h M, itreturnstherealandimaginarypartsoftheendpointcorrectionsinequation(13.9.13),andthefactor W(θ). The code is turgid, but only because the formulas above are complicated. The formulashavecancellationstohighpowersof θ. Itisthereforenecessarytocomputetheright- hand sidesindouble precision, evenwhen thecorrections aredesiredonly tosingleprecision.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. SUBROUTINE dftcor(w,delta,a,b,endpts,corre,corim,corfac) REAL a,b,corfac,corim,corre,delta,w,endpts(8) For an integral approximated by a discrete Fourier transform, this routine computes thecorrection factor that multiplies the DFT and the endpoint correction to be added. Inputis the angular frequency w,s t e p s i z e delta , lower and upper limits of the integral aandb, 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 andcorim . REAL a0i,a0r,a1i,a1r,a2i,a2r,a3i,a3r,arg,c,cl,cr,s,sl,sr,t, * t2,t4,t6 DOUBLE PRECISION cth,ctth,spth2,sth,sth4i,stth,th,th2,th4, * tmth2,tth4i th=w*delta if (a.ge.b.or.th.lt.0.d0.or.th.gt.3.1416d0) * pause ’bad arguments to dftcor’ if(abs(th).lt.5.d-2)then Use series. t=tht2=t*tt4=t2*t2 t6=t4*t2 corfac=1.-(11./720.)*t4+(23./15120.)*t6a0r=(-2./3.)+t2/45.+(103./15120.)*t4-(169./226800.)*t6 a1r=(7./24.)-(7./180.)*t2+(5./3456.)*t4-(7./259200.)*t6 13.9ComputingFourierIntegralsUsingtheFFT 581Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).a2r=(-1./6.)+t2/45.-(5./6048.)*t4+t6/64800. a3r=(1./24.)-t2/180.+(5./24192.)*t4-t6/259200. a0i=t*(2./45.+(2./105.)*t2-(8./2835.)*t4+(86./467775.)*t6) a1i=t*(7./72.-t2/168.+(11./72576.)*t4-(13./5987520.)*t6)a2i=t*(-7./90.+t2/210.-(11./90720.)*t4+(13./7484400.)*t6) a3i=t*(7./360.-t2/840.+(11./362880.)*t4-(13./29937600.)*t6) else Use trigonometric formulas in double precision. cth=cos(th)sth=sin(th) ctth=cth**2-sth**2 stth=2.d0*sth*cthth2=th*thth4=th2*th2 tmth2=3.d0-th2 spth2=6.d0+th2sth4i=1./(6.d0*th4) tth4i=2.d0*sth4i corfac=tth4i*spth2*(3.d0-4.d0*cth+ctth)a0r=sth4i*(-42.d0+5.d0*th2+spth2*(8.d0*cth-ctth))a0i=sth4i*(th*(-12.d0+6.d0*th2)+spth2*stth) a1r=sth4i*(14.d0*tmth2-7.d0*spth2*cth) a1i=sth4i*(30.d0*th-5.d0*spth2*sth)a2r=tth4i*(-4.d0*tmth2+2.d0*spth2*cth) a2i=tth4i*(-12.d0*th+2.d0*spth2*sth) a3r=sth4i*(2.d0*tmth2-spth2*cth)a3i=sth4i*(6.d0*th-spth2*sth) endif 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*srreturnEND Since the use of dftcorcan be confusing, we also give an illustrative program dftint which uses dftcorto compute equation (13.9.1) for general a, b, ω,andh(t). Several points within this program bear mentioning: The parameters MandNDFTcorrespond to MandN in the above discussion. On successive calls, we recompute the Fourier transform only if a orbhas changed. (We should also recompute if h(t)has changed, but FORTRAN doesn’t provide a way for us to test this.) Since dftintis 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 ofdegree MPOLon the DFTspectrum. You should be warned that a large factor of oversampling (N/greatermuchM) is required for this interpolation to be accurate. After interpolation, we add the endpoint corrections from dftcor, which can be evaluated for any ω. While dftcoris good at what it does, dftintis illustrative only. It is not a general purpose program, because it does not adapt its parameters M,NDFT,MPOL, or its interpolation scheme,toanyparticularfunction h(t). Youwillhavetoexperimentwithyourownapplication. SUBROUTINE dftint(func,a,b,w,cosint,sinint) INTEGER M,NDFT,MPOLREAL a,b,cosint,sinint,w,func,TWOPIPARAMETER (M=64,NDFT=1024,MPOL=6,TWOPI=2.*3.14159265) EXTERNAL func C USES dftcor,func,polint,realft 582 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Example program illustrating how to use the routine dftcor . The user supplies an external function func that returns the quantity h(t). The routine then returns/integraltextb acos(ωt)h(t)dt ascosint and/integraltextb asin(ωt)h(t)dtassinint . Parameters: The values of M,NDFT ,a n d MPOL are merely illustrative and should be opti- mized for your particular application. Mis the number of subintervals, NDFT is the length of the FFT (a power of 2), and MPOL is the degree of polynomial interpolation used to obtain the desired frequency from the FFT. INTEGER init,j,nn REAL aold,bold,c,cdft,cerr,corfac,corim,corre,delta,en,s, * sdft,serr,cpol(MPOL),data(NDFT),endpts(8),spol(MPOL),* xpol(MPOL) SAVE init,aold,bold,delta,data,endpts DATA init/0/,aold/-1.e30/,bold/-1.e30/if (init.ne.1.or.a.ne.aold.or.b.ne.bold) then Do we need to initialize, or is only ω changed? init=1 aold=a bold=bdelta=(b-a)/M do 11j=1,M+1 Load the function values into the data array. data(j)=func(a+(j-1)*delta) enddo 11 do12j=M+2,NDFT Zero pad the rest of the data array. data(j)=0. enddo 12 do13j=1,4 Load the endpoints. endpts(j)=data(j) endpts(j+4)=data(M-3+j) enddo 13 call realft(data,NDFT,1) realft returns the unused value corresponding to ωN/2indata(2) . We actually want this element to contain the imaginary part corresponding to ω0,w h i c hi sz e r o . data(2)=0. endif Now interpolate on the DFT result for the desired frequency. If the frequency is an ωn, i.e., the quantity enis an integer, then cdft=data(2*en-1) ,sdft=data(2*en) ,a n dy o uc o u l d omit the interpolation. en=w*delta*NDFT/TWOPI+1. nn=min(max(int(en-0.5*MPOL+1.),1),NDFT/2-MPOL+1) Leftmost point for the interpola- tion. do14j=1,MPOL cpol(j)=data(2*nn-1) spol(j)=data(2*nn) xpol(j)=nnnn=nn+1 enddo 14 call polint(xpol,cpol,MPOL,en,cdft,cerr) call polint(xpol,spol,MPOL,en,sdft,serr)call dftcor(w,delta,a,b,endpts,corre,corim,corfac) Now get the endpoint cor- rection and the multiplica- tive factor W(θ).cdft=cdft*corfac+corre sdft=sdft*corfac+corimc=delta*cos(w*a) Finally multiply by ∆andexp(iωa). s=delta*sin(w*a) cosint=c*cdft-s*sdft sinint=s*cdft+c*sdftreturn END Sometimes one is interested only in the discrete frequencies ωmof equation (13.9.5), the ones that have integral numbers of periods in the interval [a, b]. For smooth h(t), the value of Itends to be much smaller in magnitude at these ω’s than at values in between, since the integral half-periods tend to cancel precisely. (That iswhy one must oversample forinterpolation to be accurate: I(ω)is oscillatory with small magnitude near the ω m’s.) If you want these ωm’s without messy (and possibly inaccurate) interpolation, you have to set Nto 13.9ComputingFourierIntegralsUsingtheFFT 583Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).a multiple of M(compare equations 13.9.5 and 13.9.12). In the method implemented above, however, Nmust beat least M+1,so thesmallestsuch multipleis 2M, resultingin afactor ∼2 unnecessary computing. Alternatively, one can derive a formula like equation (13.9.13), but with the last sample function hM=h(b)omitted from the DFT, but included entirely in theendpoint correctionfor hM. Then onecan set M=N(anintegerpower of2)and getthe specialfrequenciesofequation (13.9.5)withnoadditional overhead. Themodifiedformulais I(ωm)=∆ eiω ma/braceleftbigg W(θ)[DFT(h0...h M−1)]m +α0(θ)h0+α1(θ)h1+α2(θ)h2+α3(θ)h3 +eiω(b−a)/bracketleftBig A(θ)hM+α* 1(θ)hM−1+α* 2(θ)hM−2+α* 3(θ)hM−3/bracketrightBig/bracerightbigg(13.9.14 ) where θ≡ωm∆andA(θ)is given by A(θ)=−α0(θ)( 13.9.15 ) for the trapezoidal case, or A(θ)=(−6+1 1 θ2)+( 6+ θ2)c o s2 θ 6θ4−iIm[α0(θ)] ≈1 3+1 45θ2−8 945θ4+11 14175θ6−iIm[α0(θ)](13.9.16 ) for the cubic case. Factorslike W(θ)arisenaturallywheneveronecalculatesFouriercoefficientsofsmooth functions, and they are sometimes called attenuation factors [1]. However, the endpoint corrections are equally important in obtaining accurate values of integrals. Narasimhanand Karthikeyan [2]have given a formula that is algebraically equivalent to our trapezoidal formula. However, their formula requires the evaluation of twoFFTs, which is unnecessary. The basic idea used here goes back at least to Filon [3]in 1928 (before the FFT!). He used Simpson’s rule (quadratic interpolation). Since this interpolation is not left-right symmetric,two Fouriertransforms arerequired. An alternative algorithmfor equation (13.9.14) hasbeengiven by Lyness in [4]; for related references, see [5]. To our knowledge, the cubic-order formulas derived here have not previously appeared in the literature. Calculating Fouriertransforms when therange of integrationis (−∞,∞)can be tricky. If the function falls off reasonably quickly at infinity, you can split the integral at a largeenough value of t. For example, the integration to +∞can be written /integraldisplay ∞ aeiωth(t)dt=/integraldisplayb aeiωth(t)dt+/integraldisplay∞ beiωth(t)dt =/integraldisplayb aeiωth(t)dt−h(b)eiωb iω+h/prime(b)eiωb (iω)2−··· (13.9.17 ) The splitting point bmust be chosen large enough that the remaining integral over (b,∞)is small. Successive terms in its asymptotic expansion are found by integrating by parts. Theintegral over (a, b)can be done using dftint. You keep as many terms in the asymptotic expansion as you can easily compute. See [6]for some examples of this idea. More powerful methods, which work well for long-tailed functions but which do not use the FFT,are described in [7-9]. CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), p. 88. [1] Narasimhan, M.S. and Karthi keyan, M. 1984, IEEE Transactions on Antennas & Propagation , vol. 32, pp. 404–408. [2] Filon, L.N.G. 1928, Proceedings of the RoyalSociety of Edinburgh , vol. 49, pp. 38–47. [3] 584 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Giunta,G.andMurli,A.1987, ACMTransactionsonMathematicalSoftware ,vol.13,pp.97–107. [4] Lyness, J.N. 1987, in Numerical Integration , P. Keast and G. Fairweather, eds. (Dordrecht: Reidel). [5] Pantis, G. 1975, Journal of Computational Physics , vol. 17, pp. 229–233. [6] Blakemore, M., Evans, G.A., and Hyslop, J. 1976, Journal of Computational Physics , vol. 22, pp. 352–376. [7] Lyness,J.N., andKaper,T.J.1987, SIAM JournalonScientific andStatistical Computing ,vol.8, pp. 1005–1011. [8] Thakkar,A.J., andSmith,V.H. 1975, ComputerPhysicsCommunications , vol.10,pp.73–79.[9] 13.10 Wavelet Transforms Like thefast Fouriertransform(FFT),thediscrete wavelettransform(DWT)is afast,linearoperationthatoperatesonadatavectorwhoselengthisanintegerpower of two, transforming it into a numerically different vector of the same length. AlsoliketheFFT,thewavelettransformisinvertibleandinfactorthogonal—theinverse transform, when viewed as a big matrix, is simply the transpose of the transform. Both FFT and DWT, therefore, can be viewed as a rotation in function space, from the input space (or time) domain, where the basis functions are the unit vectors e i, or Dirac delta functions in the continuumlimit, to a different domain. For the FFT,this new domain has basis functions that are the familiar sines and cosines. In the wavelet domain, the basis functions are somewhat more complicated and have the fanciful names “mother functions” and “wavelets.” Ofcoursethereareaninfinityofpossiblebasesforfunctionspace,almostallof themuninteresting! Whatmakesthewaveletbasisinterestingisthat, unlikesinesand cosines, individual wavelet functions are quite localized in space; simultaneously, likesines and cosines, individual wavelet functions are quite localized in frequency or (more precisely) characteristic scale. As we will see below, the particular kindof dual localization achieved by wavelets renders large classes of functions and operatorssparse,orsparsetosomehighaccuracy,whentransformedintothewavelet domain. Analogously with the Fourier domain, where a class of computations, likeconvolutions, become computationally fast, there is a large class of computations — those that can take advantage of sparsity — that become computationally fast in the wavelet domain [1]. Unlike sines and cosines, which define a unique Fourier transform, there is not one single unique set of wavelets; in fact, there are infinitely many possiblesets. Roughly, the different sets of wavelets make different trade-offs between how compactly they are localized in space and how smooth they are. (There are further fine distinctions.) Daubechies Wavelet FilterCoefficients A particular set of wavelets is specified by a particular set of numbers, called wavelet filter coefficients . Here, we will largely restrict ourselves to wavelet filters in a class discovered by Daubechies [2]. This class includes members ranging from highlylocalizedtohighlysmooth. Thesimplest(andmostlocalized)member,often calledDAUB4, has onlyfourcoefficients, c0,...,c 3. For the momentwe specialize to this case for ease of notation.