c12-5 (779561), страница 3
Текст из файла (страница 3)
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).processing on it, e.g., filtering, and then takes the inverse transform.
Figure 12.5.2shows an example of the use of this kind of code: A sharp image becomes blurrywhen its high-frequency spatial components are suppressed by the factor (here)max (1 − 6f 2 /fc2 , 0). The second program example illustrates a three-dimensionaltransform, where the three dimensions have different lengths. The third programexample is an example of convolution, as it might occur in a program to computethe potential generated by a three-dimensional distribution of sources.12.5 Fourier Transforms of Real Data in Two and Three Dimensions531int j;float fac,r,i,***data1,***data2,**speq1,**speq2,*sp1,*sp2;data1=f3tensor(1,N,1,N,1,N);data2=f3tensor(1,N,1,N,1,N);speq1=matrix(1,N,1,2*N);speq2=matrix(1,N,1,2*N);...*/rlft3(data1,speq1,N,N,N,1);FFT both input arrays.rlft3(data2,speq2,N,N,N,1);fac=2.0/(N*N*N);Factor needed to get normalized inverse.sp1 = &data1[1][1][1];sp2 = &data2[1][1][1];for (j=1;j<=N*N*N/2;j++) {Note how this can be made a single for-loop inr = sp1[0]*sp2[0] - sp1[1]*sp2[1];stead of three nested ones by usingi = sp1[0]*sp2[1] + sp1[1]*sp2[0];the pointers sp1 and sp2.sp1[0] = fac*r;sp1[1] = fac*i;sp1 += 2;sp2 += 2;}sp1 = &speq1[1][1];sp2 = &speq2[1][1];for (j=1;j<=N*N;j++) {r = sp1[0]*sp2[0] - sp1[1]*sp2[1];i = sp1[0]*sp2[1] + sp1[1]*sp2[0];sp1[0] = fac*r;sp1[1] = fac*i;sp1 += 2;sp2 += 2;}rlft3(data1,speq1,N,N,N,-1);Inverse FFT the product of the two FFTs./*...*/free_matrix(speq2,1,N,1,2*N);free_matrix(speq1,1,N,1,2*N);free_f3tensor(data2,1,N,1,N,1,N);free_f3tensor(data1,1,N,1,N,1,N);return 0;}To extend rlft3 to four dimensions, you simply add an additional (outer) nestedfor loop in i0, analogous to the present i1.
(Modifying the routine to do an arbitrarynumber of dimensions, as in fourn, is a good programming exercise for the reader.)CITED REFERENCES AND FURTHER READING:Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall).Swartztrauber, P. N. 1986, Mathematics of Computation, vol. 47, pp. 323–346.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)./*532Chapter 12.Fast Fourier Transform12.6 External Storage or Memory-Local FFTs#include <stdio.h>#include <math.h>#include "nrutil.h"#define KBF 128void fourfs(FILE *file[5], unsigned long nn[], int ndim, int isign)One- or multi-dimensional Fourier transform of a large data set stored on external media.
Oninput, ndim is the number of dimensions, and nn[1..ndim] contains the lengths of each dimension (number of real and imaginary value pairs), which must be powers of two. file[1..4]contains the stream pointers to 4 temporary files, each large enough to hold half of the data.The four streams must be opened in the system’s “binary” (as opposed to “text”) mode. Theinput data must be in C normal order, with its first half stored in file file[1], its secondhalf in file[2], in native floating point form. KBF real numbers are processed per bufferedread or write. isign should be set to 1 for the Fourier transform, to −1 for its inverse. Onoutput, values in the array file may have been permuted; the first half of the result is stored infile[3], the second half in file[4].
N.B.: For ndim > 1, the output is stored by columns,i.e., not in C normal order; in other words, the output is the transpose of that which would havebeen produced by routine fourn.{void fourew(FILE *file[5], int *na, int *nb, int *nc, int *nd);unsigned long j,j12,jk,k,kk,n=1,mm,kc=0,kd,ks,kr,nr,ns,nv;int cc,na,nb,nc,nd;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).Sometime in your life, you might have to compute the Fourier transform of a reallylarge data set, larger than the size of your computer’s physical memory.
In such a case,the data will be stored on some external medium, such as magnetic or optical tape or disk.Needed is an algorithm that makes some manageable number of sequential passes throughthe external data, processing it on the fly and outputting intermediate results to other externalmedia, which can be read on subsequent passes.In fact, an algorithm of just this description was developed by Singleton [1] very soonafter the discovery of the FFT. The algorithm requires four sequential storage devices, eachcapable of holding half of the input data. The first half of the input data is initially on onedevice, the second half on another.Singleton’s algorithm is based on the observation that it is possible to bit-reverse 2Mvalues by the following sequence of operations: On the first pass, values are read alternatelyfrom the two input devices, and written to a single output device (until it holds half the data),and then to the other output device.
On the second pass, the output devices become inputdevices, and vice versa. Now, we copy two values from the first device, then two valuesfrom the second, writing them (as before) first to fill one output device, then to fill a second.Subsequent passes read 4, 8, etc., input values at a time. After completion of pass M − 1,the data are in bit-reverse order.Singleton’s next observation is that it is possible to alternate the passes of essentiallythis bit-reversal technique with passes that implement one stage of the Danielson-Lanczoscombination formula (12.2.3). The scheme, roughly, is this: One starts as before with halfthe input data on one device, half on another. In the first pass, one complex value is readfrom each input device.
Two combinations are formed, and one is written to each of twooutput devices. After this “computing” pass, the devices are rewound, and a “permutation”pass is performed, where groups of values are read from the first input device and alternatelywritten to the first and second output devices; when the first input device is exhausted, thesecond is similarly processed. This sequence of computing and permutation passes is repeatedM − K − 1 times, where 2K is the size of internal buffer available to the program. Thesecond phase of the computation consists of a final K computation passes.
What distinguishesthe second phase from the first is that, now, the permutations are local enough to do in placeduring the computation. There are thus no separate permutation passes in the second phase.In all, there are 2M − K − 2 passes through the data.Here is an implementation of Singleton’s algorithm, based on [1]:.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
