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

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).void cosft2(float y[], int n, int isign)Calculates the “staggered” cosine transform of a set y[1..n] of real-valued data points. Thetransformed 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 outputarray should be multiplied by 2/n.{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.12.4 FFT in Two or More Dimensions521}}}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 onAcoustics, 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. COM25, pp. 1004–1009. [1]12.4 FFT in Two or More DimensionsGiven a complex function h(k1 , k2 ) defined over the two-dimensional grid0 ≤ k1 ≤ N1 − 1, 0 ≤ k2 ≤ N2 − 1, we can define its two-dimensional discreteFourier transform as a complex function H(n1 , n2 ), defined over the same grid,H(n1 , n2 ) ≡NX2 −1 N1 −1Xexp(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 reversingthe 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 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)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).An alternative way of implementing this algorithm is to form an auxiliaryfunction by copying the even elements of fj into the first N/2 locations, and theodd elements into the next N/2 elements in reverse order. However, it is not easyto implement the alternative algorithm without a temporary storage array and weprefer the above in-place algorithm.Finally, we mention that there exist fast cosine transforms for small N that donot rely on an auxiliary function or use an FFT routine.

Instead, they carry out thetransform directly, often coded in hardware for fixed N of small dimension [1]..