f12-6

12.6ExternalStorageorMemory-LocalFFTs 525Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).12.6 External Storage or Memory-Local FFTs Sometime in your life, you might have to compute the Fourier transform of a really largedata 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 soon after 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 one device, the second half on another. Singleton’s algorithm is based on the observation that it is possible to bit-reverse 2 M 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 inputdevices, and vice versa. Now, we copy twovalues from the first device, then twovalues 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-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, thesecondissimilarlyprocessed. ThissequenceofcomputingandpermutationpassesisrepeatedM−K−1times, where 2 Kis the size of internal buffer available to the program. The secondphaseofthecomputationconsistsofafinal Kcomputationpasses. Whatdistinguishes the 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−2passes through the data. Here is an implementation of Singleton’s algorithm, based on [1]: SUBROUTINE fourfs(iunit,nn,ndim,isign) INTEGER ndim,nn(ndim),isign,iunit(4),KBFPARAMETER (KBF=128) C USES fourew One- or multi-dimensional Fourier transform of a large data set stored on external media.On input, ndimis the number of dimensions, and nn(1:ndim) contains the lengths of each dimension (number of complex values), which must be powers of two. iunit(1:4) contains theunit numbers of4 sequential files,each largeenough to hold halfofthe data.The four units must be opened for FORTRANunformatted access. The input data must be inFORTRANnormal order, with its first half stored on unit iunit(1), its second half on iunit(2), in unformatted form, with KBFreal numbers per record. isignshould be set to 1 forthe Fourier transform, to −1forits inverse. Onoutput, valuesinthearray iunit may have been permuted; the first half of the result is stored on iunit(3), the second halfon iunit(4). N.B.: For ndim >1,theoutputisstoredbyrows,i.e., notinFORTRAN normal order; in other words, the output is the transpose of that which would have been produced by routine fourn. INTEGER j,j12,jk,k,kk,n,mm,kc,kd,ks,kr,nr,ns,nv,jx, * mate(4),na,nb,nc,nd REAL tempr,tempi,afa(KBF),afb(KBF),afc(KBF) DOUBLE PRECISION wr,wi,wpr,wpi,wtemp,thetaSAVE mate DATA mate /2,1,4,3/ 526 Chapter12. FastFourierTransformSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).n=1 do11j=1,ndim n=n*nn(j) if (nn(j).le.1) * pause ’invalid dimension or wrong ndim in fourfs’ enddo 11 nv=ndim jk=nn(nv)mm=n ns=n/KBF nr=ns/2kc=0kd=KBF/2 ks=n call fourew(iunit,na,nb,nc,nd) The first phase of the transform starts here. 1 continue Start of the computing pass. theta=3.141592653589793d0/(isign*n/mm)wpr=-2.d0*sin(0.5d0*theta)**2wpi=sin(theta) wr=1.d0 wi=0.d0mm=mm/2 do 13j12=1,2 kr=0 2 continue read (iunit(na)) (afa(jx),jx=1,KBF) read (iunit(nb)) (afb(jx),jx=1,KBF) do12j=1,KBF,2 tempr=sngl(wr)*afb(j)-sngl(wi)*afb(j+1) tempi=sngl(wi)*afb(j)+sngl(wr)*afb(j+1) afb(j)=afa(j)-temprafa(j)=afa(j)+temprafb(j+1)=afa(j+1)-tempi afa(j+1)=afa(j+1)+tempi enddo 12 kc=kc+kdif (kc.eq.mm) then kc=0 wtemp=wrwr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi endifwrite (iunit(nc)) (afa(jx),jx=1,KBF)write (iunit(nd)) (afb(jx),jx=1,KBF) kr=kr+1 if (kr.lt.nr) goto 2if(j12.eq.1.and.ks.ne.n.and.ks.eq.KBF) then na=mate(na) nb=na endifif (nr.eq.0) goto 3 enddo 13 3 call fourew(iunit,na,nb,nc,nd) Start of the permutation pass. jk=jk/2 4 if (jk.eq.1) then mm=nnv=nv-1jk=nn(nv) goto 4 endifks=ks/2if (ks.gt.KBF) then do 16j12=1,2 12.6ExternalStorageorMemory-LocalFFTs 527Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do15kr=1,ns,ks/KBF do14k=1,ks,KBF read (iunit(na)) (afa(jx),jx=1,KBF) write (iunit(nc)) (afa(jx),jx=1,KBF) enddo 14 nc=mate(nc) enddo 15 na=mate(na) enddo 16 call fourew(iunit,na,nb,nc,nd)goto 1 else if (ks.eq.KBF) then nb=na goto 1 endif continue j=1 The second phase of the transform starts here. Now, the remaining permutations are suffi-ciently local to be done in place. 5 continue theta=3.141592653589793d0/(isign*n/mm) wpr=-2.d0*sin(0.5d0*theta)**2wpi=sin(theta) wr=1.d0 wi=0.d0mm=mm/2ks=kd kd=kd/2 do 18j12=1,2 do17kr=1,ns read (iunit(na)) (afc(jx),jx=1,KBF) kk=1k=ks+1 6 continue tempr=sngl(wr)*afc(kk+ks)-sngl(wi)*afc(kk+ks+1) tempi=sngl(wi)*afc(kk+ks)+sngl(wr)*afc(kk+ks+1)afa(j)=afc(kk)+temprafb(j)=afc(kk)-tempr afa(j+1)=afc(kk+1)+tempi afb(j+1)=afc(kk+1)-tempij=j+2 kk=kk+2 if (kk.lt.k) goto 6kc=kc+kdif (kc.eq.mm) then kc=0 wtemp=wrwr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi endifkk=kk+ksif (kk.le.KBF) then k=kk+ks goto 6 endif if (j.gt.KBF) then write (iunit(nc)) (afa(jx),jx=1,KBF)write (iunit(nd)) (afb(jx),jx=1,KBF)j=1 endif enddo 17 na=mate(na) enddo 18 call fourew(iunit,na,nb,nc,nd) 528 Chapter12. Fast FourierTransformSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).jk=jk/2 if (jk.gt.1) goto 5 mm=n 7 if (nv.gt.1) then nv=nv-1 jk=nn(nv) if (jk.eq.1) goto 7goto 5 endif return END SUBROUTINE fourew(iunit,na,nb,nc,nd) INTEGER na,nb,nc,nd,iunit(4),ii Utility used by fourfs. Rewinds and renumbers the four files. do11ii=1,4 rewind(unit=iunit(ii)) enddo 11 ii=iunit(2) iunit(2)=iunit(4)iunit(4)=ii ii=iunit(1) iunit(1)=iunit(3)iunit(3)=ii na=3 nb=4nc=1nd=2 return END For one-dimensional data, Singleton’s algorithm produces output in exactly the same orderasastandardFFT(e.g., four1). Formultidimensionaldata,theoutputisthe transpose of theconventionalarrangement(e.g.,theoutputof fourn). Thispeculiarity,whichisintrinsicto the method, is generally only a minor inconvenience. For convolutions, one simply computesthe 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 differentdimensions arenotallthesame,thenyou mustreversetheorder ofthevaluesin nn(1:ndim) (thus giving the transpose dimensions) before performing the inverse transform. Note alsothat, 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 thatsuch reordering is ever really needed. You will likely want to modify fourfsto fit your particular application. For example, as written, KBF≡2 Kplays 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 allowedby your machine’s I/O facility. It is a simple matter to perform multiple reads for a muchlarger KBF, thus reducing the number of passes by a few. Another modification of fourfswould be for the case where your virtual memory machine has sufficient address space, but not sufficient physical memory, to do an efficientFFT by the conventional algorithm (whose memory references are extremely nonlocal). Inthat case, you will need to replace the reads, writes, and rewinds by mappings of the arraysafa,afb, and afcinto your address space. In other words, these arrays are replaced by references to a single data array, with offsets that get modified wherever fourfsperforms an I/O operation. The resulting algorithm will have its memory references local within blocksof 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. 12.6ExternalStorageorMemory-LocalFFTs 529Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Singleton, R.C. 1967, IEEE Transactions onAudioand Electroacoustics , vol. AU-15, pp.91–97. [1] Oppenheim, A.V., andSchafer, R.W. 1989, Discrete-Time Signal Processing (EnglewoodCliffs, NJ: Prentice-Hall), Chapter 9.