c12-3 (779559), страница 3
Текст из файла (страница 3)
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).1yj = sin(jπ/N )(fj + fN−j ) + (fj − fN−j )251712.3 FFT of Real Functions, Sine and Cosine Transforms#include <math.h>The sine transform, curiously, is its own inverse. If you apply it twice, you get theoriginal data, but multiplied by a factor of N/2.The other common boundary condition for differential equations is that thederivative of the function is zero at the boundary.
In this case the natural transformis 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 arrayto create an even array of double the length, and/or from whether the extended arraycontains 2N − 1, 2N , or some other number of points. In practice, only two of thenumerous 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 =N−1X1[f0 + (−1)k fN ] +fj cos(πjk/N )2j=1(12.3.17)It results from extending the given array to an even array about j = N , withf2N−j = fj ,j = 0, . .
. , N − 1(12.3.18)If you substitute this extended array into equation (12.3.9), and follow steps analogousto those leading up to equation (12.3.11), you will find that the Fourier transform isSample 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).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 arrayshould 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<=(n>>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.}}518Chapter 12.Fast Fourier Transformjust 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 is1(fj + fN−j ) − sin(jπ/N )(fj − fN−j )2j = 0, . . . , N − 1 (12.3.19)Instead of equation (12.3.15), realft now givesF2k = RkF2k+1 = F2k−1 + Ikk = 0, .
. . , (N/2 − 1)(12.3.20)The starting value for the recursion for odd k in this case isN−1X1F1 = (f0 − fN ) +fj cos(jπ/N )2j=1(12.3.21)This sum does not appear naturally among the Rk and Ik , and so we accumulate itduring the generation of the array yj .Once again this transform is its own inverse, and so the following routineworks for both directions of the transformation. Note that although this form ofthe cosine transform has N + 1 input and output values, it passes an array onlyof length N to realft.#include <math.h>#define PI 3.141592653589793void cosft1(float y[], int n)Calculates the cosine transform of a set y[1..n+1] of real-valued data points. The transformeddata replace the original data in array y. n must be a power of 2.
This program, withoutchanges, also calculates the inverse cosine transform, but in this case the output array shouldbe 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;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<=(n>>1);j++) {wr=(wtemp=wr)*wpr-wi*wpi+wr;wi=wi*wpr+wtemp*wpi+wi;y1=0.5*(y[j]+y[n2-j]);y2=(y[j]-y[n2-j]);y[j]=y1-wi*y2;y[n2-j]=y1+wi*y2;sum += wr*y2;}Initialize the recurrence.j=n/2+1 unnecessary since y[n/2+1] unchanged.Carry out the recurrence.Calculate the auxiliary function.The values for j and N − j are related.Carry along this sum for later use in unfolding the transform.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).yj =12.3 FFT of Real Functions, Sine and Cosine Transforms519Calculate the transform of the auxiliary function.sum is the value of F1 in equation (12.3.21).realft(y,n,1);y[n+1]=y[2];y[2]=sum;for (j=4;j<=n;j+=2) {sum += y[j];y[j]=sum;}Equation (12.3.20).The second important form of the cosine transform is defined byFk =N−1Xfj cosj=0πk(j + 12 )N(12.3.22)with inverseN−1πk(j + 12 )2 X0fj =Fk cosNN(12.3.23)k=0Here the prime on the summation symbol means that the term for k = 0 has acoefficient of 12 in front.
This form arises by extending the given data, defined forj = 0, . . . , N − 1, to j = N, . . . , 2N − 1 in such a way that it is even about the pointN − 12 and periodic. (It is therefore also even about j = − 12 .) The form (12.3.23)is related to Gauss-Chebyshev quadrature (see equation 4.5.19), to Chebyshevapproximation (§5.8, equation 5.8.7), and Clenshaw-Curtis quadrature (§5.9).This form of the cosine transform is useful when solving differential equationson “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 =π(j + 12 )1(fj + fN−j−1 ) − sin(fj − fN−j−1 )2Nj = 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πkπkRk − sinIkNNπkπkRk + cosIk + F2k+1F2k−1 = sinNNF2k = cos(12.3.25)(12.3.26)Note that equation (12.3.26) givesFN−1 =1RN/22(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 stepsto find the inverse.
Here is the routine: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).}520Chapter 12.Fast Fourier Transform#include <math.h>#define PI 3.141592653589793theta=0.5*PI/n;wr1=cos(theta);wi1=sin(theta);wpr = -2.0*wi1*wi1;wpi=sin(2.0*theta);if (isign == 1) {for (i=1;i<=n/2;i++) {y1=0.5*(y[i]+y[n-i+1]);y2=wi1*(y[i]-y[n-i+1]);y[i]=y1+y2;y[n-i+1]=y1-y2;wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1;wi1=wi1*wpr+wtemp*wpi+wi1;}realft(y,n,1);for (i=3;i<=n;i+=2) {wr=(wtemp=wr)*wpr-wi*wpi+wr;wi=wi*wpr+wtemp*wpi+wi;y1=y[i]*wr-y[i+1]*wi;y2=y[i+1]*wr+y[i]*wi;y[i]=y1;y[i+1]=y2;}sum=0.5*y[2];for (i=n;i>=2;i-=2) {sum1=sum;sum += y[i];y[i]=sum1;}} else if (isign == -1) {ytemp=y[n];for (i=n;i>=4;i-=2) y[i]=y[i-2]-y[i];y[2]=2.0*ytemp;for (i=3;i<=n;i+=2) {wr=(wtemp=wr)*wpr-wi*wpi+wr;wi=wi*wpr+wtemp*wpi+wi;y1=y[i]*wr+y[i+1]*wi;y2=y[i+1]*wr-y[i]*wi;y[i]=y1;y[i+1]=y2;}realft(y,n,-1);for (i=1;i<=n/2;i++) {y1=y[i]+y[n-i+1];y2=(0.5/wi1)*(y[i]-y[n-i+1]);y[i]=0.5*(y1+y2);y[n-i+1]=0.5*(y1-y2);wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1;wi1=wi1*wpr+wtemp*wpi+wi1;Initialize the recurrences.Forward transform.Calculate the auxiliary function.Carry out the recurrence.Transform the auxiliary function.Even terms.Initialize recurrence for odd termswith 12 RN/2 .Carry out recurrence for odd terms.Inverse transform.Form difference of odd terms.Calculate Rk and Ik .Invert auxiliary array.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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
