510 Chapter 12.Fast Fourier Transform integer arithmetic modulo some large prime N+1,and the Nth 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 asa“frequency”spectrum. 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] 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.IA.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- 9 tober). Press. 9 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 fi,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 satisfies FN*FSince this complex-valued array has real values for Fo 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, 兴彩 Numerical Recipes 10.621 the use of the full complex FFT algorithm for real data is inefficient.both in execution 43106 time 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 separate real functions into the input array in such a way that their individual transforms can (outside be separated from the result.This is implemented in the program twofft below. North 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
510 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 machinereadable 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). integer arithmetic modulo some large prime N+1, and the Nth 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. CITED 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 (New York: 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 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 F n, n = 0 ...N − 1 satisfies 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 inefficient, both in execution time 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 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 f;are real,the components of the discrete Fourier transform satisfy 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 gi'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. Americ computer, void twofft(f1 oat data1☐，f1 oat data2[☐，f1 oat fft1[0,f1 oat fft2[], unsigned long n) Given two real input arrays data1[1..n]and data2[1..n],this routine calls four1 and 9 Programs 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 SCIENTIFIC 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,aimi nn3=1+(nn2=2+n+n); 1920 COMPUTING(ISBN for（j=1,jj=2;j<=n;j+,jj+=2)( Pack the two real arrays into one com- fft1[jj-1]=datal[j]; plex array. ffti[jj]=data2[j]; Numerical Recipes 021 four1(ffti,n,1); Transform the complex array. fft2[1]=fft1[2]; 43108 fft1[2]=fft2[2]=0.0; for(j=3;j<=n+1;j+=2){ rep=0.5*(fft1[j]+fft1[nn2-j]); Use symmetries to separate the two trans- (outside rem=0.5*(fft1[j]-fft1[nn2-j]); forms. Software. aip=0.5*(fft1[i+1]+fft1[nn3-i]); aim=0.5*(fft1[j+1]-fft1[nn3-j]); ffti[j]-rep; Ship them out in two complex arrays. fft1[i+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;
12.3 FFT of Real Functions, Sine and Cosine Transforms 511 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 machinereadable 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). Transform of Two Real Functions Simultaneously First we show how to exploit the symmetry of the transform F n to handle two real functions at once: Since the input data fj are real, the components of the discrete Fourier transform satisfy 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 g j ’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<=n;j++,jj+=2) { Pack the two real arrays into one comfft1[jj-1]=data1[j]; plex array. fft1[jj]=data2[j]; } four1(fft1,n,1); Transform the complex array. 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]); Use symmetries to separate the two transrem=0.5*(fft1[j]-fft1[nn2-j]); forms. aip=0.5*(fft1[j+1]+fft1[nn3-j]); aim=0.5*(fft1[j+1]-fft1[nn3-j]); fft1[j]=rep; Ship them out in two complex arrays. 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; } }
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 first transform plus i times the second.Invert using fouri with isign=-1.The real and imaginary parts of the resulting complex array are the 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 虽2 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 fi as RECIPESI one data set,and the odd-numbered fi as the other.The beauty of this is that we can take the original real array and treat it as a complex array h;of half the 9 length.The first 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 arrayHF+iF,n=0....,N/2-1 with g3 9 N/2-1 F= f2ke2rikn/(N/2) OF SCIENTIFIC( (12.3.3) 、，文 61 N/2-1 Fo= f2k+1e2mikn/(N/2) The discussion of program twofft tells you how to separate the two transforms F and F out of Hn.How do you work them into the transform Fn of the original 10621 data set f;?Simply glance back at equation (12.2.3): Numerica Fn =Fe+e2rin/N Fo n=0,...,N-1 43106 (12.3.4) Expressed directly in terms of the transform Hn of our real(masquerading as (outside 腿 complex)data set,the result is North Fn(H+Hx/-n)-(H-HN/a-n")e2sin/N 1 2 n=0,N-1 (12.3.5) A few remarks: .Since FN-n*=Fn there is no point in saving the entire spectrum.The positive frequency half is sufficient 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 Ho
512 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 machinereadable 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). 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 first transform plus i times the second. Invert using four1 with isign = −1. The real and imaginary parts of the resulting complex array are the 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 f j 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 first 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 = Fe n + iFo n, n = 0, . . ., N/2 − 1 with Fe n = N/ 2−1 k=0 f2k e2πikn/(N/2) Fo n = N/ 2−1 k=0 f2k+1 e2πikn/(N/2) (12.3.3) The discussion of program twofft tells you how to separate the two transforms Fe n and Fo n out of Hn. How do you work them into the transform F n of the original data set fj ? Simply glance back at equation (12.2.3): Fn = Fe n + e2πin/N Fo n 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 Fn = 1 2 (Hn + HN/2−n*) − i 2 (Hn − HN/2−n*)e2πin/N n = 0,...,N − 1 (12.3.5) A few remarks: • Since FN−n* = Fn there is no point in saving the entire spectrum. The positive frequency half is sufficient 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 Fo 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 Fo. Despite its complicated form,the process above is invertible.First peel FN/2 out of Fo.Then construct 1 F%=2(B+Fe-n） n=0,,N/2-1(12.3.6) (F-F) and use four1 to find the inverse transform of H=) 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. Press. THE The real-valued first 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. Programs void four1(float data[],unsigned long nn,int isign); ms1ged1ong1,11,12,13,14,np3; SCIENTIFIC float c1=0.5,c2,hir,hii,h2r,h2i; double wr,wi,wpr,wpi,wtemp,theta; Double precision for the trigonomet- ric recurrences. theta=3.141592653589793/(doub1e)(n>>1); Initialize the recurrence. 1f(1s1gm=1){ c2=-0.5; four1(data,n>>1,1); The forward transform is here 1920 COMPUTING (ISBN else c2=0.5; Otherwise set up for an inverse trans- theta =-theta; form. 10-621 wtemp=sin(0.5*theta); wpr =-2.0*vtemp*wtemp; Numerical Recipes -43106 wpi=sin(theta): r=1.0+pr; wi-wpi; (outside np3=n+3: for(i=2;1>2);i++) Software. Case i=1 done separately below. 14=1+(13=np3-(12=1+(11=1+1-1))); hir-c1*(data[i1]+data[i3]); The two separate transforms are sep- hii=c1*(data[i2]-data[i4]); arated out of data h2r =-c2*(data[i2]+data[i4]); h2i=c2*(data[i1]-data[i3]); data[i1]=hir+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 wiewi*wpr+wtemp*wpi+wi; 1f(is1gn1){
12.3 FFT of Real Functions, Sine and Cosine Transforms 513 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 machinereadable 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). • 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 Fe n = 1 2 (Fn + F* N/2−n) Fo n = 1 2 e−2πin/N (Fn − F* N/2−n) n = 0, . . ., N/2 − 1 (12.3.6) and use four1 to find the inverse transform of Hn = F(1) n + iF(2) n . 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 first 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 trigonometric 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 transtheta = -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 seph1i=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 the true transform of the original real data. data[i2]=h1i+wr*h2i+wi*h2r; data[i3]=h1r-wr*h2r+wi*h2i; 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]=(hir=data[1])+data[2]; Squeeze the first and last data to- data[2]hir-data[2]; gether to get them all within the else original array. data[1]=c1*((hir=data[1])+data[2]); data[2]=c1*(hir-data[2]); fouri(data,n>>1,-1); This is the inverse transform for the case isign=-1. Fast Sine and Cosine Transforms Among their other uses,the Fourier transforms of functions can be used to solve 83g differential equations (see 819.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 first instance,the natural transform to use is the sine transform,given by N-1 sine transform (12.3.7) 3 2 1 where fj,j=0,...,N-1 is the data array,and fo =0. Press. At first 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 9 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 SCIENTIFIC( 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-fit"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三-f方j=0,..,N-1 (12.3.8) Numerical 10521 4310 Consider the FFT of this extended function: Recipes 2W-1 (outside F压=∑je2eW (12.3.9) j=0 North The half of this sum from j=N to j=2N-1 can be rewritten with the substitution j'=2N-j 2N-1 N 五方e2n2M=hN-e2x2w-2M j'=1 (12.3.10) N-1 =-∑fe-2ay2M '=0
514 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 machinereadable 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). data[1] = (h1r=data[1])+data[2]; Squeeze the first and last data together to get them all within the original array. data[2] = h1r-data[2]; } else { 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. } 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 first instance, the natural transform to use is the sine transform, given by Fk = N −1 j=1 fj sin(πjk/N) sine transform (12.3.7) where fj , j = 0,...,N − 1 is the data array, and f0 ≡ 0. At first 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-fit” 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: Fk = 2 N−1 j=0 fj e2πijk/(2N) (12.3.9) The half of this sum from j = N to j = 2N − 1 can be rewritten with the substitution j = 2N − j 2 N−1 j=N fj e2πijk/(2N) = N j=1 f2N−je2πi(2N−j )k/(2N) = − N −1 j=0 fje−2πij k/(2N) (12.3.10)
12.3 FFT of Real Functions,Sine and Cosine Transforms 515 a 0 (b) 0 (including this one) internet 3 tusers to make one paper 1988-1992 by Cambridge University Press. from NUMERICAL RECIPES IN C: THE (c) 0 ART -1 Programs 0 2π ictly proh Figure 12.3.1.Basis functions used by the Fourier transform(a),sine transform(b),and cosine transform (c),are plotted.The first five basis functions are shown in each case.(For the Fourier transform,the rea and imaginary parts of the basis functions are both shown.)While some basis functions occur in more to dir 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 SCIENTIFIC COMPUTING (ISBN of the respective basis functions,namely zero function values for sine,zero derivatives for cosine. A是 1988-19920 so that 10-521 无，[e2k/2N-e2m/2N】 Numerical Recipes 43108 F=】 =0 (12.3.11)） N-1 (outside =2i>jsin(πjk/N) Software. j=0 Thus,up to a factor 2i we get the sine transform from the FFT ofthe extended function. This method introduces a factor of two inefficiency into the computation by extending the data.This inefficiency 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 inefficiency are well rewarded
12.3 FFT of Real Functions, Sine and Cosine Transforms 515 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 machinereadable 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). (a) +1 0 −1 +1 0 −1 +1 0 −1 (b) (c) 0 2π 5 4 2 1 3 3 2 1 4 5 1 2 3 4 5 Figure 12.3.1. Basis functions used by the Fourier transform (a), sine transform (b), and cosine transform (c), are plotted. The first five 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 Fk = N −1 j=0 fj e2πijk/(2N) − e−2πijk/(2N) = 2i N −1 j=0 fj sin(πjk/N) (12.3.11) 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 inefficiency into the computation by extending the data. This inefficiency 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 inefficiency are well rewarded.
516 Chapter 12.Fast Fourier Transform From the original real data array fi we will construct an auxiliary array y;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 0=0 (12.3.12) 斯=sin(jπ/Nj+fN-)+5(j-fw-)j=1,,N-1 This array is of the same dimension as the original.Notice that the first term is symmetric about j=N/2 and the second is antisymmetric.Consequently,when realft is applied to yi,the result has real parts R and imaginary parts Ik given by N-1 Rk= ∑cos(2rjk/N 1800 j=0 872 N- ∑G)+fN-)sim(Gπ/)cos(2πjk/N) to any j=1 N- ∑25si(Gm/Ncos2mjk/N = server computer, (North America e University Press. THE j= ART 2k+1i还-sm2k-1j西 Programs send =F2k+1-F2k-1 (12.3.13) N- Copyright (C) Ik= ∑斯sin(2iN to dir j= ∑5-fN-5sin2mjkN rectcustser 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 0-521- N-1 f方si(2πjk/N) j=0 -43108-5 =F2k (12.3.14) (outside Therefore F can be determined as follows: North Software. F2k Ik F2k+1=F2k-1+Rk k=0,.,(N/2-1) (12.3.15) The even terms of F&are thus determined very directly.The odd terms require visit website f machine a recursion,the starting point of which follows from settingk=0 in equation (12.3.15)and using F1=-F-1: F1=Ro (12.3.16) The implementing program is
516 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 machinereadable 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). 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 yj = sin(jπ/N)(fj + fN−j ) + 1 2 (fj − fN−j) j = 1,...,N − 1 (12.3.12) This array is of the same dimension as the original. Notice that the first 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 Rk = N −1 j=0 yj cos(2πjk/N) = N −1 j=1 (fj + fN−j ) sin(jπ/N) cos(2πjk/N) = N −1 j=0 2fj sin(jπ/N) cos(2πjk/N) = N −1 j=0 fj sin (2k + 1)jπ N − sin (2k − 1)jπ N = F2k+1 − F2k−1 (12.3.13) Ik = N −1 j=0 yj sin(2πjk/N) = N −1 j=1 (fj − fN−j ) 1 2 sin(2πjk/N) = N −1 j=0 fj sin(2πjk/N) = 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: F1 = 1 2 R0 (12.3.16) 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. 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,wtempi Double precision in the trigono- metric recurrences. thetas=3.14159265358979/(doub1e)n; Initialize the recurrence wtemp=sin(0.5*theta); wpr =-2.0*wtemp*wtemp; wpi=sin(theta); y[1]=0.0: for(j=2;j1)+1;j+)( wr=(wtemp=wr)*wpr-wi*wpi+wr; Calculate the sine for the auxiliary array. wiswi*wpr+wtemp*wpi+wi; The cosine is needed to continue the recurrence. y1=wi*(y[j]+y[n2-j]); Construct the auxiliary array. RECIPES 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. w》 Press. y[1]*=0.5; Initialize the sum used for odd terms below. sum=y[2]=0.0; for(j=1;:j<=n-1;j+=2)[ sum +y[j]; y[j]y[j+1]; Even terms determined directly. y[j+1]=sum; Odd terms determined by this running sum. IENTIFIC to dir The sine transform,curiously,is its own inverse.If you apply it twice,you get the original data,but multiplied by a factor of N/2 1920 The other common boundary condition for differential equations is that the derivative of the function is zero at the boundary.In this case the natural transform is the cosine transform.There are several possible ways of defining the transform. Numerica 10621 Each can be thought of as resulting from a different way of extending a given array 43106 to create an even array of double the length,and/or from whether the extended array Recipes contains 2N-1,2N,or some other number of points.In practice,only two of the (outside numerous possibilities are useful so we will restrict ourselves to just these two. The first form of the cosine transform uses N+1 data points: North Software. N-1 F=2L+(-1)*fN+∑f方cos(xjk/N) (12.3.17) j=1 It results from extending the given array to an even array about j=N,with f2N-j=f方，j=0,.,N-1 (12.3.18) Ifyou substitute this extended array into equation(12.3.9),and follow steps analogous to those leading up to equation(12.3.11),you will find that the Fourier transform is
12.3 FFT of Real Functions, Sine and Cosine Transforms 517 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 machinereadable 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). #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. { 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 trigonometric 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<=n-1;j+=2) { sum += y[j]; y[j]=y[j+1]; Even terms determined directly. y[j+1]=sum; Odd terms determined by this running sum. } } The sine transform, curiously, is its own inverse. If you apply it twice, you get the original data, but multiplied by a factor of N/2. The other common boundary condition for differential equations is that the derivative of the function is zero at the boundary. In this case the natural transform is the cosine transform. There are several possible ways of defining the transform. Each can be thought of as resulting from a different way of extending a given array to create an even array of double the length, and/or from whether the extended array contains 2N − 1, 2N, or some other number of points. In practice, only two of the numerous possibilities are useful so we will restrict ourselves to just these two. The first form of the cosine transform uses N + 1 data points: Fk = 1 2 [f0 + (−1)kfN ] + N −1 j=1 fj cos(πjk/N) (12.3.17) It results from extending the given array to an even array about j = N, with f2N−j = fj , j = 0,...,N − 1 (12.3.18) If you substitute this extended array into equation (12.3.9),and follow steps analogous to those leading up to equation (12.3.11), you will find that the Fourier transform is
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 S4.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 inefficiency. In this case the auxiliary function is 5=2G+fN-)-sim(Gπ/N)(5-fN-）) j=0,..,N-1(123.19) Instead of equation (12.3.15),realft now gives F2k=Rk F2k+1=F2k-1+Ikk=0,.,(N/2-1) (12.3.20) The starting value for the recursion for odd k in this case is N-1 =2(6-fw)+ fj cos(iπ/N) (12.3.21) RECIPES 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 会 Press. 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 Programs of length N to realft. SCIENTIFIC #include #def1nePI3.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 1920 COMPUTING (ISBN be multiplied by 2/n. void realft(float data[],unsigned long n,int isign); 10-621 int j,n2; float sum,y1,y2; double theta,wi=0.0,wpi,wpr,wr=1.0,wtemp; Numerical Recipes 43108 Double precision for the trigonometric recurrences. theta=PI/n; Initialize the recurrence. (outside 膜 wtemp=sin(0.5*theta); Software. wpr =-2.0*vtemp*wtemp; wpi=sin(theta); sum=0.5*(y[1]-y[n+1]); ying of y[1]=0.5*(y[1]+y[m+1]); n2=n+2; for(j=2;j1);j+) j=n/2+1 unnecessary since y [n/2+1]unchanged. wr=(wtemp=wr)*wpr-wi*wpi+wr; Carry out the recurrence wiswi*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-1*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
518 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 machinereadable 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). 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 inefficiency. In this case the auxiliary function is yj = 1 2 (fj + fN−j) − sin(jπ/N)(fj − fN−j) j = 0,...,N − 1 (12.3.19) 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 F1 = 1 2 (f0 − fN ) + N −1 j=1 fj cos(jπ/N) (12.3.21) 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 F in equation (12.3.21). for(j=4;j<=n;j+=2){ sum +y[j]; Equation (12.3.20). y[j]=sum; The second important form of the cosine transform is defined by N-1 Fk= cos πk(j+) (12.3.22) .com or j=0 with inverse N-1 2 Fk.cos k(行+) fj= N k=0 W (12.3.23) RECIPESI 、” 2 Here the prime on the summation symbol means that the term for k =0 has a coefficient of in front.This form arises by extending the given data,defined for j=0,...,N-1,to j N,...,2N-1 in such a way that it is even about the point Press. N-and periodic.(It is therefore also even aboutj=-.)The form (12.3.23) is related to Gauss-Chebyshev quadrature (see equation 4.5.19),to Chebyshev approximation(85.8,equation 5.8.7),and Clenshaw-Curtis quadrature(85.9). This form of the cosine transform is useful when solving differential equations SCIENTIFIC( on"staggered"grids,where the variables are centered midway between mesh points. It is also the standard form in the field of data compression and image processing. 61 The auxiliary function used in this case is similar to equation(12.3.19): 斯=+N-+血6+6-N-) j=0,..,N-1 9 (12.3.24) Recipes Numerical 10621 Carrying out the steps similar to those used to get from(12.3.12)to (12.3.15),we find 43106 A=A-m装 (12.3.25) F2k-1 sin- +cosFake nk (12.3.26) Note that equation (12.3.26)gives FN-1=RN/2 (12.3.27) Thus the even components are found directly from(12.3.25),while the odd com- ponents are found by recursing(12.3.26)down from k=N/2-1,using(12.3.27) to start. Since the transform is not self-inverting,we have to reverse the above steps to find the inverse.Here is the routine:
12.3 FFT of Real Functions, Sine and Cosine Transforms 519 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 machinereadable 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). realft(y,n,1); Calculate the transform of the auxiliary funcy[n+1]=y[2]; tion. y[2]=sum; sum is the value of F1 in equation (12.3.21). for (j=4;j<=n;j+=2) { sum += y[j]; Equation (12.3.20). y[j]=sum; } } The second important form of the cosine transform is defined by Fk = N −1 j=0 fj cos πk(j + 1 2 ) N (12.3.22) with inverse fj = 2 N N −1 k=0 Fk cos πk(j + 1 2 ) N (12.3.23) Here the prime on the summation symbol means that the term for k = 0 has a coefficient of 1 2 in front. This form arises by extending the given data, defined for j = 0,...,N − 1, to j = N,..., 2N − 1 in such a way that it is even about the point N − 1 2 and periodic. (It is therefore also even about j = − 1 2 .) The form (12.3.23) is related to Gauss-Chebyshev quadrature (see equation 4.5.19), to Chebyshev approximation (§5.8, equation 5.8.7), and Clenshaw-Curtis quadrature (§5.9). This form of the cosine transform is useful when solving differential equations on “staggered” grids, where the variables are centered midway between mesh points. It is also the standard form in the field of data compression and image processing. The auxiliary function used in this case is similar to equation (12.3.19): yj = 1 2 (fj + fN−j−1) + sin π(j + 1 2 ) N (fj − fN−j−1) j = 0,...,N − 1 (12.3.24) Carrying out the steps similar to those used to get from (12.3.12) to (12.3.15), we find F2k = cos πk N Rk − sin πk N Ik (12.3.25) F2k−1 = sin πk N Rk + cos πk N Ik + F2k+1 (12.3.26) Note that equation (12.3.26) gives FN−1 = 1 2 RN/2 (12.3.27) Thus the even components are found directly from (12.3.25), while the odd components are found by recursing (12.3.26) down from k = N/2 − 1, using (12.3.27) to start. Since the transform is not self-inverting, we have to reverse the above steps to find the inverse. Here is the routine: