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 P1, P2, P3,... denote the vectors (N1, 0, 0,...), (0, N2, 0,...), (0, 0, N3,...), and so forth, then H(n ± Pj ) = 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)
526 Chapter 12.Fast Fourier Transform In fact,this set of values is overcomplete,because there are additional symmetry relations among the transform values that have nL=0 and nL =NL/2.However these symmetries are complicated and their use becomes extremely confusing. Therefore,we will compute our FFT on the lattice subset of equation(12.5.3), even though this requires a small amount of extra storage for the answer,i.e.,the transform is not quite"in place."(Although an in-place transform is in fact possible. we have found it virtually impossible to explain to any user how to unscramble its output,i.e.,where to find the real and imaginary components of the transform at some particular frequency!) 三 We will implement the multidimensional real Fourier transform for the three dimensional case L=3,with the input data stored as a real,three-dimensional array data[1..nn1][1..nn2][1..nn3].This scheme will allow two-dimensional data to be processed with effectively no loss of efficiency simply by choosing nn1 =1. (Note that it must be the first dimension that is set to 1.The output spectrum comes back packaged,logically at least,as a complex,three-dimensional array that we can call SPEC[1..nn1][1..nn2][1..nn3/2+1](cf.equation 12.5.3).In the first two of its three dimensions,the respective frequency values fi or f2 are stored in wrap- around order,that is with zero frequency in the first index value,the smallest positive 2 frequency in the second index value,the smallest negative frequency in the last index value,and so on(cf.the discussion leading up to routines four1 and fourn).The third of the three dimensions returns only the positive half of the frequency spectrum Figure 12.5.1 shows the logical storage scheme.The returned portion of the complex output spectrum is shown as the unshaded part of the lower figure. 9 The physical,as opposed to logical,packaging of the output spectrum is neces- 三堡色∽ sarily a bit different from the logical packaging,because C does not have a convenient. OF SCIENTIFIC portable mechanism for equivalencing real and complex arrays.The subscript range SPEC[1..nn1][1..nn2][1..nn3/2]is returned in the input array data[1..nn1] [1..nn2][1..nn3],with the correspondence Re(SPEc[i1][i2][i3])=data[i1][i2][2*i3-1] (12.5.4) Im(SPEC [i1][i2][i3])=data[i1][i2][2*i3] The remaining "plane"of values,SPEC[1..nn1][1..nn2][nn3/2+1],is returned in the two-dimensional float array speq[1..nn1][1..2*nn2],with the corre- 有 Numerica 10621 431 spondence Recipes Re(SPEC[i1][i2][nn3/2+1])=speq[i1][2*i2-1] 腿 (12.5.5) North Im(SPEc[i1][i2][nn3/2+1])=speq[i1][2*i2] Note that speq contains frequency components whose third component f3 is at the Nyquist critical frequency +fe.In some applications these values will in fact be ignored or set to zero,since they are intrinsically aliased between positive and negative frequencies. With this much introduction,the implementing procedure,called rlft3,is something of an anticlimax.Look in the innermost loop in the procedure,and you will see equation (12.3.5)implemented on the last transform index.The case of i3=1 is coded separately,to account for the fact that speq is to be filled instead of
526 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). In fact, this set of values is overcomplete, because there are additional symmetry relations among the transform values that have nL = 0 and nL = NL/2. However these symmetries are complicated and their use becomes extremely confusing. Therefore, we will compute our FFT on the lattice subset of equation (12.5.3), even though this requires a small amount of extra storage for the answer, i.e., the transform is not quite “in place.” (Although an in-place transform is in fact possible, we have found it virtually impossible to explain to any user how to unscramble its output, i.e., where to find the real and imaginary components of the transform at some particular frequency!) We will implement the multidimensional real Fourier transform for the three dimensional case L = 3, with the input data stored as a real, three-dimensional array data[1..nn1][1..nn2][1..nn3]. This scheme will allow two-dimensional data to be processed with effectively no loss of efficiency simply by choosing nn1 = 1. (Note that it must be the first dimension that is set to 1.) The output spectrum comes back packaged, logically at least, as a complex, three-dimensional array that we can call SPEC[1..nn1][1..nn2][1..nn3/2+1] (cf. equation 12.5.3). In the first two of its three dimensions, the respective frequency values f 1 or f2 are stored in wraparound order, that is with zero frequency in the first index value, the smallest positive frequency in the second index value, the smallest negative frequency in the last index value, and so on (cf. the discussion leading up to routines four1 and fourn). The third of the three dimensions returns only the positive half of the frequency spectrum. Figure 12.5.1 shows the logical storage scheme. The returned portion of the complex output spectrum is shown as the unshaded part of the lower figure. The physical, as opposed to logical, packaging of the output spectrum is necessarily a bit different from the logical packaging, because C does not have a convenient, portable mechanism for equivalencing real and complex arrays. The subscript range SPEC[1..nn1][1..nn2][1..nn3/2] is returned in the input array data[1..nn1] [1..nn2][1..nn3], with the correspondence Re(SPEC[i1][i2][i3]) = data[i1][i2][2*i3-1] Im(SPEC[i1][i2][i3]) = data[i1][i2][2*i3] (12.5.4) The remaining “plane” of values, SPEC[1..nn1][1..nn2][nn3/2+1], is returned in the two-dimensional float array speq[1..nn1][1..2*nn2],with the correspondence Re(SPEC[i1][i2][nn3/2+1]) = speq[i1][2*i2-1] Im(SPEC[i1][i2][nn3/2+1]) = speq[i1][2*i2] (12.5.5) Note that speq contains frequency components whose third component f 3 is at the Nyquist critical frequency ±fc. In some applications these values will in fact be ignored or set to zero, since they are intrinsically aliased between positive and negative frequencies. With this much introduction, the implementing procedure, called rlft3, is something of an anticlimax. Look in the innermost loop in the procedure, and you will see equation (12.3.5) implemented on the last transform index. The case of i3=1 is coded separately, to account for the fact that speq is to be filled instead of
12.5 Fourier Transforms of Real Data in Two and Three Dimensions 527 nn2,1 nn2,nn3 Input data array float data [1..nnl][1..nn2][1..nn3] Output spectrum arrays (complex) http://ww.nr.com or call 1-800-872- Permission is read able files (including this one) nn2,1 nn2,nn3/2 from NUMERICAL RECIPES IN C: 1,1 1,nn3 -7423(North America 2. 及=6 only), granted for internet users to make one paper copy for their send 州PM电r:to to 1988-1992 by Cambridge University Press.Programs Copyright(C) to dir 1,1 5=0 1,nn3/2 v@cambridge.org 1988-1992 by Numerical Recipes THE ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) 5=- (outside North America) Software. visit website Figure 12.5.1.Input and output data arrangement for rlft3.All arrays shown are presumed to have a first (leftmost)dimension of range [1..nn1],coming out of the page.The input data array is a real,three-dimensional array data[1..nn1][1..nn2][1..nn3].(For two-dimensional data,one sets nn1 1.)The output data can be viewed as a single complex array with dimensions [1..nn1][1..nn2][1..nn3/2+1](cf.equation 12.5.3),corresponding to the frequency components fi and f2 being stored in wrap-around order,but only positive f values being stored (others being obtainable by symmetry).The output data is actually returned mostly in the input array data,but partly stored in the real array speq[1..nn1][1..2*nn2].See text for details
12.5 Fourier Transforms of Real Data in Two and Three Dimensions 527 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). 1,nn3/2 Input data array Output spectrum arrays (complex) nn2, 1 nn2, nn3 float data[1..nn1][1..nn2][1..nn3] 1,1 1, nn3 nn2,nn3/2 f3 = 0 f2 = 0 f 3 = – fc f2 = −fc 1,1 nn2,1 returned in spec[1..nn1][1..2*nn2] returned in data[1..nn1][1..nn2][1..nn3] f2 = fc Figure 12.5.1. Input and output data arrangement for rlft3. All arrays shown are presumed to have a first (leftmost) dimension of range [1..nn1], coming out of the page. The input data array is a real, three-dimensional array data[1..nn1][1..nn2][1..nn3]. (For two-dimensional data, one sets nn1 = 1.) The output data can be viewed as a single complex array with dimensions [1..nn1][1..nn2][1..nn3/2+1] (cf. equation 12.5.3), corresponding to the frequency components f1 and f2 being stored in wrap-around order, but only positive f3 values being stored (others being obtainable by symmetry). The output data is actually returned mostly in the input array data, but partly stored in the real array speq[1..nn1][1..2*nn2]. See text for details
528 Chapter 12.Fast Fourier Transform overwriting the input array of data.The three enclosing for loops (indices 12,i3. and i1,from inside to outside)could in fact be done in any order-their actions all commute.We chose the order shown because of the following considerations:(i)i3 should not be the inner loop,because if it is,then the recurrence relations on wr and wi become burdensome.(ii)On virtual-memory machines,i1 should be the outer loop.because (with C order of array storage)this results in the array data,which might be very large,being accessed in block sequential order. Note that the work done in rlft3 is quite(logarithmically)small,compared to the associated complex FFT.fourn.Since C does not have a convenient complex 5常 type,the operations are carried out explicitly below in terms of real and imaginary 8 parts.The routine rlft3 is based on an earlier routine by G.B.Rybicki. #include <math.h> nted for void rlft3(float ***data,float **speq,unsigned long nn1,unsigned long nn2, unsigned long nn3,int isign) Given a three-dimensional real array data[1..nn1][1..nn2][1..nn3](where nn1 =1 for the case of a logically two-dimensional array),this routine returns (for isign=1)the complex ERICI fast Fourier transform as two complex arrays:On output,data contains the zero and positive frequency values of the third frequency component,while speq [1..nn1][1..2*nn2]contains the Nyquist critical frequency values of the third frequency component.First (and second) 令 frequency components are stored for zero,positive,and negative frequencies,in standard wrap- around order.See text for description of how complex values are arranged.For isign=-1,the inverse transform (times nn1*nn2*nn3/2 as a constant multiplicative factor)is performed, Press. with output data (viewed as a real array)deriving from input data (viewed as complex)and ART speq.For inverse transforms on data not generated first by a forward transform,make sure the complex input data array satisfies property (12.5.2).The dimensions nn1,nn2,nn3 must Programs always be integer powers of 2. void fourn(float data[],unsigned long nn[],int ndim,int isign); SCIENTIFIC void nrerror(char error_text[]); unsigned long i1,12,i3,j1,j2,j3,nn[4],ii3; double theta,wi,wpi,wpr,wr,wtempi float c1,c2,hir,hii,h2r,h2i; if (1+&data [nn1][nn2][nn3]-&data[1][1][1]!nn1*nn2*nn3) nrerror("rlft3:problem with dimensions or contiguity of data array\n"); 19881992 COMPUTING(ISBN c1=0.5; c2=-0.5*1s1gn; thetas=is1gn*(6.28318530717959/nn3); wtemp=sin(0.5*theta); wpr =-2.0*wtemp*wtemp; Numerical Recipes -43108 wpi=sin(theta); nn[1]=nn1: nn[2]=nn2; nn[3]=nn3>>1; (outside 膜 1f(1s1gm=1)[ Case of forward transform. oftware. fourn(&data[1][1][1]-1,nn,3,isign); Here is where most all of the com- North for(11=1;11<=nn1;11++) pute time is spent. for(12=1,j2=0;i2<nn2;12+)[ Extend data periodically into speq. Ame speg[i1][+j2]=data[i1][i2][1]: speq[i1][+j2]=data[i1][i2][2]; for(11=1;11<=nn1;11++){ j1=(11!=1？nn1-11+2:1); Zero frequency is its own reflection,otherwise locate corresponding negative frequency in wrap-around order. wr=1.0; Initialize trigonometric recurrence. w1=0.0; for(113=1,13=1;13<=(nn3>>2)+1;13++,1i3+=2)
528 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). overwriting the input array of data. The three enclosing for loops (indices i2, i3, and i1, from inside to outside) could in fact be done in any order — their actions all commute. We chose the order shown because of the following considerations: (i) i3 should not be the inner loop, because if it is, then the recurrence relations on wr and wi become burdensome. (ii) On virtual-memory machines, i1 should be the outer loop, because (with C order of array storage) this results in the array data, which might be very large, being accessed in block sequential order. Note that the work done in rlft3 is quite (logarithmically) small, compared to the associated complex FFT, fourn. Since C does not have a convenient complex type, the operations are carried out explicitly below in terms of real and imaginary parts. The routine rlft3 is based on an earlier routine by G.B. Rybicki. #include <math.h> void rlft3(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, int isign) Given a three-dimensional real array data[1..nn1][1..nn2][1..nn3] (where nn1 = 1 for the case of a logically two-dimensional array), this routine returns (for isign=1) the complex fast Fourier transform as two complex arrays: On output, data contains the zero and positive frequency values of the third frequency component, while speq[1..nn1][1..2*nn2] contains the Nyquist critical frequency values of the third frequency component. First (and second) frequency components are stored for zero, positive, and negative frequencies, in standard wraparound order. See text for description of how complex values are arranged. For isign=-1, the inverse transform (times nn1*nn2*nn3/2 as a constant multiplicative factor) is performed, with output data (viewed as a real array) deriving from input data (viewed as complex) and speq. For inverse transforms on data not generated first by a forward transform, make sure the complex input data array satisfies property (12.5.2). The dimensions nn1, nn2, nn3 must always be integer powers of 2. { void fourn(float data[], unsigned long nn[], int ndim, int isign); void nrerror(char error_text[]); unsigned long i1,i2,i3,j1,j2,j3,nn[4],ii3; double theta,wi,wpi,wpr,wr,wtemp; float c1,c2,h1r,h1i,h2r,h2i; if (1+&data[nn1][nn2][nn3]-&data[1][1][1] != nn1*nn2*nn3) nrerror("rlft3: problem with dimensions or contiguity of data array\n"); c1=0.5; c2 = -0.5*isign; theta=isign*(6.28318530717959/nn3); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); nn[1]=nn1; nn[2]=nn2; nn[3]=nn3 >> 1; if (isign == 1) { Case of forward transform. fourn(&data[1][1][1]-1,nn,3,isign); Here is where most all of the comfor (i1=1;i1<=nn1;i1++) pute time is spent. for (i2=1,j2=0;i2<=nn2;i2++) { Extend data periodically into speq. speq[i1][++j2]=data[i1][i2][1]; speq[i1][++j2]=data[i1][i2][2]; } } for (i1=1;i1<=nn1;i1++) { j1=(i1 != 1 ? nn1-i1+2 : 1); Zero frequency is its own reflection, otherwise locate corresponding negative frequency in wrap-around order. wr=1.0; Initialize trigonometric recurrence. wi=0.0; for (ii3=1,i3=1;i3<=(nn3>>2)+1;i3++,ii3+=2) {
12.5 Fourier Transforms of Real Data in Two and Three Dimensions 529 http://www.nr. read able files .com or call 1-800 (including this one) granted fori internet to any Figure 12.5.2.(a)A two-dimensional image with intensities either purely black or purely white.(b)The tusers to make one paper 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN C: same image,after it has been low-pass filtered using rlft3.Regions with fine-scale features become gray. for(12=1:12<=nn2;12++){ 1f(13==1) Equation (12.3.5). (North America server computer, THE j2=(121=1？((nn2-12)<1)+3:1); ART h1r=c1*(data[11][i2][1]+speq[1][j2]): h1i=c1*(data[i1][i2][2]-speq[j1][j2+1]); 是 h2i=c2*(data[i1][i2][1]-speq[j1][j2]); copy for thei Programs h2r=-c2*(data[11][i2][2]+speq[j1][j2+1]): data[i1][i2][1]-hir+h2r; send data[11][i2][2]=h1i+h2i; speq[j1][j2]=hir-h2r; Copyright (C) speq[j1][j2+1]=h2i-h1i; else j2=(121=17nn2-12+2:1); j3=nn3+3-(13<<1); email to directcustsen hir=c1*(data[i1][i2][ii3]+data[j1][j2][j3]); h1i=c1*(data[i1][i2][ii3+1]-data[j1][j2][j3+1]); OF SCIENTIFIC COMPUTING(ISBN 0-521- h21=c2*(data[i1][i2][ii3]-data[j1][j2][j3]); h2r=-c2*(data[i1][i2][ii3+1]+data[j1][j2][j3+1]); data[i1][i2][ii3]=hir+wr*h2r-wi*h2i; data[i1][i2][ii3+1]-h1i+wr*h2i+wi*h2r; data[j1][j2][j3]=hir-wr*h2r+wi*h2i; @cambridge.org 1988-1992 by Numerical Recipes data[j1][j2][j3+1]=-h1i+wr*h2i+wi*h2r; -431085 wr=(wtemp=wr)*wpr-wi*wpi+wr; Do the recurrence. wi=wi*wpr+wtemp*wpi+wi; Software. d 1f(1s1gn=-1) (outside North America). Case of reverse transform. fourn(&data[1][1][1]-1,nn,3,isign); visit website machine We now give some fragments from notional calling programs,to clarify the use of rlft3 for two-and three-dimensional data.Note again that the routine does not actually distinguish between two and three dimensions;two is treated like three,but with the first dimension having length 1.Since the first dimension is the outer loop. virtually no inefficiency is introduced. The first program fragment FFTs a two-dimensional data array,allows for some
12.5 Fourier Transforms of Real Data in Two and Three Dimensions 529 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). Figure 12.5.2. (a) A two-dimensional image with intensities either purely black or purely white. (b) The same image, after it has been low-pass filtered using rlft3. Regions with fine-scale features become gray. for (i2=1;i2<=nn2;i2++) { if (i3 == 1) { Equation (12.3.5). j2=(i2 != 1 ? ((nn2-i2)<<1)+3 : 1); h1r=c1*(data[i1][i2][1]+speq[j1][j2]); h1i=c1*(data[i1][i2][2]-speq[j1][j2+1]); h2i=c2*(data[i1][i2][1]-speq[j1][j2]); h2r= -c2*(data[i1][i2][2]+speq[j1][j2+1]); data[i1][i2][1]=h1r+h2r; data[i1][i2][2]=h1i+h2i; speq[j1][j2]=h1r-h2r; speq[j1][j2+1]=h2i-h1i; } else { j2=(i2 != 1 ? nn2-i2+2 : 1); j3=nn3+3-(i3<<1); h1r=c1*(data[i1][i2][ii3]+data[j1][j2][j3]); h1i=c1*(data[i1][i2][ii3+1]-data[j1][j2][j3+1]); h2i=c2*(data[i1][i2][ii3]-data[j1][j2][j3]); h2r= -c2*(data[i1][i2][ii3+1]+data[j1][j2][j3+1]); data[i1][i2][ii3]=h1r+wr*h2r-wi*h2i; data[i1][i2][ii3+1]=h1i+wr*h2i+wi*h2r; data[j1][j2][j3]=h1r-wr*h2r+wi*h2i; data[j1][j2][j3+1]= -h1i+wr*h2i+wi*h2r; } } wr=(wtemp=wr)*wpr-wi*wpi+wr; Do the recurrence. wi=wi*wpr+wtemp*wpi+wi; } } if (isign == -1) Case of reverse transform. fourn(&data[1][1][1]-1,nn,3,isign); } We now give some fragments from notional calling programs, to clarify the use of rlft3 for two- and three-dimensional data. Note again that the routine does not actually distinguish between two and three dimensions; two is treated like three, but with the first dimension having length 1. Since the first dimension is the outer loop, virtually no inefficiency is introduced. The first program fragment FFTs a two-dimensional data array, allows for some