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

For a two-dimensional array, this isequivalent to storing the array by rows. If isign is input as −1, data is replaced by its inversetransform times the product of the lengths of all dimensions.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).Re Imrow of 2N2 float numbers524Chapter 12.Fast Fourier Transform{}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).int idim;unsigned long i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;unsigned long ibit,k1,k2,n,nprev,nrem,ntot;float tempi,tempr;double theta,wi,wpi,wpr,wr,wtemp;Double precision for trigonometric recurrences.for (ntot=1,idim=1;idim<=ndim;idim++)Compute total number of complex valntot *= nn[idim];ues.nprev=1;for (idim=ndim;idim>=1;idim--) {Main loop over the dimensions.n=nn[idim];nrem=ntot/(n*nprev);ip1=nprev << 1;ip2=ip1*n;ip3=ip2*nrem;i2rev=1;for (i2=1;i2<=ip2;i2+=ip1) {This is the bit-reversal section of theif (i2 < i2rev) {routine.for (i1=i2;i1<=i2+ip1-2;i1+=2) {for (i3=i1;i3<=ip3;i3+=ip2) {i3rev=i2rev+i3-i2;SWAP(data[i3],data[i3rev]);SWAP(data[i3+1],data[i3rev+1]);}}}ibit=ip2 >> 1;while (ibit >= ip1 && i2rev > ibit) {i2rev -= ibit;ibit >>= 1;}i2rev += ibit;}ifp1=ip1;Here begins the Danielson-Lanczos secwhile (ifp1 < ip2) {tion of the routine.ifp2=ifp1 << 1;theta=isign*6.28318530717959/(ifp2/ip1);Initialize for the trig.

recurwtemp=sin(0.5*theta);rence.wpr = -2.0*wtemp*wtemp;wpi=sin(theta);wr=1.0;wi=0.0;for (i3=1;i3<=ifp1;i3+=ip1) {for (i1=i3;i1<=i3+ip1-2;i1+=2) {for (i2=i1;i2<=ip3;i2+=ifp2) {k1=i2;Danielson-Lanczos formula:k2=k1+ifp1;tempr=(float)wr*data[k2]-(float)wi*data[k2+1];tempi=(float)wr*data[k2+1]+(float)wi*data[k2];data[k2]=data[k1]-tempr;data[k2+1]=data[k1+1]-tempi;data[k1] += tempr;data[k1+1] += tempi;}}wr=(wtemp=wr)*wpr-wi*wpi+wr; Trigonometric recurrence.wi=wi*wpr+wtemp*wpi+wi;}ifp1=ifp2;}nprev *= n;}52512.5 Fourier Transforms of Real Data in Two and Three DimensionsCITED REFERENCES AND FURTHER READING:Nussbaumer, H.J.

1982, Fast Fourier Transform and Convolution Algorithms (New York: SpringerVerlag).Two-dimensional FFTs are particularly important in the field of image processing. An image is usually represented as a two-dimensional array of pixel intensities,real (and usually positive) numbers. One commonly desires to filter high, or low,frequency spatial components from an image; or to convolve or deconvolve theimage with some instrumental point spread function.

Use of the FFT is by far themost efficient technique.In three dimensions, a common use of the FFT is to solve Poisson’s equationfor a potential (e.g., electromagnetic or gravitational) on a three-dimensional latticethat represents the discretization of three-dimensional space. Here the source terms(mass or charge distribution) and the desired potentials are also real. In two andthree dimensions, with large arrays, memory is often at a premium. It is thereforeimportant to perform the FFTs, insofar as possible, on the data “in place.” Wewant a routine with functionality similar to the multidimensional FFT routine fourn(§12.4), but which operates on real, not complex, input data.

We give such aroutine in this section. The development is analogous to that of §12.3 leading tothe one-dimensional routine realft. (You might wish to review that material atthis point, particularly equation 12.3.5.)It is convenient to think of the independent variables n1 , . . . , nL in equation(12.4.3) as representing an L-dimensional vector ~n in wave-number space, withvalues on the lattice of integers. The transform H(n1 , .

. . , nL ) is then denotedH(~n).It is easy to see that the transform H(~n) is periodic in each of its L dimensions.~2, P~ 3, . . . denote the vectors (N1 , 0, 0, . . .), (0, N2 , 0, . . .),Specifically, if P~1 , P(0, 0, N3 , . . .), and so forth, thenH(~n ± P~j ) = H(~n)j = 1, .

. . , L(12.5.1)Equation (12.5.1) holds for any input data, real or complex. When the data is real,we have the additional symmetryH(−~n) = H(~n)*(12.5.2)Equations (12.5.1) and (12.5.2) imply that the full transform can be trivially obtainedfrom the subset of lattice values ~n that have0 ≤ n1 ≤ N1 − 10 ≤ n2 ≤ N2 − 1(12.5.3)···0 ≤ nL ≤NL2Sample 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).12.5 Fourier Transforms of Real Data in Twoand Three Dimensions.