c13-9 (Numerical Recipes in C), страница 2

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).Here θ ≡ ω∆, and the functions W (θ) and αj (θ) are defined byZ ∞W (θ) ≡ds eiθs ψ(s)−∞Z ∞αj (θ) ≡ds eiθs ϕj (s − j)13.9 Computing Fourier Integrals Using the FFT587Cubic order:W (θ) =6 + θ23θ4(3 − 4 cos θ + cos 2θ) ≈ 1 −11 423 6θ +θ72015120(−42 + 5θ2 ) + (6 + θ2 )(8 cos θ − cos 2θ)(−12θ + 6θ3 ) + (6 + θ2 ) sin 2θ+i46θ6θ41 22 22103 4169 68 4862θ +θ −θ + iθ+θ −θ +θ6≈− +34515120226800451052835467775α0 (θ) =−4(3 − θ2 ) + 2(6 + θ2 ) cos θ−12θ + 2(6 + θ2 ) sin θ+i43θ3θ41 21 215 41711 413θ −θ +θ6 + iθ −+θ −θ +θ6≈− +64560486480090210907207484400α2 (θ) =2(3 − θ2 ) − (6 + θ2 ) cos θ6θ − (6 + θ2 ) sin θ+i46θ6θ41 21 215111137−θ +θ4 −θ6 + iθ−θ +θ4 −θ6≈241802419225920036084036288029937600α3 (θ) =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 thefactor W (θ). The code is turgid, but only because the formulas above are complicated. Theformulas 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-overvalue of θ depends on your machine’s wordlength, but you can always find it experimentallyas 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 theangular frequency w, stepsize delta, lower and upper limits of the integral a and b, while thearray endpts contains the first 4 and last 4 function values. The correction factor W (θ) isreturned as corfac, while the real and imaginary parts of the endpoint correction are returnedas 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.01416e0) nrerror("bad arguments to dftcor");if (fabs(th) < 5.0e-2) {Use series.t=th;t2=t*t;t4=t2*t2;t6=t4*t2;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).14(3 − θ2 ) − 7(6 + θ2 ) cos θ30θ − 5(6 + θ2 ) sin θ+i6θ46θ47 21 275 4711 41376−θ +θ −θ + iθ−θ +θ −θ6≈24180345625920072168725765987520α1 (θ) =588Chapter 13.Fourier and Spectral Applications}Since the use of dftcor can be confusing, we also give an illustrative program dftintwhich uses dftcor to compute equation (13.9.1) for general a, b, ω, and h(t). Several pointswithin this program bear mentioning: The parameters M and NDFT correspond to M and Nin the above discussion.

On successive calls, we recompute the Fourier transform only ifa or b or h(t) has changed.Since dftint is designed to work for any value of ω satisfying ω∆ < π, not just thespecial values returned by the DFT (equation 13.9.12), we do polynomial interpolation ofdegree 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 theendpoint 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 generalpurpose program, because it does not adapt its parameters M, NDFT, MPOL, or its interpolationscheme, to any particular function h(t).

You will have to experiment with your ownapplication.#include <math.h>#include "nrutil.h"#define M 64#define NDFT 1024#define MPOL 6#define TWOPI (2.0*3.14159265)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).*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;13.9 Computing Fourier Integrals Using the FFT589The values of M, NDFT, and MPOL are merely illustrative and should be optimized for yourparticular application.

M is the number of subintervals, NDFT is the length of the FFT (a powerof 2), and MPOL is the degree of polynomial interpolation used to obtain the desired frequencyfrom the FFT.cpol=vector(1,MPOL);spol=vector(1,MPOL);xpol=vector(1,MPOL);if (init != 1 || a != aold || b != bold || func != funcold) {Do we need to initialize, or is only ω changed?init=1;aold=a;bold=b;funcold=func;delta=(b-a)/M;Load the function values into the data array.for (j=1;j<=M+1;j++)data[j]=(*func)(a+(j-1)*delta);for (j=M+2;j<=NDFT;j++)Zero pad the rest of the data array.data[j]=0.0;for (j=1;j<=4;j++) {Load the endpoints.endpts[j]=data[j];endpts[j+4]=data[M-3+j];}realft(data,NDFT,1);realft returns the unused value corresponding to ωN/2 in data[2].

We actually wantthis element to contain the imaginary part corresponding to ω0 , which is zero.data[2]=0.0;}Now interpolate on the DFT result for the desired frequency. If the frequency is an ωn ,i.e., the quantity en is an integer, then cdft=data[2*en-1], sdft=data[2*en], and youcould omit the interpolation.en=w*delta*NDFT/TWOPI+1.0;nn=IMIN(IMAX((int)(en-0.5*MPOL+1.0),1),NDFT/2-MPOL+1); Leftmost point for thefor (j=1;j<=MPOL;j++,nn++) {interpolation.cpol[j]=data[2*nn-1];spol[j]=data[2*nn];xpol[j]=nn;}polint(xpol,cpol,MPOL,en,&cdft,&cerr);polint(xpol,spol,MPOL,en,&sdft,&serr);dftcor(w,delta,a,b,endpts,&corre,&corim,&corfac);Now get the endpoint corcdft *= corfac;rection and the mulsdft *= corfac;tiplicative factor W (θ).cdft += corre;sdft += corim;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).void dftint(float (*func)(float), float a, float b, float w, float *cosint,float *sinint)Example program illustrating how to use the routine dftcor. The user supplies an externalRfunction func that returns the quantity h(t).

The routine then returns ab cos(ωt)h(t) dt asRbcosint and a sin(ωt)h(t) dt as sinint.{void dftcor(float w, float delta, float a, float b, float endpts[],float *corre, float *corim, float *corfac);void polint(float xa[], float ya[], int n, float x, float *y, float *dy);void realft(float data[], unsigned long n, int isign);static int init=0;int j,nn;static float aold = -1.e30,bold = -1.e30,delta,(*funcold)(float);static float data[NDFT+1],endpts[9];float c,cdft,cerr,corfac,corim,corre,en,s;float sdft,serr,*cpol,*spol,*xpol;590Chapter 13.Fourier and Spectral ApplicationsFinally multiply by ∆ and exp(iωa).c=delta*cos(w*a);s=delta*sin(w*a);*cosint=c*cdft-s*sdft;*sinint=s*cdft+c*sdft;free_vector(cpol,1,MPOL);free_vector(spol,1,MPOL);free_vector(xpol,1,MPOL);Sometimes one is interested only in the discrete frequencies ωm of equation (13.9.5),the ones that have integral numbers of periods in the interval [a, b]. For smooth h(t), thevalue of I tends to be much smaller in magnitude at these ω’s than at values in between,since the integral half-periods tend to cancel precisely.

(That is why one must oversample forinterpolation to be accurate: I(ω) is oscillatory with small magnitude near the ωm ’s.) If youwant these ωm ’s without messy (and possibly inaccurate) interpolation, you have to set N toa multiple of M (compare equations 13.9.5 and 13.9.12). In the method implemented above,however, N must be at least M + 1, so the smallest such multiple is 2M , resulting in a factor∼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 inthe endpoint correction for hM .

Then one can set M = N (an integer power of 2) and get thespecial frequencies of equation (13.9.5) with no additional overhead. The modified formula isI(ωm ) = ∆eiωm a W (θ)[DFT(h0 . . . hM −1 )]m+ α0(θ)h0 + α1 (θ)h1 + α2 (θ)h2 + α3 (θ)h3(13.9.14)hi+ eiω(b−a) A(θ)hM + α*1 (θ)hM −1 + α*2 (θ)hM −2 + α*3 (θ)hM −3where θ ≡ ωm ∆ and A(θ) is given byA(θ) = −α0 (θ)(13.9.15)for the trapezoidal case, or(−6 + 11θ2 ) + (6 + θ2 ) cos 2θ− i Im[α0(θ)]6θ418 411 61 2≈ +θ −θ +θ − i Im[α0 (θ)]34594514175A(θ) =(13.9.16)for the cubic case.Factors like W (θ) arise naturally whenever one calculates Fourier coefficients of smoothfunctions, and they are sometimes called attenuation factors [1].