12.2 Fast Fourier Transform(FFT) 507 length of the real array (data[1..2*nn])is 2 times nn.with each complex value occupying two consecutive locations.In other words,data[]is the real part of fo,data[2]is the imaginary part of fo,and so on up to data[2*nn-1],which is the real part of fN-1,and data[2*nn],which is the imaginary part of fN-1. The FFT routine gives back the Fn's packed in exactly the same fashion,as nn complex numbers. The real and imaginary parts ofthe zero frequency component Foare in data[1] and data[2];the smallest nonzero positive frequency has real and imaginary parts in data[3]and data [4]:the smallest (in magnitude)nonzero negative frequency has real and imaginary parts in data[2*nn-1]and data[2*nn].Positive frequencies 8 increasing in magnitude are stored in the real-imaginary pairs data [5],data [6] up to data[nn-1],data [nn].Negative frequencies of increasing magnitude are stored in data[2*nn-3],data[2*nn-2]down to data [nn+3],data [nn+4] Finally,the pair data[nn+1],data [nn+2]contain the real and imaginary parts of the one aliased point that contains the most positive and the most negative frequency. You should try to develop a familiarity with this storage arrangement of complex from NUMERICAL RECIPES I 18881892 spectra,also shown in Figure 12.2.2,since it is the practical standard. This is a good place to remind you that you can also use a routine like four1 without modification even if your input data array is zero-offset,that is has the range data[o..2*nn-1].In this case,simply decrement the pointer to data by one when Press. THE four1 is invoked,e.g.,four1(data-1,1024,1);.The real part of fo will now be ART returned in data[o],the imaginary part in data[1],and so on.See 81.2. 9 #include <math.h> Program #define SWAP(a,b)tempr=(a);(a)=(b);(b)=tempr void fouri(float data[],unsigned long nn,int isign) Replaces data[1..2*nn]by its discrete Fourier transform,if isign is input as 1;or replaces data[1..2*nn]by nn times its inverse discrete Fourier transform,if isign is input as -1. data is a complex array of length nn or,equivalently,a real array of length 2*nn.nn MUST be an integer power of 2 (this is not checked for!). OF SCIENTIFIC COMPUTING(ISBN unsigned long n,mmax,m,j,istep,i; 198918920 double wtemp,wr,wpr,wpi,wi,theta; Double precision for the trigonomet- float tempr,tempi; ric recurrences. Numerical Recipes 10-621 n=nn<<1; j=1; for(1=1;1<n;i+=2)[ This is the bit-reversal section of the 43108 1f(j>1)[ routine. SWAP(data[i],data[i]); Exchange the two complex numbers. SWAP(data[j+1],data[i+1]); (outside 膜 Software. menn: wh11e(m>=2&j>m){ j-=m; m>>=1; j+=m; Here begins the Danielson-Lanczos section of the routine. mmax=2: while (n mmax){ Outer loop executed log2 nn times. istep=mmax <1; theta=is1gn*(6.28318530717959/mmax); Initialize the trigonometric recurrence. wtemp=sin(0.5*theta); wpr =-2.0*wtemp*wtemp;12.2 Fast Fourier Transform (FFT) 507 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). length of the real array (data[1..2*nn]) is 2 times nn, with each complex value occupying two consecutive locations. In other words, data[1] is the real part of f0, data[2] is the imaginary part of f0, and so on up to data[2*nn-1], which is the real part of fN−1, and data[2*nn], which is the imaginary part of fN−1. The FFT routine gives back the Fn’s packed in exactly the same fashion, as nn complex numbers. The real and imaginary parts of the zero frequency component F 0 are in data[1] and data[2]; the smallest nonzero positive frequency has real and imaginary parts in data[3] and data[4]; the smallest (in magnitude) nonzero negative frequency has real and imaginary parts in data[2*nn-1] and data[2*nn]. Positive frequencies increasing in magnitude are stored in the real-imaginary pairs data[5], data[6] up to data[nn-1], data[nn]. Negative frequencies of increasing magnitude are stored in data[2*nn-3], data[2*nn-2] down to data[nn+3], data[nn+4]. Finally, the pair data[nn+1], data[nn+2] contain the real and imaginary parts of the one aliased point that contains the most positive and the most negative frequency. You should try to develop a familiarity with this storage arrangement of complex spectra, also shown in Figure 12.2.2, since it is the practical standard. This is a good place to remind you that you can also use a routine like four1 without modification even if your input data array is zero-offset, that is has the range data[0..2*nn-1]. In this case, simply decrement the pointer to data by one when four1 is invoked, e.g., four1(data-1,1024,1);. The real part of f0 will now be returned in data[0], the imaginary part in data[1], and so on. See §1.2. #include <math.h> #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void four1(float data[], unsigned long nn, int isign) Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as 1; or replaces data[1..2*nn] by nn times its inverse discrete Fourier transform, if isign is input as −1. data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn MUST be an integer power of 2 (this is not checked for!). { unsigned long n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; Double precision for the trigonomet￾float tempr,tempi; ric recurrences. n=nn << 1; j=1; for (i=1;i<n;i+=2) { This is the bit-reversal section of the if (j > i) { routine. SWAP(data[j],data[i]); Exchange the two complex numbers. SWAP(data[j+1],data[i+1]); } m=nn; while (m >= 2 && j > m) { j -= m; m >>= 1; } j += m; } Here begins the Danielson-Lanczos section of the routine. mmax=2; while (n > mmax) { Outer loop executed log2 nn times. istep=mmax << 1; theta=isign*(6.28318530717959/mmax); Initialize the trigonometric recurrence. wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp;