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

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).float data[1..nn1][1..nn2][1..nn3]528Chapter 12.Fast Fourier Transform#include <math.h>void rlft3(float ***data, float **speq, unsigned long nn1, unsigned long nn2,unsigned long nn3, int isign)Given a three-dimensional real array data[1..nn1][1..nn2][1..nn3] (where nn1 = 1 forthe case of a logically two-dimensional array), this routine returns (for isign=1) the complexfast Fourier transform as two complex arrays: On output, data contains the zero and positivefrequency values of the third frequency component, while speq[1..nn1][1..2*nn2] containsthe Nyquist critical frequency values of the third frequency component.

First (and second)frequency components are stored for zero, positive, and negative frequencies, in standard wraparound order. See text for description of how complex values are arranged. For isign=-1, theinverse transform (times nn1*nn2*nn3/2 as a constant multiplicative factor) is performed,with output data (viewed as a real array) deriving from input data (viewed as complex) andspeq. For inverse transforms on data not generated first by a forward transform, make surethe complex input data array satisfies property (12.5.2). The dimensions nn1, nn2, nn3 mustalways be integer powers of 2.{void fourn(float data[], unsigned long nn[], int ndim, int isign);void nrerror(char error_text[]);unsigned long i1,i2,i3,j1,j2,j3,nn[4],ii3;double theta,wi,wpi,wpr,wr,wtemp;float c1,c2,h1r,h1i,h2r,h2i;if (1+&data[nn1][nn2][nn3]-&data[1][1][1] != nn1*nn2*nn3)nrerror("rlft3: problem with dimensions or contiguity of data array\n");c1=0.5;c2 = -0.5*isign;theta=isign*(6.28318530717959/nn3);wtemp=sin(0.5*theta);wpr = -2.0*wtemp*wtemp;wpi=sin(theta);nn[1]=nn1;nn[2]=nn2;nn[3]=nn3 >> 1;if (isign == 1) {Case of forward transform.fourn(&data[1][1][1]-1,nn,3,isign);Here is where most all of the comfor (i1=1;i1<=nn1;i1++)pute time is spent.for (i2=1,j2=0;i2<=nn2;i2++) {Extend data periodically into speq.speq[i1][++j2]=data[i1][i2][1];speq[i1][++j2]=data[i1][i2][2];}}for (i1=1;i1<=nn1;i1++) {j1=(i1 != 1 ? nn1-i1+2 : 1);Zero frequency is its own reflection, otherwise locate corresponding negative frequencyin wrap-around order.wr=1.0;Initialize trigonometric recurrence.wi=0.0;for (ii3=1,i3=1;i3<=(nn3>>2)+1;i3++,ii3+=2) {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).overwriting the input array of data. The three enclosing for loops (indices i2, i3,and i1, from inside to outside) could in fact be done in any order — their actions allcommute.

We chose the order shown because of the following considerations: (i) i3should not be the inner loop, because if it is, then the recurrence relations on wr andwi become burdensome. (ii) On virtual-memory machines, i1 should be the outerloop, because (with C order of array storage) this results in the array data, whichmight be very large, being accessed in block sequential order.Note that the work done in rlft3 is quite (logarithmically) small, compared tothe associated complex FFT, fourn. Since C does not have a convenient complextype, the operations are carried out explicitly below in terms of real and imaginaryparts.

The routine rlft3 is based on an earlier routine by G.B. Rybicki.12.5 Fourier Transforms of Real Data in Two and Three Dimensions529for (i2=1;i2<=nn2;i2++) {if (i3 == 1) {Equation (12.3.5).j2=(i2 != 1 ? ((nn2-i2)<<1)+3 : 1);h1r=c1*(data[i1][i2][1]+speq[j1][j2]);h1i=c1*(data[i1][i2][2]-speq[j1][j2+1]);h2i=c2*(data[i1][i2][1]-speq[j1][j2]);h2r= -c2*(data[i1][i2][2]+speq[j1][j2+1]);data[i1][i2][1]=h1r+h2r;data[i1][i2][2]=h1i+h2i;speq[j1][j2]=h1r-h2r;speq[j1][j2+1]=h2i-h1i;} else {j2=(i2 != 1 ? nn2-i2+2 : 1);j3=nn3+3-(i3<<1);h1r=c1*(data[i1][i2][ii3]+data[j1][j2][j3]);h1i=c1*(data[i1][i2][ii3+1]-data[j1][j2][j3+1]);h2i=c2*(data[i1][i2][ii3]-data[j1][j2][j3]);h2r= -c2*(data[i1][i2][ii3+1]+data[j1][j2][j3+1]);data[i1][i2][ii3]=h1r+wr*h2r-wi*h2i;data[i1][i2][ii3+1]=h1i+wr*h2i+wi*h2r;data[j1][j2][j3]=h1r-wr*h2r+wi*h2i;data[j1][j2][j3+1]= -h1i+wr*h2i+wi*h2r;}}wr=(wtemp=wr)*wpr-wi*wpi+wr;Do the recurrence.wi=wi*wpr+wtemp*wpi+wi;}}if (isign == -1)fourn(&data[1][1][1]-1,nn,3,isign);Case of reverse transform.}We now give some fragments from notional calling programs, to clarify theuse of rlft3 for two- and three-dimensional data.

Note again that the routine doesnot actually distinguish between two and three dimensions; two is treated like three,but with the first dimension having length 1. Since the first dimension is the outerloop, virtually no inefficiency is introduced.The first program fragment FFTs a two-dimensional data array, allows for someSample 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).Figure 12.5.2.

(a) A two-dimensional image with intensities either purely black or purely white. (b) Thesame image, after it has been low-pass filtered using rlft3. Regions with fine-scale features become gray.530Chapter 12.Fast Fourier Transform#include <stdlib.h>#include "nrutil.h"#define N2 256#define N3 256Note that the first component must be set to 1.int main(void)/* example1 */This fragment shows how one might filter a 256 by 256 digital image.{void rlft3(float ***data, float **speq, unsigned long nn1,unsigned long nn2, unsigned long nn3, int isign);float ***data, **speq;data=f3tensor(1,1,1,N2,1,N3);speq=matrix(1,1,1,2*N2);/*...*/Here the image would be loaded into data.rlft3(data,speq,1,N2,N3,1);/*...*/Here the arrays data and speq would be multiplied by arlft3(data,speq,1,N2,N3,-1);suitable filter function (of frequency)./*...*/Here the filtered image would be unloaded from data.free_matrix(speq,1,1,1,2*N2);free_f3tensor(data,1,1,1,N2,1,N3);return 0;}#define N1 32#define N2 64#define N3 16int main(void)/* example2 */This fragment shows how one might FFT a real three-dimensional array of size 32 by 64 by 16.{void rlft3(float ***data, float **speq, unsigned long nn1,unsigned long nn2, unsigned long nn3, int isign);int j;float ***data,**speq;data=f3tensor(1,N1,1,N2,1,N3);speq=matrix(1,N1,1,2*N2);/*...*/Here load data.rlft3(data,speq,N1,N2,N3,1);/*...*/Here unload data and speq.free_matrix(speq,1,N1,1,2*N2);free_f3tensor(data,1,N1,1,N2,1,N3);return 0;}#define N 32int main(void)/* example3 */This fragment shows how one might convolve two real, three-dimensional arrays of size 32 by32 by 32, replacing the first array by the result.{void rlft3(float ***data, float **speq, unsigned long nn1,unsigned long nn2, unsigned long nn3, int isign);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.