Thompson - Computing for Scientists and Engineers (523188), страница 61
Текст из файла (страница 61)
Consider, for example, y4, for which its index210satisfies 4 - 1 = 3 = 2 0 + 2 1 + 2 1, which is 011 in binary representation.After the transform is finished y4 has been moved into the position occupied by y8.210Upon noting that 8 - 1 = 7 = 2 1 + 2 1 + 2 1, which is 110 in binary, it seemsthat we have bit-reversed 0 11. We can formalize this example, as in the followingexercise.Exercise 9.11Consider the binary representation of the integers between 0 and 2v - 1, therange of the index k - 1 for the FFT coefficients Ck. Such an integer can berepresented by(9.35)in which each ji is either 0 or 1. The bit-reversed integer is then jR, where(9.36)(a) Show that jR can be generated similarly to Horner’s algorithm (Section 4.1)as follows:j T = j; j R = 0;Iterate the next steps in integer arithmeticFor i = 1 to v{jD = jT/2;ji -1 = jT - 2jD;JR = 2jR + ji-1 ;jT = jD;}End of for loop(b) Verify this algorithm for the example k = 4, considered above.
n9.3THE FAST FOURIER TRANSFORM ALGORITHM333The final step in an FFT algorithm is therefore to bit-reverse the labels on the coefficients according to this scheme. Note that it is k - 1, rather than k, that has tobe reversed. This is because we used a labeling scheme for the frequencies in whichthey range from the first harmonic (k = 1) to the (N - 1) th harmonic, as in (9.26).Other derivations of the FFT may use a labeling from k = 0 to k = N - 1, thusavoiding this slight awkwardness. A related presentation of the DFT, lending itselfto derivation of the FFT, is given in the article by Peters.Having derived an FFT algorithm, was it worth our effort and will it be worththe effort to code it?Efficiency of FFT and conventional transformsWe now discuss the efficiency of the FFT by comparing the times for calculating thediscrete Fourier transform of N real data points using the fast (FFT) and conventional (CFT) algorithms.
Our comparison is intended to estimate the dependence ofthe times on N, independent of details of the computer system used.For the FFT the total number of arithmetic steps of the form (9.30) is roughlyproportional to Nv = N log2 N. A more-detailed analysis, as presented by Brigham, shows that the time, T (FFT ), for the fast Fourier transform of N points inour radix-2 algorithm is about(9.37)where Ta is the computer time for a typical operation in real arithmetic, such as multiplication or addition.For the CFT algorithm the time, T(CFT), to compute N coefficients involvingN operations must scale as N2.
It is estimated by Brigham to be about(9.38)Graphical comparison of the two algorithms is made in Figure 9.7.Exercise 9.12(a) On my workstation I programmed and ran the FFT algorithm and programgiven in Section 9.7. I found that for N = 128 the running time wasT (FFT) = 2.1 seconds, for N = 256 T (FFT) = 4.4 seconds, while forN = 1024 T (FFT) = 20.1 seconds. From these data, estimate a range forTa, and show that the FFT time for a data sample of N = 8192 = 213 points isexpected to be about 3.5 minutes.(b) Use the estimates of Ta from (a) and (9.37) and (9.38) to show that a CFTis estimated to take the impractically long time of about 47 hours of continuousrunning to transform 8192 data points.
Thus, the CFT would be about 800times slower than the FFT. nThis exercise should help you to appreciate the practicality of the FFT for largeN values, together with the impracticality of the CFT for the same values of N. For334DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESFIGURE 9.7 Comparative times of FFT and CFT according to (9.37) and (9.38).small values of N, which we omit from the graph, there will be more startup time forthe FFT, and some time (proportional to N) will be required for the bit reversal. Already for N = 64 the figure shows that the FFT is predicted to be one order of magnitude faster than the CFT, according to the estimates (9.37) and (9.38).
You canuse the program in Section 9.7 to compare the timing for the FFT on your owncomputer.Further reading on the FFT algorithm, including its generalization to transformsin two dimensions, which are widely used in imaging science, is provided by themonograph of Brigham, in Oppenheim and Schafer, in the book on signal processing by Roberts and Mullis. Transformation into real sums of cosines and sines,which is particularly convenient for data, is called a Hartley transform, as discussedin Bracewell’s book of this title.
We apply the FFT algorithm in Section 9.7 to develop a program which is then applied in Section 9.8 to the Fourier analysis of anelectroencephalogram.9.4 FOURIER SERIES: HARMONIC APPROXIMATIONSIt is often convenient to obtain Fourier expansion coefficients by the procedure of integration, rather than by the summation used in the discrete Fourier transform considered in Sections 9.2 and 9.3. The main reason is that integrals are often easier toevaluate analytically than are sums. Integration requires, however, that the functionbe defined at all points in the range of integration.
Therefore, Fourier series aremost appropriate for use with y defined by a function rather than the yj being discrete data.From discrete transforms to seriesThe transition from the discrete Fourier transform to the Fourier series is made byletting the number of x values in the DFT, N, become indefinitely large in such away that9.4 FOURIER SERIES: HARMONIC APPROXIMATIONS335(9.39)and that the x increments become infinitesimalThis makes the range of x from 0 toThe Fourier expansion (9.7) thus becomes(9.41)The Fourier series coefficients are obtained by applying the same limiting process to(9.6). Thus, for k = -N to N, the Fourier series coefficients are(9.42)Exercise 9.13Work through in detail the algebraic steps from the discrete Fourier transform tothe Fourier series.
nThe more usual trigonometric form of the Fourier series is obtained by expanding the exponential in the preceding formulas into sines and cosines. Thus, theFourier series expansion in terms of cosines and sines is written(9.43)You will find directly that the Fourier coefficients of the cosine terms are(9.44)The Fourier coefficients of the sine terms are given by(9.45)We have made the conventional division of the normalizing factors, putting theuniformly in the coefficients, (9.44) and (9.45), and therefore omitting it from theexpansion (9.43).336DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESExercise 9.14Derive (9.44) and (9.45) from the corresponding exponential form of the Fourierseries, (9.41) and (9.42) by following the indicated steps. nNow that we have the formulas for the Fourier series amplitudes, how are theyto be interpreted?Interpreting Fourier coefficientsConsider equations (9.44) and (9.45) for the Fourier coefficients.
The integrals forak and bk indicate to what degree y(x) varies like the oscillatory functions cos (kx)and sin (kx), respectively. Both of these functions have a period of oscillation ofThus, as k increases, the more rapidly varying parts of y(x) are the most important in determining the corresponding ak and bk. In musical notation, k = 1gives the fundamental frequency, and all the higher frequencies are its harmonics.An alternative name for the Fourier series expansion is the “harmonic expansion.”Since, as (9.44) shows, the coefficient a0 is just the average value of y over theinterval 0 toit is often assumed that this average value has been subtracted outbefore y is Fourier analyzed.
Thus, the formulas for the ak and bk become symmetric, and there need be no special treatment for the DC (k = 0) amplitude. The treatment of Fourier series then becomes closer to the usual treatment of the discreteFourier transform.Exercise 9.15Show that for the Fourier series (9.43), with the expansion in the range x = 0to x =one predicts a periodicity of y ofthat is,= y(x). nThis exercise shows that the Fourier series expansion that we have derived abovegenerates a function with the same period as the range of x values used to generatethe coefficients. This is the same property as found for the discrete Fourier transform in Section 9.2, where N data values produced N coefficients. The generalization to Fourier series for arbitrary intervals is made in the next subsection.Fourier series for arbitrary intervalsSuppose that we have a function of x and want to compute its Fourier series by using the range from x- to x+.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
