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

nn MUSTbe an integer power of 2 (this is not checked for!).{unsigned long n,mmax,m,j,istep,i;double wtemp,wr,wpr,wpi,wi,theta;Double precision for the trigonometfloat tempr,tempi;ric recurrences.n=nn << 1;j=1;for (i=1;i<n;i+=2) {This is the bit-reversal section of theif (j > i) {routine.SWAP(data[j],data[i]);Exchange the two complex numbers.SWAP(data[j+1],data[i+1]);}m=n >> 1;while (m >= 2 && j > m) {j -= m;m >>= 1;}j += m;}Here begins the Danielson-Lanczos section of the routine.mmax=2;while (n > mmax) {Outer loop executed log2 nn times.istep=mmax << 1;theta=isign*(6.28318530717959/mmax);Initialize the trigonometric recurrence.wtemp=sin(0.5*theta);wpr = -2.0*wtemp*wtemp;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).length of the real array (data[1..2*nn]) is 2 times nn, with each complex valueoccupying two consecutive locations. In other words, data[1] is the real part off0 , data[2] is the imaginary part of f0 , and so on up to data[2*nn-1], whichis the real part of fN−1 , and data[2*nn], which is the imaginary part of fN−1 .The FFT routine gives back the Fn ’s packed in exactly the same fashion, as nncomplex numbers.The real and imaginary parts of the zero frequency component F0 are in data[1]and data[2]; the smallest nonzero positive frequency has real and imaginary parts indata[3] and data[4]; the smallest (in magnitude) nonzero negative frequency hasreal and imaginary parts in data[2*nn-1] and data[2*nn].

Positive frequenciesincreasing in magnitude are stored in the real-imaginary pairs data[5], data[6]up to data[nn-1], data[nn]. Negative frequencies of increasing magnitude arestored in data[2*nn-3], data[2*nn-2] down to data[nn+3], data[nn+4].Finally, the pair data[nn+1], data[nn+2] contain the real and imaginary parts ofthe one aliased point that contains the most positive and the most negative frequency.You should try to develop a familiarity with this storage arrangement of complexspectra, also shown in Figure 12.2.2, since it is the practical standard.This is a good place to remind you that you can also use a routine like four1without modification even if your input data array is zero-offset, that is has the rangedata[0..2*nn-1].

In this case, simply decrement the pointer to data by one whenfour1 is invoked, e.g., four1(data-1,1024,1);. The real part of f0 will now bereturned in data[0], the imaginary part in data[1], and so on. See §1.2.508Chapter 12.2imag3real4imagt=0t=∆2N − 3real2N − 2imag2N − 1real2N(a)imag1real2imag3real4imagN−1realNimagN+1realN+2imagN+3realN+4imag2N − 1real2Nimagf=0f=1N∆f=N/2 − 1N∆f=±1(combination)2∆f= −N/2 − 1N∆f=−1N∆t = (N − 2)∆t = (N − 1)∆(b)Figure 12.2.2.Input and output arrays for FFT. (a) The input array contains N (a power of 2)complex time samples in a real array of length 2N , with real and imaginary parts alternating.

(b) Theoutput array contains the complex Fourier spectrum at N values of frequency. Real and imaginary partsagain alternate. The array starts with zero frequency, works up to the most positive frequency (whichis ambiguous with the most negative frequency). Negative frequencies follow, from the second-mostnegative up to the frequency just below zero.wpi=sin(theta);wr=1.0;wi=0.0;for (m=1;m<mmax;m+=2) {for (i=m;i<=n;i+=istep) {j=i+mmax;tempr=wr*data[j]-wi*data[j+1];tempi=wr*data[j+1]+wi*data[j];data[j]=data[i]-tempr;data[j+1]=data[i+1]-tempi;data[i] += tempr;data[i+1] += tempi;}wr=(wtemp=wr)*wpr-wi*wpi+wr;wi=wi*wpr+wtemp*wpi+wi;}mmax=istep;Here are the two nested inner loops.This is the Danielson-Lanczos formula:Trigonometric recurrence.}}(A double precision version of four1, named dfour1, is used by the routine mpmulin §20.6.

You can easily make the conversion, or else get the converted routinefrom the Numerical Recipes diskette.)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).real array of length 2 Nrealreal array of length 2 N1Fast Fourier Transform12.2 Fast Fourier Transform (FFT)509Other FFT AlgorithmsSample 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).We should mention that there are a number of variants on the basic FFT algorithmgiven above. As we have seen, that algorithm first rearranges the input elementsinto bit-reverse order, then builds up the output transform in log2 N iterations.

Inthe literature, this sequence is called a decimation-in-time or Cooley-Tukey FFTalgorithm. It is also possible to derive FFT algorithms that first go through a set oflog2 N iterations on the input data, and rearrange the output values into bit-reverseorder. These are called decimation-in-frequency or Sande-Tukey FFT algorithms. Forsome applications, such as convolution (§13.1), one takes a data set into the Fourierdomain and then, after some manipulation, back out again. In these cases it is possibleto avoid all bit reversing. You use a decimation-in-frequency algorithm (without itsbit reversing) to get into the “scrambled” Fourier domain, do your operations there,and then use an inverse algorithm (without its bit reversing) to get back to the timedomain. While elegant in principle, this procedure does not in practice save muchcomputation time, since the bit reversals represent only a small fraction of an FFT’soperations count, and since most useful operations in the frequency domain requirea knowledge of which points correspond to which frequencies.Another class of FFTs subdivides the initial data set of length N not all theway down to the trivial transform of length 1, but rather only down to some othersmall power of 2, for example N = 4, base-4 FFTs, or N = 8, base-8 FFTs.

Thesesmall transforms are then done by small sections of highly optimized coding whichtake advantage of special symmetries of that particular small N . For example, forN = 4, the trigonometric sines and cosines that enter are all ±1 or 0, so manymultiplications are eliminated, leaving largely additions and subtractions. Thesecan be faster than simpler FFTs by some significant, but not overwhelming, factor,e.g., 20 or 30 percent.There are also FFT algorithms for data sets of length N not a power oftwo. They work by using relations analogous to the Danielson-Lanczos Lemma tosubdivide the initial problem into successively smaller problems, not by factors of2, but by whatever small prime factors happen to divide N . The larger that thelargest prime factor of N is, the worse this method works.

If N is prime, then nosubdivision is possible, and the user (whether he knows it or not) is taking a slowFourier transform, of order N 2 instead of order N log2 N . Our advice is to stay clearof such FFT implementations, with perhaps one class of exceptions, the WinogradFourier transform algorithms. Winograd algorithms are in some ways analogous tothe base-4 and base-8 FFTs. Winograd has derived highly optimized codings fortaking small-N discrete Fourier transforms, e.g., for N = 2, 3, 4, 5, 7, 8, 11, 13, 16.The algorithms also use a new and clever way of combining the subfactors.