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 the522 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 machine￾readable 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 L￾dimensions, 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 multidi￾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 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