c12-3 (Numerical Recipes in C)

510Chapter 12.Fast Fourier TransformCITED REFERENCES AND FURTHER READING:Nussbaumer, H.J. 1982, Fast Fourier Transform and Convolution Algorithms (New York: SpringerVerlag).Elliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (NewYork: Academic Press).Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall). [1]Bloomfield, 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 (October).12.3 FFT of Real Functions, Sine and CosineTransformsIt happens frequently that the data whose FFT is desired consist of real-valuedsamples fj , j = 0 . . . N − 1. To use four1, we put these into a complex arraywith all imaginary parts set to zero. The resulting transform Fn , n = 0 . . . N − 1satisfies FN−n * = Fn .

Since this complex-valued array has real values for F0and FN/2 , and (N/2) − 1 other independent values F1 . . . FN/2−1 , it has the same2(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 inefficient, both in executiontime and in storage required. You would think that there is a better way.There are two better ways. The first is “mass production”: Pack two separatereal functions into the input array in such a way that their individual transforms canbe 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 twoof an item when you only need one.

However, remember that for correlations andconvolutions the Fourier transforms of two functions are involved, and this is ahandy way to do them both at once. The second method is to pack the real inputarray cleverly, without extra zeros, into a complex array of half its length. One thenperforms a complex FFT on this shorter length; the trick is then to get the requiredanswer out of the result. This is done in the program realft below.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).integer arithmetic modulo some large prime N +1, and the N th root of 1 by themodulo arithmetic equivalent. Strictly speaking, these are not Fourier transformsat all, but the properties are quite similar and computational speed can be farsuperior. On the other hand, their use is somewhat restricted to quantities likecorrelations and convolutions since the transform itself is not easily interpretableas a “frequency” spectrum.12.3 FFT of Real Functions, Sine and Cosine Transforms511Transform of Two Real Functions SimultaneouslyFirst we show how to exploit the symmetry of the transform Fn to handletwo real functions at once: Since the input data fj are real, the components of thediscrete Fourier transform satisfy(12.3.1)where the asterisk denotes complex conjugation.

By the same token, the discreteFourier 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 oflength N simultaneously by packing the two data arrays as the real and imaginaryparts, respectively, of the complex input array of four1. Then the resulting transformarray 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 andreturns two complex output arrays, fft1[1..2n] and fft2[1..2n], each of complex lengthn (i.e., real length 2*n), which contain the discrete Fourier transforms of the respective dataarrays.

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<=n;j++,jj+=2) {fft1[jj-1]=data1[j];fft1[jj]=data2[j];}four1(fft1,n,1);fft2[1]=fft1[2];fft1[2]=fft2[2]=0.0;for (j=3;j<=n+1;j+=2) {rep=0.5*(fft1[j]+fft1[nn2-j]);rem=0.5*(fft1[j]-fft1[nn2-j]);aip=0.5*(fft1[j+1]+fft1[nn3-j]);aim=0.5*(fft1[j+1]-fft1[nn3-j]);fft1[j]=rep;fft1[j+1]=aim;fft1[nn2-j]=rep;fft1[nn3-j] = -aim;fft2[j]=aip;fft2[j+1] = -rem;fft2[nn2-j]=aip;fft2[nn3-j]=rem;}}Pack the two real arrays into one complex array.Transform the complex array.Use symmetries to separate the two transforms.Ship them out in two complex arrays.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).FN−n = (Fn )*512Chapter 12.Fast Fourier TransformFFT of Single Real FunctionTo implement the second method, which allows us to perform the FFT ofa single real function without redundancy, we split the data set in half, therebyforming two real arrays of half the size.

We can apply the program above to thesetwo, but of course the result will not be the transform of the original data. It willbe a schizophrenic combination of two transforms, each of which has half of theinformation 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 asone data set, and the odd-numbered fj as the other. The beauty of this is thatwe can take the original real array and treat it as a complex array hj of half thelength.

The first data set is the real part of this array, and the second is theimaginary part, as prescribed for twofft. No repacking is required. In other wordshj = f2j + if2j+1 , j = 0, . . . , N/2 − 1. We submit this to four1, and it will giveback a complex array Hn = Fne + iFno , n = 0, . . . , N/2 − 1 withXN/2−1Fne =f2k e2πikn/(N/2)k=0X(12.3.3)N/2−1Fno =f2k+1 e2πikn/(N/2)k=0The discussion of program twofft tells you how to separate the two transformsFne and Fno out of Hn . How do you work them into the transform Fn of the originaldata set fj ? Simply glance back at equation (12.2.3):Fn = Fne + e2πin/N Fnon = 0, .

. . , N − 1(12.3.4)Expressed directly in terms of the transform Hn of our real (masquerading ascomplex) data set, the result isFn =1i(Hn + HN/2−n *) − (Hn − HN/2−n *)e2πin/N22n = 0, . . . , N − 1(12.3.5)A few remarks:• Since FN−n * = Fn there is no point in saving the entire spectrum. Thepositive frequency half is sufficient and can be stored in the same array asthe original data.

The operation can, in fact, be done in place.• Even so, we need values Hn , n = 0, . . . , N/2 whereas four1 gives onlythe values n = 0, . . . , N/2 − 1. Symmetry to the rescue, HN/2 = H0 .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).What about the reverse process? Suppose you have two complex transformarrays, each of which has the symmetry (12.3.1), so that you know that the inversesof both transforms are real functions.

Can you invert both in a single FFT? This iseven easier than the other direction. Use the fact that the FFT is linear and formthe sum of the first transform plus i times the second. Invert using four1 withisign = −1. The real and imaginary parts of the resulting complex array are thetwo desired real functions.51312.3 FFT of Real Functions, Sine and Cosine Transforms• The values F0 and FN/2 are real and independent.