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 <math.h> 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<=(n>>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 <math.h> 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<=(n>>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)