Fast Fourier Transform part 4

Chia sẻ: Dasdsadasd Edwqdqd | Ngày: | Loại File: PDF | Số trang:12

Thêm vào BST

Báo xấu

93
lượt xem 5
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

integer arithmetic modulo some large prime N +1, and the N th root of 1 by the modulo arithmetic equivalent. Strictly speaking, these are not Fourier transforms at all, but the properties are quite similar and computational speed can be far superior. On the other hand, their use is somewhat restricted to quantities like correlations and convolutions since the transform itself is not easily interpretable as a “frequency” spectrum.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Fast Fourier Transform part 4

510 Chapter 12. Fast Fourier Transform integer arithmetic modulo some large prime N +1, and the N th root of 1 by the modulo arithmetic equivalent. Strictly speaking, these are not Fourier transforms at all, but the properties are quite similar and computational speed can be far superior. On the other hand, their use is somewhat restricted to quantities like correlations and convolutions since the transform itself is not easily interpretable as a “frequency” spectrum. visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) CITED REFERENCES AND FURTHER READING: Nussbaumer, H.J. 1982, Fast Fourier Transform and Convolution Algorithms (New York: Springer- Verlag). Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (New York: Academic Press). Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall). [1] Bloomﬁeld, P. 1976, Fourier Analysis of Time Series – An Introduction (New York: Wiley). Van Loan, C. 1992, Computational Frameworks for the Fast Fourier Transform (Philadelphia: S.I.A.M.). Beauchamp, K.G. 1984, Applications of Walsh Functions and Related Functions (New York: Academic Press) [non-Fourier transforms]. Heideman, M.T., Johnson, D.H., and Burris, C.S. 1984, IEEE ASSP Magazine, pp. 14–21 (Oc- tober). 12.3 FFT of Real Functions, Sine and Cosine Transforms It happens frequently that the data whose FFT is desired consist of real-valued samples fj , j = 0 . . . N − 1. To use four1, we put these into a complex array with all imaginary parts set to zero. The resulting transform Fn , n = 0 . . . N − 1 satisﬁes FN−n * = Fn . Since this complex-valued array has real values for F0 and FN/2 , and (N/2) − 1 other independent values F1 . . . FN/2−1 , it has the same 2(N/2 − 1) + 2 = N “degrees of freedom” as the original, real data set. However, the use of the full complex FFT algorithm for real data is inefﬁcient, both in execution time and in storage required. You would think that there is a better way. There are two better ways. The ﬁrst is “mass production”: Pack two separate real functions into the input array in such a way that their individual transforms can be separated from the result. This is implemented in the program twofft below. This may remind you of a one-cent sale, at which you are coerced to purchase two of an item when you only need one. However, remember that for correlations and convolutions the Fourier transforms of two functions are involved, and this is a handy way to do them both at once. The second method is to pack the real input array cleverly, without extra zeros, into a complex array of half its length. One then performs a complex FFT on this shorter length; the trick is then to get the required answer out of the result. This is done in the program realft below.
12.3 FFT of Real Functions, Sine and Cosine Transforms 511 Transform of Two Real Functions Simultaneously First we show how to exploit the symmetry of the transform Fn to handle two real functions at once: Since the input data fj are real, the components of the discrete Fourier transform satisfy visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) FN−n = (Fn )* (12.3.1) where the asterisk denotes complex conjugation. By the same token, the discrete Fourier transform of a purely imaginary set of gj ’s has the opposite symmetry. GN−n = −(Gn )* (12.3.2) Therefore we can take the discrete Fourier transform of two real functions each of length N simultaneously by packing the two data arrays as the real and imaginary parts, respectively, of the complex input array of four1. Then the resulting transform array can be unpacked into two complex arrays with the aid of the two symmetries. Routine twofft works out these ideas. void twofft(float data1[], float data2[], float fft1[], float fft2[], unsigned long n) Given two real input arrays data1[1..n] and data2[1..n], this routine calls four1 and returns two complex output arrays, fft1[1..2n] and fft2[1..2n], each of complex length n (i.e., real length 2*n), which contain the discrete Fourier transforms of the respective data arrays. n MUST be an integer power of 2. { void four1(float data[], unsigned long nn, int isign); unsigned long nn3,nn2,jj,j; float rep,rem,aip,aim; nn3=1+(nn2=2+n+n); for (j=1,jj=2;j
512 Chapter 12. Fast Fourier Transform What about the reverse process? Suppose you have two complex transform arrays, each of which has the symmetry (12.3.1), so that you know that the inverses of both transforms are real functions. Can you invert both in a single FFT? This is even easier than the other direction. Use the fact that the FFT is linear and form the sum of the ﬁrst transform plus i times the second. Invert using four1 with isign = −1. The real and imaginary parts of the resulting complex array are the visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) two desired real functions. FFT of Single Real Function To implement the second method, which allows us to perform the FFT of a single real function without redundancy, we split the data set in half, thereby forming two real arrays of half the size. We can apply the program above to these two, but of course the result will not be the transform of the original data. It will be a schizophrenic combination of two transforms, each of which has half of the information we need. Fortunately, this schizophrenia is treatable. It works like this: The right way to split the original data is to take the even-numbered fj as one data set, and the odd-numbered fj as the other. The beauty of this is that we can take the original real array and treat it as a complex array hj of half the length. The ﬁrst data set is the real part of this array, and the second is the imaginary part, as prescribed for twofft. No repacking is required. In other words hj = f2j + if2j+1 , j = 0, . . . , N/2 − 1. We submit this to four1, and it will give back a complex array Hn = Fn + iFn , n = 0, . . . , N/2 − 1 with e o N/2−1 e Fn = f2k e2πikn/(N/2) k=0 (12.3.3) N/2−1 o Fn = f2k+1 e2πikn/(N/2) k=0 The discussion of program twofft tells you how to separate the two transforms e o Fn and Fn out of Hn . How do you work them into the transform Fn of the original data set fj ? Simply glance back at equation (12.2.3): e o Fn = Fn + e2πin/N Fn n = 0, . . . , N − 1 (12.3.4) Expressed directly in terms of the transform Hn of our real (masquerading as complex) data set, the result is 1 i Fn = (Hn + HN/2−n *) − (Hn − HN/2−n *)e2πin/N n = 0, . . . , N − 1 2 2 (12.3.5) A few remarks: • Since FN−n * = Fn there is no point in saving the entire spectrum. The positive frequency half is sufﬁcient and can be stored in the same array as the original data. The operation can, in fact, be done in place. • Even so, we need values Hn , n = 0, . . . , N/2 whereas four1 gives only the values n = 0, . . . , N/2 − 1. Symmetry to the rescue, HN/2 = H0 .
12.3 FFT of Real Functions, Sine and Cosine Transforms 513 • The values F0 and FN/2 are real and independent. In order to actually get the entire Fn in the original array space, it is convenient to put FN/2 into the imaginary part of F0 . • Despite its complicated form, the process above is invertible. First peel FN/2 out of F0 . Then construct visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) e 1 * Fn = (Fn + FN/2−n ) 2 n = 0, . . . , N/2 − 1 (12.3.6) 1 Fn = e−2πin/N (Fn − FN/2−n ) o * 2 (1) (2) and use four1 to ﬁnd the inverse transform of Hn = Fn + iFn . Surprisingly, the actual algebraic steps are virtually identical to those of the forward transform. Here is a representation of what we have said: #include void realft(float data[], unsigned long n, int isign) Calculates the Fourier transform of a set of n real-valued data points. Replaces this data (which is stored in array data[1..n]) by the positive frequency half of its complex Fourier transform. The real-valued ﬁrst and last components of the complex transform are returned as elements data[1] and data[2], respectively. n must be a power of 2. This routine also calculates the inverse transform of a complex data array if it is the transform of real data. (Result in this case must be multiplied by 2/n.) { void four1(float data[], unsigned long nn, int isign); unsigned long i,i1,i2,i3,i4,np3; float c1=0.5,c2,h1r,h1i,h2r,h2i; double wr,wi,wpr,wpi,wtemp,theta; Double precision for the trigonomet- ric recurrences. theta=3.141592653589793/(double) (n>>1); Initialize the recurrence. if (isign == 1) { c2 = -0.5; four1(data,n>>1,1); The forward transform is here. } else { c2=0.5; Otherwise set up for an inverse trans- theta = -theta; form. } wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0+wpr; wi=wpi; np3=n+3; for (i=2;i>2);i++) { Case i=1 done separately below. i4=1+(i3=np3-(i2=1+(i1=i+i-1))); h1r=c1*(data[i1]+data[i3]); The two separate transforms are sep- h1i=c1*(data[i2]-data[i4]); arated out of data. h2r = -c2*(data[i2]+data[i4]); h2i=c2*(data[i1]-data[i3]); data[i1]=h1r+wr*h2r-wi*h2i; Here they are recombined to form data[i2]=h1i+wr*h2i+wi*h2r; the true transform of the origi- data[i3]=h1r-wr*h2r+wi*h2i; nal real data. data[i4] = -h1i+wr*h2i+wi*h2r; wr=(wtemp=wr)*wpr-wi*wpi+wr; The recurrence. wi=wi*wpr+wtemp*wpi+wi; } if (isign == 1) {
514 Chapter 12. Fast Fourier Transform data[1] = (h1r=data[1])+data[2]; Squeeze the ﬁrst and last data to- data[2] = h1r-data[2]; gether to get them all within the } else { original array. data[1]=c1*((h1r=data[1])+data[2]); data[2]=c1*(h1r-data[2]); four1(data,n>>1,-1); This is the inverse transform for the } case isign=-1. } visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) Fast Sine and Cosine Transforms Among their other uses, the Fourier transforms of functions can be used to solve differential equations (see §19.4). The most common boundary conditions for the solutions are 1) they have the value zero at the boundaries, or 2) their derivatives are zero at the boundaries. In the ﬁrst instance, the natural transform to use is the sine transform, given by N−1 Fk = fj sin(πjk/N ) sine transform (12.3.7) j=1 where fj , j = 0, . . . , N − 1 is the data array, and f0 ≡ 0. At ﬁrst blush this appears to be simply the imaginary part of the discrete Fourier transform. However, the argument of the sine differs by a factor of two from the value that would make this so. The sine transform uses sines only as a complete set of functions in the interval from 0 to 2π, and, as we shall see, the cosine transform uses cosines only. By contrast, the normal FFT uses both sines and cosines, but only half as many of each. (See Figure 12.3.1.) The expression (12.3.7) can be “force-ﬁt” into a form that allows its calculation via the FFT. The idea is to extend the given function rightward past its last tabulated value. We extend the data to twice their length in such a way as to make them an odd function about j = N , with fN = 0, f2N−j ≡ −fj j = 0, . . . , N − 1 (12.3.8) Consider the FFT of this extended function: 2N−1 Fk = fj e2πijk/(2N) (12.3.9) j=0 The half of this sum from j = N to j = 2N − 1 can be rewritten with the substitution j = 2N − j 2N−1 N fj e2πijk/(2N) = f2N−j e2πi(2N−j )k/(2N) j=N j =1 (12.3.10) N−1 =− fj e−2πij k/(2N) j =0
12.3 FFT of Real Functions, Sine and Cosine Transforms 515 1 +1 4 2 (a) 0 5 3 visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) −1 3 2 1 +1 (b) 0 4 5 −1 1 +1 2 3 (c) 0 4 −1 5 0 2π Figure 12.3.1. Basis functions used by the Fourier transform (a), sine transform (b), and cosine transform (c), are plotted. The ﬁrst ﬁve basis functions are shown in each case. (For the Fourier transform, the real and imaginary parts of the basis functions are both shown.) While some basis functions occur in more than one transform, the basis sets are distinct. For example, the sine transform functions labeled (1), (3), (5) are not present in the Fourier basis. Any of the three sets can expand any function in the interval shown; however, the sine or cosine transform best expands functions matching the boundary conditions of the respective basis functions, namely zero function values for sine, zero derivatives for cosine. so that N−1 Fk = fj e2πijk/(2N) − e−2πijk/(2N) j=0 (12.3.11) N−1 = 2i fj sin(πjk/N ) j=0 Thus, up to a factor 2i we get the sine transform from the FFT of the extended function. This method introduces a factor of two inefﬁciency into the computation by extending the data. This inefﬁciency shows up in the FFT output, which has zeros for the real part of every element of the transform. For a one-dimensional problem, the factor of two may be bearable, especially in view of the simplicity of the method. When we work with partial differential equations in two or three dimensions, though, the factor becomes four or eight, so efforts to eliminate the inefﬁciency are well rewarded.
516 Chapter 12. Fast Fourier Transform From the original real data array fj we will construct an auxiliary array yj and apply to it the routine realft. The output will then be used to construct the desired transform. For the sine transform of data fj , j = 1, . . . , N − 1, the auxiliary array is y0 = 0 1 (12.3.12) visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) yj = sin(jπ/N )(fj + fN−j ) + (fj − fN−j ) j = 1, . . . , N − 1 2 This array is of the same dimension as the original. Notice that the ﬁrst term is symmetric about j = N/2 and the second is antisymmetric. Consequently, when realft is applied to yj , the result has real parts Rk and imaginary parts Ik given by N−1 Rk = yj cos(2πjk/N ) j=0 N−1 = (fj + fN−j ) sin(jπ/N ) cos(2πjk/N ) j=1 N−1 = 2fj sin(jπ/N ) cos(2πjk/N ) j=0 N−1 (2k + 1)jπ (2k − 1)jπ = fj sin − sin N N j=0 = F2k+1 − F2k−1 (12.3.13) N−1 Ik = yj sin(2πjk/N ) j=0 N−1 1 = (fj − fN−j ) sin(2πjk/N ) j=1 2 N−1 = fj sin(2πjk/N ) j=0 = F2k (12.3.14) Therefore Fk can be determined as follows: F2k = Ik F2k+1 = F2k−1 + Rk k = 0, . . . , (N/2 − 1) (12.3.15) The even terms of Fk are thus determined very directly. The odd terms require a recursion, the starting point of which follows from setting k = 0 in equation (12.3.15) and using F1 = −F−1 : 1 F1 = R0 (12.3.16) 2 The implementing program is
12.3 FFT of Real Functions, Sine and Cosine Transforms 517 #include void sinft(float y[], int n) Calculates the sine transform of a set of n real-valued data points stored in array y[1..n]. The number n must be a power of 2. On exit y is replaced by its transform. This program, without changes, also calculates the inverse sine transform, but in this case the output array should be multiplied by 2/n. { visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) void realft(float data[], unsigned long n, int isign); int j,n2=n+2; float sum,y1,y2; double theta,wi=0.0,wr=1.0,wpi,wpr,wtemp; Double precision in the trigono- metric recurrences. theta=3.14159265358979/(double) n; Initialize the recurrence. wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); y[1]=0.0; for (j=2;j>1)+1;j++) { wr=(wtemp=wr)*wpr-wi*wpi+wr; Calculate the sine for the auxiliary array. wi=wi*wpr+wtemp*wpi+wi; The cosine is needed to continue the recurrence. y1=wi*(y[j]+y[n2-j]); Construct the auxiliary array. y2=0.5*(y[j]-y[n2-j]); y[j]=y1+y2; Terms j and N − j are related. y[n2-j]=y1-y2; } realft(y,n,1); Transform the auxiliary array. y[1]*=0.5; Initialize the sum used for odd terms below. sum=y[2]=0.0; for (j=1;j
518 Chapter 12. Fast Fourier Transform just twice the cosine transform (12.3.17). Another way of thinking about the formula (12.3.17) is to notice that it is the Chebyshev Gauss-Lobatto quadrature formula (see §4.5), often used in Clenshaw-Curtis adaptive quadrature (see §5.9, equation 5.9.4). Once again the transform can be computed without the factor of two inefﬁciency. In this case the auxiliary function is visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) 1 yj = (fj + fN−j ) − sin(jπ/N )(fj − fN−j ) j = 0, . . . , N − 1 (12.3.19) 2 Instead of equation (12.3.15), realft now gives F2k = Rk F2k+1 = F2k−1 + Ik k = 0, . . . , (N/2 − 1) (12.3.20) The starting value for the recursion for odd k in this case is N−1 1 F1 = (f0 − fN ) + fj cos(jπ/N ) (12.3.21) 2 j=1 This sum does not appear naturally among the Rk and Ik , and so we accumulate it during the generation of the array yj . Once again this transform is its own inverse, and so the following routine works for both directions of the transformation. Note that although this form of the cosine transform has N + 1 input and output values, it passes an array only of length N to realft. #include #define PI 3.141592653589793 void cosft1(float y[], int n) Calculates the cosine transform of a set y[1..n+1] of real-valued data points. The transformed data replace the original data in array y. n must be a power of 2. This program, without changes, also calculates the inverse cosine transform, but in this case the output array should be multiplied by 2/n. { void realft(float data[], unsigned long n, int isign); int j,n2; float sum,y1,y2; double theta,wi=0.0,wpi,wpr,wr=1.0,wtemp; Double precision for the trigonometric recurrences. theta=PI/n; Initialize the recurrence. wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); sum=0.5*(y[1]-y[n+1]); y[1]=0.5*(y[1]+y[n+1]); n2=n+2; for (j=2;j>1);j++) { j=n/2+1 unnecessary since y[n/2+1] unchanged. wr=(wtemp=wr)*wpr-wi*wpi+wr; Carry out the recurrence. wi=wi*wpr+wtemp*wpi+wi; y1=0.5*(y[j]+y[n2-j]); Calculate the auxiliary function. y2=(y[j]-y[n2-j]); y[j]=y1-wi*y2; The values for j and N − j are related. y[n2-j]=y1+wi*y2; sum += wr*y2; Carry along this sum for later use in unfold- } ing the transform.
12.3 FFT of Real Functions, Sine and Cosine Transforms 519 realft(y,n,1); Calculate the transform of the auxiliary func- y[n+1]=y[2]; tion. y[2]=sum; sum is the value of F1 in equation (12.3.21). for (j=4;j
520 Chapter 12. Fast Fourier Transform #include #define PI 3.141592653589793 void cosft2(float y[], int n, int isign) Calculates the “staggered” cosine transform of a set y[1..n] of real-valued data points. The transformed data replace the original data in array y. n must be a power of 2. Set isign to +1 for a transform, and to −1 for an inverse transform. For an inverse transform, the output array should be multiplied by 2/n. visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) { void realft(float data[], unsigned long n, int isign); int i; float sum,sum1,y1,y2,ytemp; double theta,wi=0.0,wi1,wpi,wpr,wr=1.0,wr1,wtemp; Double precision for the trigonometric recurrences. theta=0.5*PI/n; Initialize the recurrences. wr1=cos(theta); wi1=sin(theta); wpr = -2.0*wi1*wi1; wpi=sin(2.0*theta); if (isign == 1) { Forward transform. for (i=1;i=4;i-=2) y[i]=y[i-2]-y[i]; Form diﬀerence of odd terms. y[2]=2.0*ytemp; for (i=3;i
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 fj into the ﬁrst N/2 locations, and the visit 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). readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine- 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) 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 ﬁxed 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 ) deﬁned over the two-dimensional grid 0 ≤ k1 ≤ N1 − 1, 0 ≤ k2 ≤ N2 − 1, we can deﬁne its two-dimensional discrete Fourier transform as a complex function H(n1 , n2 ), deﬁned over the same grid, N2 −1 N1 −1 H(n1 , n2 ) ≡ exp(2πik2 n2 /N2 ) exp(2πik1 n1 /N1 ) h(k1 , k2 ) k2 =0 k1 =0 (12.4.1) By pulling the “subscripts 2” exponential outside of the sum over k1 , or by reversing the order of summation and pulling the “subscripts 1” outside of the sum over k2 , 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(k1 , k2)])