Heath - Scientific Computing (523150), страница 100

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 100)

Hence, we expect the work required to be O(n2 )arithmetic operations instead of the O(n3 ) that would be required to factor and solve thelinear system. With the fast Fourier transform (FFT) algorithm, however, we can do muchbetter still, bringing the work down to O(n log2 n) arithmetic operations, as we will soonsee.12.2FFT AlgorithmBy taking advantage of certain symmetries and redundancies in the definition of the DFT, ashortcut algorithm can be developed for evaluating the DFT very efficiently. For illustration,consider the case n = 4. From the definition of the DFT we haveym =3Xxk wmk ,m = 0, .

. . , 3.k=0Writing out the four equations in full, we havey0 = x0 w0 + x1 w0 + x2 w0 + x3 w0 ,y1 = x0 w0 + x1 w1 + x2 w2 + x3 w3 ,y2 = x0 w0 + x1 w2 + x2 w4 + x3 w6 ,y3 = x0 w0 + x1 w3 + x2 w6 + x3 w9 .Noting from Fig. 12.1 thatw0 = w4 = 1,w2 = w6 = −1,w9 = w1 ,and regrouping, we obtain the four equationsy0 = (x0 + w0 x2 ) + w0 (x1 + w0 x3 ),y1 = (x0 − w0 x2 ) + w1 (x1 − w0 x3 ),y2 = (x0 + w0 x2 ) + w2 (x1 + w0 x3 ),y3 = (x0 − w0 x2 ) + w3 (x1 − w0 x3 ).We now observe that the transform can be computed with only 8 additions/subtractionsand 6 multiplications, instead of the expected (4 − 1) ∗ 4 = 12 additions and 42 = 16multiplications.

Actually, even fewer multiplications are required for this small case, sincew0 = 1, but we have tried to illustrate how the algorithm works in general. The main pointis that computing the DFT of the original 4-point sequence has been reduced to computingthe DFT of its two 2-point even and odd subsequences. This property holds in general:the DFT of an n-point sequence can be computed by breaking it into two DFTs of half thelength, provided n is even.Before stating the general FFT algorithm formally, we first consider a complementarydevelopment in terms of matrices that yields additional insight.

For a given integer n, we12.2. FFT ALGORITHM373define the Fourier matrix Fn by {Fn }mk = wmk . Thus, for example,1111111i −1 −i .F1 = 1, F2 =, F4 = 1 −11 −11 −1 1 −i −1iLet P4 be the permutation matrix10P4 = 000010010000,01and D2 be the diagonal matrix1D2 =00.iThen we have1111 1 −1i −i F2=F4 P4 = F211 −1 −11 −1 −iiD2 F2,−D2 F2i.e., F4 can be rearranged so that each block is a diagonally scaled version of F2 . Such ahierarchical splitting can be carried out at each level, provided the number of points is even.In general, Pn is the permutation that groups the even-numbered columns of Fn before theodd-numbered columns, and Dn/2 = diag(1, . .

. , wn/2 ). Thus, to apply Fn to a sequence oflength n, we need merely apply Fn/2 to its even and odd subsequences and scale the results,where necessary, by Dn/2 .This recursive divide-and-conquer approach to computing the DFT is formalized in thefollowing statement of the FFT algorithm, where we assume that n is a power of two:procedure fft (x, y, n, w)if n = 1 theny[0] = x[0]elsefor k = 0 to (n/2) − 1p[k] = x[2k]{ split into even ands[k] = x[2k + 1]odd subsequences }endfft (p, q, n/2, w2 )fft (s, t, n/2, w2 )for k = 0 to n − 1y[k] = q[k mod (n/2)] + wk t[k mod (n/2)]endendWe note the following points about the preceding algorithm:374CHAPTER 12. FAST FOURIER TRANSFORM• There are log2 n levels of recursion, each of which involves O(n) arithmetic operations,so the total cost is O(n log2 n).

If the weights wk are precomputed, the total number ofreal floating-point operations required by the FFT algorithm is 5n log2 n, compared with8n2 real floating-point operations for an ordinary complex matrix-vector product.• For clarity, separate arrays are used for the subsequences, but in fact the transform canbe computed in place using no additional storage.• The input sequence is assumed to be complex. If the input sequence is real, then additional symmetries in the DFT can be exploited to reduce the storage and operation countby half.• The output sequence is not produced in the natural order. Most FFT routines automatically allow for this by appropriately rearranging either the input or output sequence.This additional step does not affect the overall computational complexity, because thenecessary rearrangement (analogous to a sort) also costs O(n log2 n).• The FFT algorithm can be formulated using iteration rather than recursion, which isoften desirable for greater efficiency or when using a programming language that doesnot support recursion.Despite its name, the fast Fourier transform is an algorithm, not a transform.

It isa particular way (or family of ways) of computing the discrete Fourier transform of asequence in a very efficient manner. As we have seen, the DFT is defined in terms ofa matrix-vector product, whose straightforward evaluation would appear to require O(n2 )arithmetic operations. Use of the FFT algorithm reduces the work to only O(n log2 n), whichmakes an enormous practical difference in the time required to transform large sequences,as illustrated in Table 12.1.Table 12.1: Complexity of FFT algorithm versus matrix-vector multiplicationn n log2 nn2643844096128896163842562048655365124608262144102410240 1048576Owing to the similar form of the DFT and its inverse (they differ only in the sign ofthe exponent), the FFT algorithm can also be used to compute the inverse DFT efficiently.This ability to transform back and forth quickly between the time domain and the frequencydomain makes it practical to perform any computations or analysis that may be requiredin whichever domain is more convenient or efficient.12.2.1Limitations of the FFTAlthough the FFT algorithm has revolutionized many aspects of numerical computation, itis not always applicable or maximally efficient.

In particular, the input sequence is assumedto be:12.3. APPLICATIONS OF DFT375• Equally spaced• Periodic• A power of two in lengthThe first two of these properties follow from the definition of the DFT, whereas the third isrequired for maximal efficiency of the FFT algorithm. For example, the interpolant givenby the DFT, as a linear combination of sines and cosines, will necessarily be periodic, whichmeans that for a sequence of length n, we must have x0 = xn or, more generally, xk = xn+kfor any integer k (note that only x0 through xn−1 are actually specified).

Thus, some caremust be taken in applying the FFT algorithm to produce the most meaningful results asefficiently as possible. For instance, transforming a sequence that is not really periodic orpadding a sequence to make its length a power of two may introduce spurious noise andcomplicate the interpretation of the results.It is possible to define a “mixed-radix” FFT algorithm that does not require the numberof points n to be a power of two. The more general algorithm is still based on divide-andconquer; the sequence is not necessarily split exactly in half at each level, however, butrather by the smallest prime factor of the remaining sequence length.

For example, asequence of length 45 would be split into three subsequences of length 15, each of whichin turn would be split into three subsequences of length five. When a subsequence can besplit no further (i.e., when its length is prime), then its transform must be computed byconventional matrix-vector multiplication. The efficiency of such an algorithm depends onwhether n is a product of small primes (ideally a power of two, of course). If this is not thecase, then much of the computational advantage of the FFT may be lost.

For example, if nitself is a prime, then the original sequence cannot be split at all, and the “fast” algorithmthen becomes standard O(n2 ) matrix-vector multiplication.12.3Applications of DFTThe DFT is often of direct interest itself and is also useful as a computational tool thatprovides an efficient means for computing other quantities. The DFT is of direct interestin detecting periodicities or cycles in discrete data. Moreover, it can be used to removeunwanted periodicities. For example, to remove high-frequency noise from a sequence, onecan compute its DFT, set the high-frequency components of the transformed sequence tozero, then compute the inverse DFT of the modified sequence to get back into the originaldomain.As another example, weather data often contain two distinct cycles, diurnal and annual.One might want to remove one of these in order to study the other in isolation. Economicdata are also often “seasonally adjusted,” removing unwanted periodicities to reveal seculartrends.

Because of such uses, the DFT is of vital importance in many aspects of signalprocessing, such as digital filtering.Some computations are simpler or more efficient in the frequency domain than in thetime domain. Examples include the discrete circular convolution of two sequences u and vof length n,n−1X{u ? v}m =vk um−k , m = 0, 1, .

. . , n − 1,k=0376CHAPTER 12. FAST FOURIER TRANSFORMand related quantities such as the cross correlation of two sequences or the autocorrelationof a sequence with itself. In each case, the equivalent operation in the frequency domain issimply pointwise multiplication (in some cases with complex conjugation).For example, convolution is equivalent to multiplication by a circulant matrix , and sucha matrix is diagonalized by the DFT. Thus, if the convolution z of u and v is given by  uun−1 un−2 · · ·u1   v0 z00u0un−1 · · ·u2   v1  z1   u 1 .   . .

Характеристики

Список файлов книги