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

In order to actuallyget the entire Fn in the original array space, it is convenient to put FN/2into the imaginary part of F0 .• Despite its complicated form, the process above is invertible. First peelFN/2 out of F0 . Then constructn = 0, . . . , N/2 − 1(1)(12.3.6)(2)and use four1 to find the inverse transform of Hn = Fn + iFn .Surprisingly, the actual algebraic steps are virtually identical to those ofthe 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 (whichis 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 elementsdata[1] and data[2], respectively. n must be a power of 2. This routine also calculates theinverse transform of a complex data array if it is the transform of real data.

(Result in this casemust 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 formdata[i2]=h1i+wr*h2i+wi*h2r;the true transform of the origidata[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.wi=wi*wpr+wtemp*wpi+wi;}if (isign == 1) {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).1*(Fn + FN/2−n)21*Fno = e−2πin/N (Fn − FN/2−n)2Fne =514Chapter 12.Fast Fourier Transformdata[1] = (h1r=data[1])+data[2];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);}Squeeze the first and last data together to get them all within theoriginal array.This is the inverse transform for thecase isign=-1.Fast Sine and Cosine TransformsAmong their other uses, the Fourier transforms of functions can be used to solvedifferential equations (see §19.4). The most common boundary conditions for thesolutions are 1) they have the value zero at the boundaries, or 2) their derivativesare zero at the boundaries.

In the first instance, the natural transform to use is thesine transform, given byFk =N−1Xfj sin(πjk/N )sine transform(12.3.7)j=1where 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 Fouriertransform. However, the argument of the sine differs by a factor of two from thevalue that would make this so. The sine transform uses sines only as a complete setof functions in the interval from 0 to 2π, and, as we shall see, the cosine transformuses cosines only.

By contrast, the normal FFT uses both sines and cosines, but onlyhalf as many of each. (See Figure 12.3.1.)The expression (12.3.7) can be “force-fit” into a form that allows its calculationvia the FFT. The idea is to extend the given function rightward past its last tabulatedvalue. We extend the data to twice their length in such a way as to make them anodd function about j = N , with fN = 0,f2N−j ≡ −fjj = 0, . .

. , N − 1(12.3.8)Consider the FFT of this extended function:Fk =2N−1Xfj e2πijk/(2N)(12.3.9)j=0The half of this sum from j = N to j = 2N − 1 can be rewritten with thesubstitution j 0 = 2N − j2N−1Xj=Nfj e2πijk/(2N) =NXf2N−j 0 e2πi(2N−j0)k/(2N)j 0 =1=−N−1Xj 0=0(12.3.10)fj 0 e−2πij0k/(2N)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).}12.3 FFT of Real Functions, Sine and Cosine Transforms1+124(a)033+121045−11+123(c)04−1502π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 realand imaginary parts of the basis functions are both shown.) While some basis functions occur in morethan 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 intervalshown; however, the sine or cosine transform best expands functions matching the boundary conditionsof the respective basis functions, namely zero function values for sine, zero derivatives for cosine.so thatFk =N−1Xhifj e2πijk/(2N) − e−2πijk/(2N)j=0= 2iN−1X(12.3.11)fj sin(πjk/N )j=0Thus, up to a factor 2i we get the sine transform from the FFT of the extendedfunction.This method introduces a factor of two inefficiency into the computation byextending the data.

This inefficiency shows up in the FFT output, which haszeros for the real part of every element of the transform. For a one-dimensionalproblem, the factor of two may be bearable, especially in view of the simplicityof the method. When we work with partial differential equations in two or threedimensions, though, the factor becomes four or eight, so efforts to eliminate theinefficiency are well rewarded.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).5−1(b)515516Chapter 12.Fast Fourier TransformFrom the original real data array fj we will construct an auxiliary array yj andapply to it the routine realft. The output will then be used to construct the desiredtransform. For the sine transform of data fj , j = 1, . . . , N − 1, the auxiliary array isy0 = 0j = 1, .

. . , N − 1(12.3.12)This array is of the same dimension as the original. Notice that the first term issymmetric about j = N/2 and the second is antisymmetric. Consequently, whenrealft is applied to yj , the result has real parts Rk and imaginary parts Ik given byRk =N−1Xyj cos(2πjk/N )j=0=N−1X(fj + fN−j ) sin(jπ/N ) cos(2πjk/N )j=1=N−1X2fj sin(jπ/N ) cos(2πjk/N )j=0=N−1Xj=0fj(2k − 1)jπ(2k + 1)jπ− sinsinNN= F2k+1 − F2k−1Ik =N−1X(12.3.13)yj sin(2πjk/N )j=0==N−1X1(fj − fN−j ) sin(2πjk/N )2j=1N−1Xfj sin(2πjk/N )j=0(12.3.14)= F2kTherefore Fk can be determined as follows:F2k = IkF2k+1 = F2k−1 + Rkk = 0, .

. . , (N/2 − 1)(12.3.15)The even terms of Fk are thus determined very directly. The odd terms requirea recursion, the starting point of which follows from setting k = 0 in equation(12.3.15) and using F1 = −F−1 :F1 =The implementing program is1R02(12.3.16)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.