532 Chapter 12.Fast Fourier Transform 12.6 External Storage or Memory-Local FFTs Sometime in your life,you might have to compute the Fourier transform of a really large 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 through the external data,processing it on the fly and outputting intermediate results to other external media,which can be read on subsequent passes. In fact,an algorithm of just this description was developed by Singleton [1]very soon after the discovery of the FFT.The algorithm requires four sequential storage devices,each capable of holding half of the input data.The first half of the input data is initially on one 8 device,the second half on another Singleton's algorithm is based on the observation that it is possible to bit-reverse2 values by the following sequence of operations:On the first pass,values are read alternately from 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 input devices,and vice versa.Now,we copy nwo values from the first device,then two values from 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 essentially this bit-reversal technique with passes that implement one stage of the Danielson-Lanczos combination formula (12.2.3).The scheme,roughly,is this:One starts as before with half 导不 Press. the input data on one device,half on another.In the first pass,one complex value is read ART from each input device.Two combinations are formed,and one is written to each of two output devices.After this“computing”pass,the devices are rewound,anda“permutation” Programs pass is performed,where groups of values are read from the first input device and alternately written to the first and second output devices;when the first input device is exhausted,the second is similarly processed.This sequence of computing and permutation passes is repeated SCIENTIFIC M-K-1 times,where 2 is the size of internal buffer available to the program.The second phase of the computation consists of a final K computation passes.What distinguishes the second phase from the first is that,now,the permutations are local enough to do in place during 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]: (ISBN #include 21 #include #include "nrutil.h" #define KBF 128 Numerical Recipes 43108 void 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.On (outside input,ndim is the number of dimensions,and nn[1..ndim]contains the lengths of each di- mension (number of real and imaginary value pairs),which must be powers of two.file[1..4] North Software. 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.The input data must be in C normal order,with its first half stored in file file[1],its second half in file[2],in native floating point form.KBF real numbers are processed per buffered read or write.isign should be set to 1 for the Fourier transform,to -1 for its inverse.On output,values in the array file may have been permuted;the first half of the result is stored in file [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 have been 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;
532 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). 12.6 External Storage or Memory-Local FFTs Sometime in your life, you might have to compute the Fourier transform of a really large 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 through the external data, processing it on the fly and outputting intermediate results to other external media, which can be read on subsequent passes. In fact, an algorithm of just this description was developed by Singleton [1] very soon after the discovery of the FFT. The algorithm requires four sequential storage devices, each capable of holding half of the input data. The first half of the input data is initially on one device, the second half on another. Singleton’s algorithm is based on the observation that it is possible to bit-reverse 2M values by the following sequence of operations: On the first pass, values are read alternately from 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 input devices, and vice versa. Now, we copy two values from the first device, then two values from 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 essentially this bit-reversal technique with passes that implement one stage of the Danielson-Lanczos combination formula (12.2.3). The scheme, roughly, is this: One starts as before with half the input data on one device, half on another. In the first pass, one complex value is read from each input device. Two combinations are formed, and one is written to each of two output 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 alternately written to the first and second output devices; when the first input device is exhausted, the second is similarly processed. This sequence of computing and permutation passes is repeated M − K − 1 times, where 2K is the size of internal buffer available to the program. The second phase of the computation consists of a final K computation passes. What distinguishes the second phase from the first is that, now, the permutations are local enough to do in place during 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]: #include #include #include "nrutil.h" #define KBF 128 void 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. On input, 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. The input data must be in C normal order, with its first half stored in file file[1], its second half in file[2], in native floating point form. KBF real numbers are processed per buffered read or write. isign should be set to 1 for the Fourier transform, to −1 for its inverse. On output, values in the array file may have been permuted; the first half of the result is stored in file[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 have been 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;
12.6 External Storage or Memory-Local FFTs 533 float tempr,tempi,*afa,*afb,*afc; double wr,wi,wpr,wpi,wtemp,theta; static int mate[5]-(0,2,1,4,3}; afa=vector(1,KBF); afb=vector(1,KBF); afc=vector(1,KBF); for (j=1;j1: granted fori kd=KBF >>1; ks=n; .com or call 1-800-872- (including this one) fourew (file,&na,&nb,&nc,&nd); internet The first phase of the transform starts here. for(;;)[ Start of the computing pass. theta=1s1gn*3.141592653589793/(n/mm); wtemp=sin(0.5*theta); wpr =-2.0*wtemp*wtemp; wpi=sin(theta); r=1.0 -7423(North America w1=0.0; to any server computer,is strictly prohibited. t users to make one paper copy for their 1988-1992 by Cambridge University Press.Programs Copyright (C)1988-1992 by Numerical Recipes from NUMERICAL RECIPES IN THE mm>>=1; for(j12=1;j12>=1; while (jk ==1){ mm=ni
12.6 External Storage or Memory-Local FFTs 533 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). float tempr,tempi,*afa,*afb,*afc; double wr,wi,wpr,wpi,wtemp,theta; static int mate[5] = {0,2,1,4,3}; afa=vector(1,KBF); afb=vector(1,KBF); afc=vector(1,KBF); for (j=1;j> 1; kd=KBF >> 1; ks=n; fourew(file,&na,&nb,&nc,&nd); The first phase of the transform starts here. for (;;) { Start of the computing pass. theta=isign*3.141592653589793/(n/mm); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; mm >>= 1; for (j12=1;j12>= 1; while (jk == 1) { mm=n;
534 Chapter 12. Fast Fourier Transform jk=nn [++nv] k8>>=1; if (ks KBF){ for(j12=1;j12>=1; ks=kd; copy for their kd>>=1: for(j12=1:j12KBF){ cc=fwrite(&afa[1],sizeof(float),KBF,file [nc]); if (cc !KBF)nrerror("write error in fourfs"); cc=fwrite(&afb[1],sizeof(float),KBF,file [nd]); if (cc !KBF)nrerror("write error in fourfs"); j=1;
534 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). jk=nn[++nv]; } ks >>= 1; if (ks > KBF) { for (j12=1;j12>= 1; ks=kd; kd >>= 1; for (j12=1;j12 KBF) break; else k=kk+ks; } if (j > KBF) { cc=fwrite(&afa[1],sizeof(float),KBF,file[nc]); if (cc != KBF) nrerror("write error in fourfs"); cc=fwrite(&afb[1],sizeof(float),KBF,file[nd]); if (cc != KBF) nrerror("write error in fourfs"); j=1; }
12.6 External Storage or Memory-Local FFTs 535 na-mate[na]; fourew(file,&na,&nb,&nc,&nd); jk>>=1; if (jk 1)continue; mm=n; do if (nv ndim)jk=nn[++nv]; else free_vector(afc,1,KBF); free_vector(afb,1,KBF); http://www.nr. Permission is read able files free_vector(afa,1,KBF); return: 2 Jh11e(jk==1); .com or call (including this one) granted for 19881992 11-800-872 #include #define SWAP(a,b)ftemp=(a);(a)=(b);(b)=ftemp n NUMERICAL RECIPES IN C: void fourew(FILE *file[5],int *na,int *nb,int *nc,int *nd) Utility used by fourfs.Rewinds and renumbers the four files. -7423(North America 电r:1s t users to make one paper by Cambridge University Press. THE int i; FILE *ftemp; 是 ART for (i=1;i<=4;i++)rewind(file[i]); Programs SWAP(file[2],file[4]); copy for thei strictly proh SwAP(f11e[1],f11e[3]); *na=3; *nb=4; Copyright (C) *nC=11 *nd=2; 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN For one-dimensional data,Singleton's algorithm produces output in exactly the same order as a standard FFT(e.g.,four1).For multidimensional data,the output is the transpose of v@cam the conventional arrangement (e.g.,the output of fourn).This peculiarity,which is intrinsic to 10-:6211 the method,is generally only a minor inconvenience.For convolutions,one simply computes the component-by-component product of two transforms in their nonstandard arrangement. 43108 and then does an inverse transform on the result.Note that,if the lengths of the different dimensions are not all the same,then you must reverse the order of the values in nn [1.ndim] (thus giving the transpose dimensions)before performing the inverse transform.Note also (outside that,just like fourn,performing a transform and then an inverse results in multiplying the original data by the product of the lengths of all dimensions. North Software. We leave it as an exercise for the reader to figure out how to reorder fourfs's output into normal order,taking additional passes through the externally stored data.We doubt that Ame such reordering is ever really needed. You will likely want to modify fourfs to fit your particular application.For example, visit website machine as written,KBF =2 plays the dual role of being the size of the internal buffers,and the record size of the unformatted reads and writes.The latter role limits its size to that allowed by your machine's 1/O facility.It is a simple matter to perform multiple reads for a much larger KBF,thus reducing the number of passes by a few. Another modification of fourfs would be for the case where your virtual memory machine has sufficient address space,but not sufficient physical memory,to do an efficient FFT by the conventional algorithm (whose memory references are extremely nonlocal).In that case,you will need to replace the reads,writes,and rewinds by mappings of the arrays
12.6 External Storage or Memory-Local FFTs 535 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). } na=mate[na]; } fourew(file,&na,&nb,&nc,&nd); jk >>= 1; if (jk > 1) continue; mm=n; do { if (nv #define SWAP(a,b) ftemp=(a);(a)=(b);(b)=ftemp void fourew(FILE *file[5], int *na, int *nb, int *nc, int *nd) Utility used by fourfs. Rewinds and renumbers the four files. { int i; FILE *ftemp; for (i=1;i<=4;i++) rewind(file[i]); SWAP(file[2],file[4]); SWAP(file[1],file[3]); *na=3; *nb=4; *nc=1; *nd=2; } For one-dimensional data, Singleton’s algorithm produces output in exactly the same order as a standard FFT (e.g., four1). For multidimensional data, the output is the transpose of the conventional arrangement (e.g., the output of fourn). This peculiarity, which is intrinsic to the method, is generally only a minor inconvenience. For convolutions, one simply computes the component-by-component product of two transforms in their nonstandard arrangement, and then does an inverse transform on the result. Note that, if the lengths of the different dimensions are not all the same, then you must reverse the order of the values in nn[1..ndim] (thus giving the transpose dimensions) before performing the inverse transform. Note also that, just like fourn, performing a transform and then an inverse results in multiplying the original data by the product of the lengths of all dimensions. We leave it as an exercise for the reader to figure out how to reorder fourfs’s output into normal order, taking additional passes through the externally stored data. We doubt that such reordering is ever really needed. You will likely want to modify fourfs to fit your particular application. For example, as written, KBF ≡ 2K plays the dual role of being the size of the internal buffers, and the record size of the unformatted reads and writes. The latter role limits its size to that allowed by your machine’s I/O facility. It is a simple matter to perform multiple reads for a much larger KBF, thus reducing the number of passes by a few. Another modification of fourfs would be for the case where your virtual memory machine has sufficient address space, but not sufficient physical memory, to do an efficient FFT by the conventional algorithm (whose memory references are extremely nonlocal). In that case, you will need to replace the reads, writes, and rewinds by mappings of the arrays
536 Chapter 12 Fast Fourier Transform afa,afb,and afc into your address space.In other words,these arrays are replaced by references to a single data array,with offsets that get modified wherever fourfs performs an I/O operation.The resulting algorithm will have its memory references local within blocks of size KBF.Execution speed is thereby sometimes increased enormously,albeit at the cost of requiring twice as much virtual memory as an in-place FFT. CITED REFERENCES AND FURTHER READING: Singleton,R.C.1967,IEEE Transactions on Audio and Electroacoustics,vol.AU-15,pp.91-97. [1] http://www.n read able files Permission is Sample page Oppenheim.A.V.,and Schafer.R.W.1989,Discrete-Time Signa/Processing (Englewood Cliffs NJ:Prentice-Hall).Chapter 9. .com or call 1-800-872-7423(North America only)orsend email to directcustserv@cambridge.org (outside North America). granted for internet users to make one paper copy for their Copyright (C)1988-1992 by Cambridge University Press.Programs Copyright(C)1988-1992 by Numerical Recipes from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) Software
536 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). afa, afb, and afc into your address space. In other words, these arrays are replaced by references to a single data array, with offsets that get modified wherever fourfs performs an I/O operation. The resulting algorithm will have its memory references local within blocks of size KBF. Execution speed is thereby sometimes increased enormously, albeit at the cost of requiring twice as much virtual memory as an in-place FFT. CITED REFERENCES AND FURTHER READING: Singleton, R.C. 1967, IEEE Transactions on Audio and Electroacoustics, vol. AU-15, pp. 91–97. [1] Oppenheim, A.V., and Schafer, R.W. 1989, Discrete-Time Signal Processing (Englewood Cliffs, NJ: Prentice-Hall), Chapter 9
