Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 48
Текст из файла (страница 48)
We discuss thislast question in the next section.In this connection, we point out that (6.49) implies(6.61)in case f(x) has k derivatives, with the kth derivative piecewise continuous.EXERCISES6.5-l Calculate the Fourier series for the 2π-periodic function f(x) given by f(x) = x on[0,2π).6.52 Verify that the 2π-periodic function f(x) whose values on [0, 2π) are given by*6.6FAST FOURIER TRANSFORMS277is continuous and has a continuous first derivative (as a 2 π−periodic function), but has jumpsin the second derivative.
Then construct the spectrum of f(x) and show that it decays like j-3(and no faster) as6.5-3 Write the Fourier series obtained in Exercise 6.5-2 in terms of sines and cosines. Whywould you expect all the aj ’s to be zero?6.5-4 If f(x) is a 2π-periodic function, then so is the function gm(x) = f(mx), for any integerm. What is the relationship between the f(j) and thefor any6.5-5 If f(x) is a 2π-periodic function, then so is the functionnumber α. What is the relationship between the f(j) and the6.5-6 Suppose that f(x) is a very smooth function of period τ.
But, in converting it to a2 π−periodic function g(x) = f(τ x/(2 π )), you mistakenly use τ´ instead of τ for someWhat is the likely effect of this mistake on the computed Fourier coefficients6.5-7 Prove that, if f(x) is a real function, thenfor an appropriate phase shift θj. (Hint: Use the fact that any complex number z can bewritten in polar form as |z|eiθ for an appropriate θ.)6.5-8 Prove (6.53). (Hint: Recall how to sum a geometric series.)6.5-9 Prove Theorem 6.5.6.5-10 Derive the addition formulas for sin(6.42) and from the law of exponents: eA+B6.5-11 Prove that if N is even and f(x) is real, thenand cos (α + β) from Euler’s formulainterpolates f(x) at the sampling points xk, all k.6.5-12 How would you construct the trigonometric interpolant to f(x) at the points α +k2 π /N, k = 0 , . .
. , N- 1, with a some positive number less than 2π /N ?*6.6 FAST FOURIER TRANSFORMSIn discrete harmonic, or Fourier, analysis, one calculates the numbers(6.62)withall k(6.63)x k = 2πk / Nin order to resolve the 2π-periodic motion described by f(x) into simpleharmonics. As we saw in Sec. 6.5,withif f(x) has a piecewise continuous kth derivative. One is interested in whichfrequencies are present in f(x) and in their strength. But, because of thealiasing effect,is useless as an approximation tofor |j| > N/2,and is usually a good approximation only for |j| much smaller than N/2.This makes it desirable to calculatefor “large” N, and so brings upthe important question of just how one is to calculateefficiently.278APPROXIMATIONIt is clear that the evaluation of any particularrequiresmultiplications and additions.
The straightforward calculation of N suchfor |j| < N/2) would therefore takenumbers (e.g., the numbersoperations. Thus, already for 1000 sample points, we would needmillions of operations, and, until recently, this was a major obstacle to theuse of discrete Fourier analysis.This situation changed dramatically when it became well known thatthe simultaneous calculation of N consecutiveneed only take(N log N) arithmetic operations because of the strong interrelationsbetween these numbers.
The key word for this has been fast Fouriertransform, or FFT, and it has made calculations with N < 1000 routine; ithas even made it possible to use N’s in the tens of thousands.We are here able only to give an indication of the basic ideas whichhave led to such a dramatic increase in efficiency. The latest word in 1978on these matters is to be found in a paper by S. Winograd [36]. Inparticular, work done before and after publication of Cooley and Tukey’sseminal article [37] has long made clear that there are many FFTs and that,for greatest efficiency, it is necessary (and profitable) to write a differentprogram for each different value of N one wishes to use.For the analysis of the computations of the numbersfor |j| <N/2 from the numbers f(x0), . .
. , f(xN-1), it is convenient to introduce thediscrete Fourier transform FN, which carries the N-vectorto the N-vectorgiven by(6.64)with ωN an Nth root of unity,The connection between the calculation ofand this discrete Fouriertransform is as follows. If we take the particular N-vectorthen|j| < N/2(6.65)Thus, it is sufficient to concentrate on the efficient calculation of thediscrete Fourier transform.We begin this discussion with the observation that as given by (6.64)is a polynomial of degree < N in the quantityhence can beevaluated in N operations, by nested multiplication. Here, we count oneaddition plus one multiplication as one operation. It would therefore takeN2 operations for the straightforward evaluation of (6.64) for all j.*6.6FAST FOURIER TRANSFORMS279The most widely known idea for an FFT has been popularized byCooley and Tukey.
It is applicable whenever N is a product of integers. Wenow discuss this idea first in the case thatN = P · QThink of the N-vector z as stored FORTRAN-fashion in a one-dimensional array. Then we can interpret the array also FORTRAN-fashion as atwo-dimensional array Z, of dimension (P,Q). This means thatCorrespondingly, we factor the suminto a double sum,Here, we have made use of the fact thatThis makes apparent thecrucial fact that the inner sum in the last right hand side is Q-periodic in v,i.e., replacing v by v + Q does not change its value, due to the fact that1.
This means that we need only calculate this sum for v = 1, . . . , Q(and each p). Thus, for each p = 1, . . . , P, we calculate from the Q -vectorZ ( p, ·) the Q-vector whose entries are the numbersi.e., we calculate the discrete Fourier transform of the Q-vectors Z( p , ·) ,p = 1, . . . , P, at a total cost of P · Q2 = N · Q operations.Now, we could store the transform of Z( p, ·) over Z( p, ·).
But, inanticipation of further developments, we choose to store the transform ofZ (p, ·) in Z1 (·, p), where Z1 is a two-dimensional array of size (Q, P),rather than (P, Q).With this, our calculation of is reduced to the evaluation of the sumHere, we have used the notation vQ to indicate the integer between 1 and Qfor which v - vQ is divisible by Q. Thus,v = vQ + Q(v´ - 1)for some integer v´ between 1 and P.
In effect,(6.66)if we interpret the vectorFORTRAN-fashion as a two-dimensional array280APPROXIMATIONZ0 of size (Q, P). With this, we must calculateHere, the right-hand side is a polynomial of degree < P in the quantityThis quantity can be generated step by step, as inthe following convenient arrangement of the calculations.(6.67)The sum in the innermost loop is, of course, to be evaluated by nestedmultiplication. The total cost of this step is then Q · P2 = NP operations(if we neglect the N multiplications needed to generate the various x’s). Inthis way, we have obtained in Z0 the discrete Fourier transform of z at acost of only N(P + Q) operations compared to the N2 operations requiredfor the naive way.If now N is the product of three or more integers greater than 1,N = P1 · · · Pmsay, then we can calculate the discrete Fourier transform of z even morecheaply, by using the second step (6.67) in a slightly more sophisticatedway.For the description, we need a bit of notation to indicate how a givenone-dimensional array is interpreted FORTRAN-fashion equivalently as atwo- or a three-dimensional array.
If Z is a one-dimensional array oflength N, then we denote by ZA the equivalent two-dimensional array ofdimension (A, N/A), and by ZA,B the equivalent three-dimensional arrayof dimension (A, B, N/(A B)). In this way,Z A , B(a, b, c) = ZA (a, b + B(c - 1)) = ZAB ( a + A( b - 1), c)=Z ( a + A(b - 1 + B(c - 1)))Let now Z be a one-dimensional array containing z, as before, and fork = 0, . . . , m, let Zk be a one-dimensional array containing the discreteFourier transform of sections of Z as follows:(6.68)with(6.68a)*6.6FAST FOURIER TRANSFORMS281Note that Z fits the role of Zm and that Z0 contains = FNz. To get fromZk to Zk-1, use the following slightly extended version of (6.67), withB, P, A as given in (6.68a):(6.69)Indeed, the algorithm producesOn the other hand, (6.68) implies thatTherefore,But now, since1, we may add to the exponent on the right handside any integer multiple of AP, and this allows the conclusion thatand so proves that Zk-1, as produced by (6.69), satisfies (6.68) (with kreplaced by k - 1).In particular, Z 0 contains the discrete Fourier transform of z.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
