12.4 FFT in Two or More Dimensions 521 An alternative way of implementing this algorithm is to form an auxiliary function by copying the even elements of fi into the first N/2 locations,and the odd elements into the next N/2 elements in reverse order.However,it is not easy to implement the alternative algorithm without a temporary storage array and we prefer the above in-place algorithm. Finally,we mention that there exist fast cosine transforms for small N that do not rely on an auxiliary function or use an FFT routine.Instead,they carry out the transform directly,often coded in hardware for fixed N of small dimension [1]. CITED REFERENCES AND FURTHER READING: Brigham,E.O.1974,The Fast Fourier Transform(Englewood Cliffs,NJ:Prentice-Hall),$10-10. Sorensen,H.V..Jones,D.L.,Heideman,M.T.,and Burris,C.S.1987,/EEE Transactions on Acoustics,Speech,and Signal Processing,vol.ASSP-35,pp.849-863. RECIPES I Hou,H.S.1987,IEEE Transactions on Acoustics,Speech,and Signal Processing,vol.ASSP-35, pp.1455-1461 [see for additional references]. Hockney,R.W.1971,in Methods in Computational Physics,vol.9(New York:Academic Press). > 令 Press. Temperton,C.1980.Journal of Computational Physics,vol.34,pp.314-329. Clarke,R.J.1985,Transform Coding of Images,(Reading,MA:Addison-Wesley). Gonzalez,R.C.,and Wintz,P.1987,Digital /mage Processing,(Reading.MA:Addison-Wesley). 9 Chen,W.,Smith,C.H.,and Fralick,S.C.1977,IEEE Transactions on Communications,vol.COM- 25,Pp.1004-1009.[1] 5会 OF SCIENTIFIC( 12.4 FFT in Two or More Dimensions COMPUTING (ISBN 188810920 Given a complex function h(,k2)defined over the two-dimensional grid 0≤k1≤Ni-l,0≤2≤N2-l,we can define its two-dimensional discrete Recipes 10.621 Fourier transform as a complex function H(n1,n2),defined over the same grid, Numerical Recipes 43106 (outside H(n1,n2) exp(2rik2n2/N2)exp(2rikin1/N1)h(k1:k2) k2=0k1=0 North (12.4.1) By pulling the"subscripts 2"exponential outside of the sum overk1,or by reversing the order of summation and pulling the"subscripts 1"outside of the sum over 2, we can see instantly that the two-dimensional FFT can be computed by taking one- dimensional FFTs sequentially on each index of the original function.Symbolically, H(n1,n2)=FFT-on-index-1(FFT-on-index-2 [h(k1,k2)]) (12.4.2) FFT-on-index-2(FFT-on-index-1 [h(1,k2)])
12.4 FFT in Two or More Dimensions 521 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). } } } An alternative way of implementing this algorithm is to form an auxiliary function by copying the even elements of fj into the first N/2 locations, and the odd elements into the next N/2 elements in reverse order. However, it is not easy to implement the alternative algorithm without a temporary storage array and we prefer the above in-place algorithm. Finally, we mention that there exist fast cosine transforms for small N that do not rely on an auxiliary function or use an FFT routine. Instead, they carry out the transform directly, often coded in hardware for fixed N of small dimension [1]. CITED REFERENCES AND FURTHER READING: Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall), §10–10. Sorensen, H.V., Jones, D.L., Heideman, M.T., and Burris, C.S. 1987, IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, pp. 849–863. Hou, H.S. 1987, IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-35, pp. 1455–1461 [see for additional references]. Hockney, R.W. 1971, in Methods in Computational Physics, vol. 9 (New York: Academic Press). Temperton, C. 1980, Journal of Computational Physics, vol. 34, pp. 314–329. Clarke, R.J. 1985, Transform Coding of Images, (Reading, MA: Addison-Wesley). Gonzalez, R.C., and Wintz, P. 1987, Digital Image Processing, (Reading, MA: Addison-Wesley). Chen, W., Smith, C.H., and Fralick, S.C. 1977, IEEE Transactions on Communications, vol. COM- 25, pp. 1004–1009. [1] 12.4 FFT in Two or More Dimensions Given a complex function h(k1, k2) defined over the two-dimensional grid 0 ≤ k1 ≤ N1 − 1, 0 ≤ k2 ≤ N2 − 1, we can define its two-dimensional discrete Fourier transform as a complex function H(n1, n2), defined over the same grid, H(n1, n2) ≡ N 2−1 k2=0 N 1−1 k1=0 exp(2πik2n2/N2) exp(2πik1n1/N1) h(k1, k2) (12.4.1) By pulling the “subscripts 2” exponential outside of the sum over k 1, or by reversing the order of summation and pulling the “subscripts 1” outside of the sum over k 2, we can see instantly that the two-dimensional FFT can be computed by taking onedimensional FFTs sequentially on each index of the original function. Symbolically, H(n1, n2) = FFT-on-index-1(FFT-on-index-2[h(k1, k2)]) = FFT-on-index-2(FFT-on-index-1[h(k1, k2)]) (12.4.2)
522 Chapter 12.Fast Fourier Transform For this to be practical,of course,both Ni and N2 should be some efficient length for an FFT,usually a power of 2.Programming a two-dimensional FFT,using (12.4.2)with a one-dimensional FFT routine,is a bit clumsier than it seems at first. Because the one-dimensional routine requires that its input be in consecutive order as a one-dimensional complex array,you find that you are endlessly copying things out of the multidimensional input array and then copying things back into it.This is not recommended technique.Rather,you should use a multidimensional FFT routine,such as the one we give below. The generalization of (12.4.1)to more than two dimensions,say to L- dimensions,is evidently NL-1N1-1 H(n1;...,nL …∑ exp(2 nikLnL/Nz)×· 虽 k=0k1=0 (12.4.3) S x exp(2rikin/N)h(k1,...,kL) where ni and ki range from 0 to N1-1,...,nL and kL range from 0 to NL -1. How many calls to a one-dimensional FFT are in(12.4.3)?Quite a few!For each value of k1,k2,...,kL-1 you FFT to transform the L index.Then for each value of 9 k1,k2,...,kL-2 and nL you FFT to transform the L-1 index.And so on.It is best to rely on someone else having done the bookkeeping for once and for all. The inverse transforms of(12.4.1)or(12.4.3)are just what you would expect them to be:Change the i's in the exponentials to -i's,and put an overall factor of 1/(N1×·×Nz)in front of the whole thing.Most other features 、之员之等今 9 of multidimensional FFTs are also analogous to features already discussed in the one-dimensional case: Frequencies are arranged in wrap-around order in the transform,but now for each separate dimension. The input data are also treated as if they were wrapped around.If they are discontinuous across this periodic identification(in any dimension)then the spectrum will have some excess power at high frequencies because of the discontinuity.The fix,if you care,is to remove multidimensional linear trends. 10621 Ifyou are doing spatial filtering and are worried about wrap-around effects, then you need to zero-pad all around the border of the multidimensional 、p Numerica 431 array.However,be sure to notice how costly zero-padding is in multidi- (outside Recipes mensional transforms.If you use too thick a zero-pad,you are going to waste a lot of storage,especially in 3 or more dimensions! Aliasing occurs as always if sufficient bandwidth limiting does not exist North along one or more of the dimensions of the transform. The routine fourn that we furnish herewith is a descendant of one written by N. M.Brenner.It requires as input(i)a scalar,telling the number of dimensions,e.g., 2;(ii)a vector,telling the length of the array in each dimension,e.g.,(32,64).Note that these lengths must all be powers of 2,and are the numbers of complex values in each direction;(iii)the usual scalar equal to+1 indicating whether you want the transform or its inverse;and,finally (iv)the array of data. A few words about the data array:fourn accesses it as a one-dimensional array of real numbers,that is,data[1..(2NiN2...NL)],of length equal to twice the
522 Chapter 12. Fast Fourier Transform 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). For this to be practical, of course, both N1 and N2 should be some efficient length for an FFT, usually a power of 2. Programming a two-dimensional FFT, using (12.4.2) with a one-dimensional FFT routine, is a bit clumsier than it seems at first. Because the one-dimensional routine requires that its input be in consecutive order as a one-dimensional complex array, you find that you are endlessly copying things out of the multidimensional input array and then copying things back into it. This is not recommended technique. Rather, you should use a multidimensional FFT routine, such as the one we give below. The generalization of (12.4.1) to more than two dimensions, say to Ldimensions, is evidently H(n1,...,nL) ≡ N L−1 kL=0 ··· N 1−1 k1=0 exp(2πikLnL/NL) ×··· × exp(2πik1n1/N1) h(k1,...,kL) (12.4.3) where n1 and k1 range from 0 to N1 − 1, ... , nL and kL range from 0 to NL − 1. How many calls to a one-dimensional FFT are in (12.4.3)? Quite a few! For each value of k1, k2,...,kL−1 you FFT to transform the L index. Then for each value of k1, k2,...,kL−2 and nL you FFT to transform the L − 1 index. And so on. It is best to rely on someone else having done the bookkeeping for once and for all. The inverse transforms of (12.4.1) or (12.4.3) are just what you would expect them to be: Change the i’s in the exponentials to −i’s, and put an overall factor of 1/(N1 ×···× NL) in front of the whole thing. Most other features of multidimensional FFTs are also analogous to features already discussed in the one-dimensional case: • Frequencies are arranged in wrap-around order in the transform, but now for each separate dimension. • The input data are also treated as if they were wrapped around. If they are discontinuous across this periodic identification (in any dimension) then the spectrum will have some excess power at high frequencies because of the discontinuity. The fix, if you care, is to remove multidimensional linear trends. • If you are doing spatial filtering and are worried about wrap-around effects, then you need to zero-pad all around the border of the multidimensional array. However, be sure to notice how costly zero-padding is in multidimensional transforms. If you use too thick a zero-pad, you are going to waste a lot of storage, especially in 3 or more dimensions! • Aliasing occurs as always if sufficient bandwidth limiting does not exist along one or more of the dimensions of the transform. The routine fourn that we furnish herewith is a descendant of one written by N. M. Brenner. It requires as input (i) a scalar, telling the number of dimensions, e.g., 2; (ii) a vector, telling the length of the array in each dimension, e.g., (32,64). Note that these lengths must all be powers of 2, and are the numbers of complex values in each direction; (iii) the usual scalar equal to ±1 indicating whether you want the transform or its inverse; and, finally (iv) the array of data. A few words about the data array: fourn accesses it as a one-dimensional array of real numbers, that is, data[1..(2N1N2 ...NL)], of length equal to twice the
12.4 FFT in Two or More Dimensions 523 -data[1] row of 2N,float numbers row I 万=0 row 2 fi2N△1 y2N1-1 row N/2 f1=N△1 8 row N /2+1 1=±2△ 虽 三g 2N1-1 row N /2+2 =- N△1 row N =-N 令 data [2N N]--- America Press. ART Figure 12.4.1.Storage arrangement of frequencies in the output H(,f2)of a two-dimensional FFT. 9 The input data is a two-dimensional M x N2 array h(t1,t2)(stored by rows of complex numbers). 9 The output is also stored by complex rows.Each row corresponds to a particular value of fi,as shown Progra in the figure.Within each row,the arrangement of frequencies is exactly as shown in Figure 12.2.2. A and A2 are the sampling intervals in the I and 2 directions,respectively.The total number of(real) CIENTIFIC array elements is 2NI N2.The program fourn can also do more than two dimensions,and the storage arrangement generalizes in the obvious way. product of the lengths of the L dimensions.It assumes that the array represents an L-dimensional complex array,with individual components ordered as follows: (i)each complex value occupies two sequential locations,real part followed by imaginary;(ii)the first subscript changes least rapidly as one goes through the array; 10-621 the last subscript changes most rapidly (that is,"store by rows,"the C norm):(iii) subscripts range from I to their maximum values (N1,N2,...,NL,respectively), 43106 rather than from 0 to N-1,N2-1,...,NL-1.Almost all failures to get fourn uction to work result from improper understanding of the above ordering of the data array, so take care!(Figure 12.4.I illustrates the format of the output array.) (outside North Software. #include #define SWAP(a,b)tempr=(a);(a)=(b);(b)-tempr B void fourn(float data[],unsigned long nn[],int ndim,int isign) Replaces data by its ndim-dimensional discrete Fourier transform,if isign is input as 1. nn [1..ndim]is an integer array containing the lengths of each dimension (number of complex values),which MUST all be powers of 2.data is a real array of length twice the product of these lengths,in which the data are stored as in a multidimensional complex array:real and imaginary parts of each element are in consecutive locations,and the rightmost index of the array increases most rapidly as one proceeds along data.For a two-dimensional array,this is equivalent to storing the array by rows.If isign is input as-1,data is replaced by its inverse transform times the product of the lengths of all dimensions
12.4 FFT in Two or More Dimensions 523 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). row of 2N2 float numbers row 1 row 2 row N1/ 2 row N1/ 2 + 1 row N1/ 2 + 2 row N1 1 N1∆1 f1 = 0 f1 = f1 = f1 = − 1⁄ 2N1 − 1 N1∆1 1 N1∆1 f1 = ± 1 2∆1 f1 = − 1⁄ 2N1 − 1 N1∆1 data [1] Re Im data [2N1N2 ] Figure 12.4.1. Storage arrangement of frequencies in the output H(f1, f2) of a two-dimensional FFT. The input data is a two-dimensional N1 × N2 array h(t1, t2) (stored by rows of complex numbers). The output is also stored by complex rows. Each row corresponds to a particular value of f1, as shown in the figure. Within each row, the arrangement of frequencies f2 is exactly as shown in Figure 12.2.2. ∆1 and ∆2 are the sampling intervals in the 1 and 2 directions, respectively. The total number of (real) array elements is 2N1N2. The program fourn can also do more than two dimensions, and the storage arrangement generalizes in the obvious way. product of the lengths of the L dimensions. It assumes that the array represents an L-dimensional complex array, with individual components ordered as follows: (i) each complex value occupies two sequential locations, real part followed by imaginary; (ii) the first subscript changes least rapidly as one goes through the array; the last subscript changes most rapidly (that is, “store by rows,” the C norm); (iii) subscripts range from 1 to their maximum values (N1, N2,...,NL, respectively), rather than from 0 to N1 − 1, N2 − 1,..., NL − 1. Almost all failures to get fourn to work result from improper understanding of the above ordering of the data array, so take care! (Figure 12.4.1 illustrates the format of the output array.) #include #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void fourn(float data[], unsigned long nn[], int ndim, int isign) Replaces data by its ndim-dimensional discrete Fourier transform, if isign is input as 1. nn[1..ndim] is an integer array containing the lengths of each dimension (number of complex values), which MUST all be powers of 2. data is a real array of length twice the product of these lengths, in which the data are stored as in a multidimensional complex array: real and imaginary parts of each element are in consecutive locations, and the rightmost index of the array increases most rapidly as one proceeds along data. For a two-dimensional array, this is equivalent to storing the array by rows. If isign is input as −1, data is replaced by its inverse transform times the product of the lengths of all dimensions
524 Chapter 12. Fast Fourier Transform int idim; unsigned long i1,12,13,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 recur- rences for (ntot=1,idim=1;idim=1;idim--){ Main loop over the dimensions. n=nn[idim]; nrem=ntot/(n*nprev); ipi=nprev 1; while (ibit >ip1 &i2rev ibit){ i2rev -ibit; 是 1b1t>>=1; send i2rev +ibit; ifpl=ipl; Here begins the Danielson-Lanczos sec- Copyright (C) while (ifp1 ip2){ tion of the routine. to dir ifp2=ifpl <1; theta=1s1g即*6.28318530717959/(1fp2/1p1); Initialize for the trig.recur- wtemp=sin(0.5*theta); rence. wpr =-2.0*wtemp*wtemp; ectcustser ART OF SCIENTIFIC COMPUTING(ISBN 0-521 wpi=sin(theta); r=1.0; w1=0.0; for(13=1;13<=1fp1;13+=1p1){ for(11=13;11<=13+ip1-2;11+=2){ for(12=i1;i2<=ip3;i2+=1fp2) Numerical Recipes books or 1988-1992 by Numerical Recipes k1=12: Danielson-Lanczos formula: -431085 k2=k1+1fp1; tempr=(float)wr*data[k2]-(float)wi*data[k2+1]; tempi=(float)wr*data[k2+1]+(float)wi*data[k2]; data [k2]=data[ki]-tempr; data[k2+1]=data[k1+1]-tempi; data[k1]+tempr; (outside North America) Software. data[k1+1]+tempi; wr=(wtemp=wr)*wpr-wi*wpi+wr; Trigonometric recurrence wi=wi*wpr+wtemp*wpi+wi; ifpi=ifp2; nprev *n;
524 Chapter 12. Fast Fourier Transform 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). { 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=1;idim--) { Main loop over the dimensions. n=nn[idim]; nrem=ntot/(n*nprev); ip1=nprev > 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; } }
12.5 Fourier Transforms of Real Data in Two and Three Dimensions 525 CITED REFERENCES AND FURTHER READING: Nussbaumer,H.J.1982.Fast Fourier Transform and Convolution Algorithms(New York:Springer- Verlag) 12.5 Fourier Transforms of Real Data in Two and Three Dimensions Two-dimensional FFTs are particularly important in the field of image process- ing.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 the nted for image with some instrumental point spread function.Use of the FFT is by far the most efficient technique. In three dimensions.a common use of the FFT is to solve Poisson's equation for a potential (e.g.,electromagnetic or gravitational)on a three-dimensional lattice that represents the discretization of three-dimensional space.Here the source terms (mass or charge distribution)and the desired potentials are also real.In two and three dimensions,with large arrays,memory is often at a premium.It is therefore important to perform the FFTs,insofar as possible,on the data"in place."We 、3&ad旧日 Press. want 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 a g% 9 routine in this section.The development is analogous to that of $12.3 leading to the one-dimensional routine realft.(You might wish to review that material at this point,particularly equation 12.3.5.) OF SCIENTIFIC( It is convenient to think of the independent variables n1,...,nL in equation 61 (12.4.3)as representing an L-dimensional vector n in wave-number space,with values on the lattice of integers.The transform H(n1,...,nL)is then denoted H(). It is easy to see that the transform H(n)is periodic in each of its L dimensions. Specifically,if B,P2,P3,...denote the vectors (N1,0,0,...),(0,N2,0,...), (0,0,N3,...),and so forth,then P Numerica 10621 H(元±P)=H()j=1,,L (12.5.1) 43106 Equation(12.5.1)holds for any input data,real or complex.When the data is real, we have the additional symmetry (outside H(-)=H()* (12.5.2) North Equations(12.5.1)and(12.5.2)imply that the full transform can be trivially obtained from the subset of lattice values n that have 0≤n1≤N1-1 0≤n2≤N2-1 EE (12.5.3) L 0<nL≤2
12.5 Fourier Transforms of Real Data in Two and Three Dimensions 525 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). CITED REFERENCES AND FURTHER READING: Nussbaumer, H.J. 1982, Fast Fourier Transform and Convolution Algorithms (New York: SpringerVerlag). 12.5 Fourier Transforms of Real Data in Two and Three Dimensions 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 the image with some instrumental point spread function. Use of the FFT is by far the most efficient technique. In three dimensions, a common use of the FFT is to solve Poisson’s equation for a potential (e.g., electromagnetic or gravitational) on a three-dimensional lattice that represents the discretization of three-dimensional space. Here the source terms (mass or charge distribution) and the desired potentials are also real. In two and three dimensions, with large arrays, memory is often at a premium. It is therefore important to perform the FFTs, insofar as possible, on the data “in place.” We want 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 a routine in this section. The development is analogous to that of §12.3 leading to the one-dimensional routine realft. (You might wish to review that material at this 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, with values on the lattice of integers. The transform H(n1,...,nL) is then denoted H( n). It is easy to see that the transform H( n) is periodic in each of its L dimensions. Specifically, if P 1, P 2, P 3,... denote the vectors (N1, 0, 0,...), (0, N2, 0,...), (0, 0, N3,...), and so forth, then H( 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 symmetry H(− n) = H( n)* (12.5.2) Equations (12.5.1) and (12.5.2) imply that the full transform can be trivially obtained from the subset of lattice values n that have 0 ≤ n1 ≤ N1 − 1 0 ≤ n2 ≤ N2 − 1 ··· 0 ≤ nL ≤ NL 2 (12.5.3)