Thompson - Computing for Scientists and Engineers (523188), страница 60
Текст из файла (страница 60)
The complex amplitude may beexpressed as(9.21)where the relation between frequencyand number k0 is(9.22)The second line in (9.21) is obtained either by applying L’Hôpital’s rule from differential calculus to the first line or from (9.14) directly. Notice that this transform hasno intrinsic dependence on the stepsize h used to compute it, except through thescale factor in (9.21), which was inserted in (9.14) only to produce the appropriatelimit for the integral transform in Section 10.2, and through the conversion (9.15)between number and frequency. Therefore, when discussing the oscillator in thefollowing we set h = 1.Exercise 9.7(a) Derive (9.21) from (9.14) by the route indicated.(b) Apply L’Hôpital’s rule to the first line of (9.21) to produce the second linein this equation.
n326DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESThe symmetry of the discrete transform for real functions, (9.10), does not holdfor the oscillator, but there is a related symmetry, namely(9.23)in which we have extended the notation so that the quantity in parentheses indicatesthe choice of “resonance” frequency. When k0 is an integer or halfinteger the oscillator transform simplifies considerably, as can be shown by substituting the valuesof the circular functions at multiples of or ofFor k0 an integer(9.24)so that the transform is discontinuous at k0.
For k0 a half-integer, such as 1/2 or15/2, we have(9.25)which is singular across k0. Thus, Re ck is symmetric about k0, while Im ck iseither zero or is antisymmetric about k0.Exercise 9.8Derive the results (9.23) through (9.25) for the discrete Fourier transform of theharmonic oscillator. nNotice that the Fourier transform from x domain at a single frequency produces inthe k domain a dependence that is strongest near the frequency k0. In Section 10.2it is shown that the Fourier integral transform has a pole at k = k0.Formula (9.21) is essentially ready to compute numerically because the real andimaginary parts can be immediately identified in the complex exponential.
The program Discrete Fourier Transform for Oscillator implements formula(9.21), allowing for any N value that is nonzero and including output to a file,DFToscr, that will record the real and imaginary transform values, Re_ck andIm_ck for each k. As coded, the discrete transform is computed for integer valuesof k from 0 to N - 1, but this is easily modified if you want to look at the behaviorof the transform near the discontinuities. When programming the formula the casek = k0 must be handled specially according to (9.24).Exercise 9.9(a) Code and test program Discrete Fourier Transform for Oscillator,adapting the input/output to your computing environment.
Spot check a fewvalues against other calculations.(b) Modify the program, or use the output in file DFToscr, to prepare graphicsof the real and imaginary parts of ck against k, as in Figures 9.4 and 9.5.9.2 DISCRETE FOURIER TRANSFORMS327(c) Compare your results with those in the figure for the indicated values of thefrequencies k0. nPROGRAM 9.1 Discrete Fourier transforms for harmonic oscillators.#include#include<stdio.h><math.h>main(){/* Discrete Fourier Transform for Oscillator */FILE *fout;double pi,kzero;double pi_kzero,sinzero,w,wtop,cfact,Re_ck,Im_ck;int N,k;char wa;pi = 4*atan(l.0);printf("DFT for Oscillator\n");printf("\nWrite over output (w) or Add on (a):\n");scanf("%s",&wa) ;fout = fopen("DFToscr",&wa);N = 2;while(N ! = O){printf("\n\nInput kzero,N (N=O to end): ");scanf("%lf%i",&kzero,&N) ;if ( N == 0 )printf("\nEndDFTforOscillator");exit(O);/* Constants for k loop */pi_kzero = pi*kzero;sinzero = sin(pi_kzero);for ( k = 0; k <= N-l; k++ ) /* integers for k */{if ( k == kzero ){Re_ck = N-0.5; Im_ck = 0;}else{w = (pi_kzero-pi*k)/N;wtop = pi_kzero-w;cfact = sinzero/sin(w);Re_ck = cos(wtop)*cfact-0.5;/* real part */328DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESIm_ck = sin(wtop)*cfact; /* imaginary part */}printf("\n%i %g %g",k,Re_ck,Im_ck);fprintf(fout,"\n%i %g %g", k,Re_ck,Im_ck);}In Figures 9.4 and 9.5 we choose N = 16 and connect the DFT values by linesegments, except across the singularity in Im ck for k0 = 7.5.
The discrete Fouriertransform of the oscillator usually does not have a singularity at k0, as the examplesin the figure show.FIGURE 9.4 Real part of the discrete Fourier transform for harmonic oscillations at frequencyk0, for k0 = 5.2 (dotted), k0 = 10.8 (dash-dotted), k0 = 7.5 (dashed, a constant value). and fork0 = 9 (solid, with a spike at k = 9)..FIGURE 9.5 Imaginary part of the discrete Fourier transform for harmonic oscillations at frequency k0, for k0 = 5.2 (dotted), k0 = 10.8 (dash-dotted), k0 = 7.5 (dashed, discontinuous), andfor k0 = 9 (solid line at zero).9.3 THE FAST FOURIER TRANSFORM ALGORITHM329The symmetry (9.23) is made evident in Figures 9.4 and 9.5 by choosing twok0 values whose sum is N, namely k0 = 5.2 and k0 = 10.8. Notice that for thesek0 values that are not integers or half integers there is a broad distribution of transform values about the oscillator frequency.
These two examples give the extremesof bounded exponential behavior, namely, pure decay and pure oscillation. Themore general case of the damped oscillator is just (9.18), and the transform will havea distribution of values about k0, rather than a discontinuity there if k = k0. The investigation of the damped oscillator DFT is an interesting project that you can easilydevelop from the general expression (9.18) and Program 9.1.9.3THE FAST FOURIER TRANSFORM ALGORITHMDiscrete Fourier transforms, whose theory is developed in Section 9.2, are widelyused in data analysis, in imaging science, and in digital signal processing, because avery efficient algorithm, the fast Fourier transform (FFT), has been developed fortheir evaluation.
The nomenclature of the FFT is misleading; the algorithm appliesto the discrete transform in Section 9.2, not to the Fourier integral transform inChapter 10. The distinction is indicated in Figure 9.1. By various approximations, however, one often estimates the integral transforms by the FFT.Deriving the FFT algorithmHere we derive a straightforward and complete derivation of the FFT algorithm appropriate for real data from which the average value of y (the DC level in electricalparlance) has been subtracted out.
Equation (9.12) may then be written as(9.26)withThen, from (9.6) we have(9.27)Thus, the y values and the c values are treated symmetrically, the only difference intheir treatment being that the exponent is complex-conjugated between their two formulas. Also, the summation range is between 1 and N for both, whereas for arbitrary data the summation range was slightly asymmetric, as shown in Section 9.2.To continue the development, recall from Exercise 2.13 that the complex Nt hroots of unity, En, are given by(9.28)330DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESWe can now write compactly(9.29)Direct computation of the ck requires a time that increases at least as fast as N2because each coefficient requires combining N values of the yj and there are N coefficients.
Let us call such a direct evaluation method a conventional Fourier transform (CFT) algorithm. If the calculation could be divided into two calculations eachof length N/2, the time would be halved, since 2(N/2)2 = N2/2. This divide-andconquer strategy of breaking the calculation of the transform into smaller tasks is thegenesis of various FFT algorithms. We now describe the simplest of these, theradix-2 FFT.Assuming that N is even, one way of dividing the calculation is to combine oddvalues j = 1,3 ,..., N - 1 (that is, j = 2r - 1, with r = 1,..., N/2), then even values j = 2,4,..., N - 2,N (j = 2r, r =, l ,...,N/2).
This pattern is illustrated in Figure 9.6 for N = 8. In this example the complex factor E8 is the eighth root ofunity, namely 8 values spaced by angles ofin the complex plane, as shown.For real yj data the location in the complex plane is appropriate for a DFT calculationof the fundamental frequency, k = 1. The suggested division is that between theouter level, labeled by i = 1, and the next level, i = 2.6FIGURE 9.6 A mandala for conceptualizing the fast Fourier transform for N = 8 data. Thepositions in the complex plane of the 8 points for k = 1 are shown.9.3THE FAST FOURIER TRANSFORM ALGORITHM331With this division the formula (9.29) for the modified Fourier amplitudes Ck canbe written compactly as(9.30)In this equation we have used the distributive property of the complex exponential(9.31)and the exponentiation property(9.32)forNote also that(9.33)Exercise 9.10Prove in detail the immediately preceding properties of the complex exponentialthat are used for developing the FFT algorithm.
(For complex exponentials seeSection 2.3.) nBy this level of division to the inner sum in (9.30), there are now two discreteFourier transforms, each of half the previous length.The process of halving each transform can be continued as long as the number ofpoints to be divided in each group at each level is divisible by 2. For example, inFigure 9.6 at level i = 3 there are 4 groups, with only two points in each.
Thesepoints are just combined with a sign difference, since they are exactly out of phasewith each other. Note that within each group along an arc the relative phase betweenpoints is the same as at the outer level, i = 1.For simplicity, we restrict ourselves to cases in which the total number of datapoints, N, is a power of 2:(9.34)where v is a positive integer. This is known as a radix-2 FFT.
In our worked example in Figure 9.6 we have v = 3. Then the DFT can be halved v times, endingwith N/2 subintervals, each of length 2, as in our example.This essentially completes our outline of the radix-2 fast Fourier transform. After a little bookkeeping we will be ready to develop a computer program for the FFT,as we do in Section 9.7.332DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESBit reversal to reorder the FFT coefficientsAt each step of the fast Fourier transform the values of the coefficients Ck becomereordered because of the splitting of the transform into odd and even terms.
This reordering is called bit reversal. We will now see how bit reversal can be efficiently3done for our radix-2 FFT. Consider the example in Figure 9.6 for N = 2 = 8.The levels involved are i = 1,2 ,..., v, with v = 3. At the innermost level in thisexample, moving inward on the semicircular arcs, the subscript of the yj values arej = 1,5,3,7,2,6,4,8. What’s the pattern here?When we sort a sequence into odd and even terms, the odd terms all have thesame least-significant bit in a binary representation, while the even terms all have thesame (other) least-significant bit.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
