# 5: Fourier Interpolant and DFT 1 The Discrete Fourier Transform (DFT) of a periodic array fj , for j = 0, 1, . . . , N 1 (corresponding to data at equally spaced points, starting at the left end point of the interval of periodicity) is evaluated via the Fast Fourier Transform (FFT) algorithm (N power of 2). Use an FFT package, i.e. an already coded FFT (the functions t and it in Matlab or numpy.t in python). 1. Which of the following expressions dene the Fourier coecients (the DFT) that your t package (name it) returns for k = 0, 1, . . . , N 1? N 1 fj ei2kj/N (1) fj ei2kj/N (2) fj ei2k(j1)/N ck = (3) j=0 N ck = j=1 N ck = j=0 2. Which of the following expressions dene the inverse DFT computed by t package? fj 2 = N fj = fj = 1 N 1 N N 1 ck ei2kj/N for j = 0, 1, . . . , N 1 (4) ck ei2kj/N for j = 0, 1, . . . , N 1 (5) k=0 N 1 k=0 N ck ei2kj/N for j = 0, 1, . . . , N 1 (6) k=1 3. Prove that if the yj , for j = 0, 1, . . . , N 1 are real numbers then c0 is real and N k cN k = ck , where the bar denotes complex conjugate. Hint: prove that N k = N . 4. Let PN (x) be the trigonometric polynomial of lowest order that interpolates the periodic array fj , j = 0, 1, . . . , N 1 at the equidistributed nodes xj = j(2/N ), 1 All course materials (class lectures and discussions, handouts, homework assignments, examinations, web materials) and the intellectual content of the course itself are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of the instructor. 1 j = 0, 1, . . . , N 1, i.e PN (x) = a0 + 2 N/21 (ak cos kx + bk sin kx) + k=1 aN/2 cos 2 N kx 2 (7) for x [0, 2], where ak 2 = N bk = 2 N N 1 fj cos kxj for k = 0, 1, . . . , N/2, (8) fj sin kxj for k = 1, . . . , N/2 1. (9) j=0 N 1 j=0 Write a formula that relates the complex Fourier coecients computed by your t package to the real Fourier coecients, ak and bk , that dene PN (x). 5. Using your t package and Prob. 4, nd P8 (x) on [0, 2] for the following periodic array: f0 = 6.000000000000000 f1 = 10.242640687119284 f2 = 2.000000000000000 f3 = 2.585786437626905 f4 = 2.000000000000000 f5 = 1.757359312880716 f6 = 6.000000000000000 f7 = 5.414213562373098 6. Let fj = esin xj , xj = j2/N for j = 0, 1, . . . , N 1. Take N = 8. Using your t package obtain P8 (x) and nd a spectral approximation of the derivative of esin x at xj for j = 0, 1, . . . , N 1 by computing P8 (xj ). Compute the actual error in the approximation. 7. Note that PN (x) is again a trigonometric polynomial of degree N/2 and whose coecients can be computed from those of PN (x) via the FFT. (a) Write a code to compute a (spectral) approximation to the derivative at xj = j(2/N ) for j = 0, 1, . . . , N 1 for the corresponding periodic array fj , j = 0, 1, . . . , N 1, using one DFT and one inverse DFT, i.e. in order N log2 N operations. (b) Test your code by comparing with your answer of the previous problem. (c) Observe the behavior of the error as N is increased to 16, and 32. Note: the contribution to the derivative from the k = N/2 node should be zero. Make sure to set the Fourier coecient for k = N/2 of the derivative equal to zero. 2 Approximation of Periodic Functions Hector D. Ceniceros 1 Approximation by Trigonometric Polynomials We now consider the problem of approximating a periodic, real-valued function f . Without loss of generality we can assume that f is of period 2 (if it p is of period p then the function F (y) = f ( 2 y) has period 2). We are going to assume that f is piece-wise continuous with finite jumps, i.e. the right and left limits (1) (2) lim f (x + h) = f ( x+ ), h0+ lim f (x h) = f ( x ) h0 exist and are finite. We would like to approximate such a periodic function f by a Trigononometric polynomial: n X 1 Sn (x) = a0 + (3) (ak cos kx + bk sin kx). 2 k=1 We can also write this in complex form, by noting eikx = cos kx + i sin kx, we have n X (4) Sn (x) = ck eikx k=n These are lecture notes for Math 104 A. These notes and all course materials are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of the instructor. 1 where 1 c 0 = a0 2 1 ck = (ak ibk ), 2 1 ck = (ak + ibk ), 2 (5) (6) (7) k = 1, . . . , n k = 1, . . . , n. Note that ck = ck , where the bar denotes the complex conjugate. We want to approximate f by Sn in the Least Squares sense 2 kf Sn k2 = (8) \u0013 21 2 \u0012Z [f (x) Sn (x)] dx = min 0 Let Jn = kf Sn k22 . Then we need to find the coefficients ck , k = 0, 1, 2, . . . n (or equivalentlya0 , a1 , ..., an , b1 , ..., bn ) which minimize Jn . We have #2 Z 2 " n X f (x) ck eikx dx Jn = 0 (9) k=n [f (x)]dx 2 = 0 + n X 2 Z Z k=n n X n X Z 2 ck f (x)eikx dx 0 2 eikx eilx dx. ck cl 0 k=n l=n This problem simplifies if we use the orthogonality of the set {1, eix , eix , . . . , einx , einx } or equivalently of the set of trigonometric functions {1, cos x, cos 2x, . . . , cos nx, sin x, sin 2x, . . . , sin nx}. As for k 6= l Z 2 Z ikx ilx e e dx = (10) 0 2 e i(k+l)x 0 2 1 i(k+l)x dx = e =0 i(k + l) 0 and for k = l Z (11) 2 e Z ikx ilx e dx = 0 dx = 2 0 2 2 Then we get Z Jn = (12) n X 2 [f (x)]dx 2 0 2 Z ikx f (x)e ck dx + 2 0 k=n n X ck ck . k=n To find the minimum we determine the critical point of Jn as a function of the ck , k = 0, 1, . . . , , n. Z 2 Jn (13) f (x)dx + 2(2)c0 = 0 = 2 c0 0 Z 2 Jn (14) f (x)eilx dx + 2(2)cl = 0, l = 1, . . . , n. = 2 cl 0 Therefore, relabeling the coefficients with k again, we get Z 2 1 ck = (15) f (x)eikx dx, k = 0, 1, . . . , n, 2 0 which are the complex Fourier coefficients of f . We can now obtain the real Fourier coefficients ak and bk by noting that (16) (17) (18) a0 = 2c0 ak = ck + ck , k = 1, . . . , n bk = i(ck ck ), k = 1, . . . , n. Hence (19) (20) Z 1 2 ak = f (x) cos kxdx, 0 Z 1 2 bk = f (x) sin kxdx, 0 k = 0, 1, . . . , n, k = 1, . . . , n. Now, if we substitute the Fourier coefficients (15) in (12) we get Z 2 n X 0 Jn = [f (x)]2 dx 2 |ck |2 0 k=n that is (21) n X 1 |ck | 2 k=n 2 3 Z 0 2 [f (x)]2 dx. This is known as Bessel's inequality. If we substitute (5)-(7) we obtain Bessel's inequality for the real Fourier coefficients (22) Z n 1 2 X 2 1 2 2 a + [f (x)]2 dx. (a + bk ) 2 0 k=1 k 0 R 2 If 0 [f (x)]2 dx is finite, as is the case for the f under consideration, then then series 1 2 X 2 a + (a + b2k ) 2 0 k=1 k converges and consequently limk ak = limk bk = 0. The convergence of a Fourier series is a delicate question. For a continuous function f , it does not follow that its Fourier series converges point-wise to it, only that it does so in the mean Z 2 (23) lim [f (x) Sn (x)]2 dx = 0. n 0 This convergence in the mean for a continuous periodic function implies that Bessel's inequality becomes the equality (24) Z 1 2 X 2 1 2 2 [f (x)]2 dx, a + (a + bk ) = 2 0 k=1 k 0 which is known as Parseval's identity. We state now a convergence result without proof. Theorem 1. Suppose that f is piecewise continuous and periodic in [0, 2] and with a piecewise continuous first derivative. Then X \u0003 1 1\u0002 a0 + (ak cos kx + bk sin kx) = f (x+ ) + f (x ) 2 2 k=1 for each x [0, 2] where ak , k = 0, 1, . . . , n, and bk k = 1, 2, . . . , n are the Fourier coefficients of f . In particular if x is a point of continuity of f X 1 a0 + (ak cos kx + bk sin kx) = f (x). 2 k=1 4 We have been working on the interval [0, 2] but we can choose any other interval of length 2. This is so because if g is 2 periodic then we have the following result Lemma 1. Z 2 (25) Z t+2 g(x)dx = 0 g(x)dx t for any real t. Proof. Define Z t+2 g(x)dx G(t) = t Then Z 0 G(t) = Z g(x)dx + t t+2 Z t+2 g(x)dx g(x)dx = 0 Z 0 t g(x)dx. 0 Then, by the Fundamental theorem of calculus G0 (t) = g(t + 2) g(t) = 0 since g is 2 periodic. Thus, G is independent of t. Example 1. f (x) = |x| on [, ]. ak = = = = Z Z 1 1 f (x) cos kxdx = |x| cos kxdx Z 2 x cos kxdx 0 u = x, dv = cos kxdx 1 du = dx, v = sin kx k \u0014 \u0015 1 Z 2 2 x sin kx sin kxdx = cos kx 2 k k 0 k 0 0 2 [(1)k 1], k = 1, ... k 2 5 and a0 = , bk = 0 for all k as the function is even. Therefore S(x) = lim Sn (x) n X 1 2 = + [(1)k 1] cos kx 2 2 k k=1 \u0014 \u0015 4 cos x cos 3x cos 5x 1 = + + + . 2 12 32 52 How do we find accurate numerical approximations to the Fourier coefficients? We know that for periodic smooth integrands, integrated over one (or multiple) period(s), the Composite Trapezoidal Rule using equidistributed nodes gives super-algebraic accuracy. Let xj = jh, j = 0, 1, . . . N , h = 2/N , then we can approximate " # N 1 X 1 1 1 ikx0 ikxN ikxj f (x0 )e . + f (xj )e + f (xN )e ck h 2 2 2 j=1 But due to periodicity f (x0 )eikx0 = f (xN )eikxN so we have (26) N N 1 1 X 1 X ikxj ck f (xj )e = f (xj )eikxj . N j=1 N j=0 and for the real Fourier coefficients we have (27) N N 1 2 X 2 X ak f (xj ) cos kxj = f (xj ) cos kxj , N j=1 N j=0 (28) N N 1 2 X 2 X bk f (xj ) sin kxj = f (xj ) sin kxj . N j=1 N j=0 2 Interpolating Fourier Polynomial Let f be a 2-periodic function and xj = j 2 , j = 0, 1, . . . N , equidistributed N nodes in [0, 2]. The interpolation problem is to find a trigonometric polynomial of lowest order Sn such that Sn (xj ) = f (xj ), for j = 0, 1, . . . , N . Because of periodicity f (x0 ) = f (xN ) so we only have N independent values. 6 If we take N = 2n then Sn has 2n + 1 = N + 1 coefficients. So we need one more condition. At the equidistributed nodes, the sine term of highest frequency vanishes as sin( N2 xj ) = sin(j) = 0 so the coefficient bN/2 is irrelevant for interpolation. We thus look for an interpolating trigonometric polynomial of the form \u0012 \u0013 N/21 X 1 1 N (29) PN (x) = a0 + (ak cos kx + bk sin kx) + aN/2 cos x . 2 2 2 k=1 The convenience of the 1/2 factor in the last term will be seen in the formulas for the coefficients below. It is conceptually and computationally simpler to work with the corresponding polynomial in complex form (30) N/2 X 0 PN (x) = ck eikx , k=N/2 where the prime in the sum means that the first and last term (for k = N/2 and k = N/2) have a factor of 1/2 and cN/2 = cN/2 , equivalent to the bN/2 = 0 condition in (29). Theorem 2. PN (x) = N/2 X 0 ck eikx k=N/2 with N 1 1 X f (xj )eikxj , ck = N j=0 k= N N ,..., 2 2 interpolates f at the equidistributed points xj = j(2/N ), j = 0, 1, . . . , N . Proof. We have PN (x) = N/2 X 0 k=N/2 ck e ikx = N 1 X j=0 N/2 1 X0 ik(xxj ) f (xj ) e . N k=N/2 Defining (31) N/2 1 X0 ik(xxj ) e lj (x) = N k=N/2 7 we have PN (x) = (32) N 1 X f (xj )lj (x).) j=0 Thus, we only need to prove that for j and m in the range 0, . . . , N 1 ( 1 for m = j, lj (xm ) = (33) 0 for m = 6 j. N/2 1 X0 ik(mj)2/N e . lj (xm ) = N k=N/2 But ei(N/2)(mj)2/N = ei(mj) = (1)(mj) so we can combine the first and the last term and remove the prime from the sum: 1 lj (xm ) = N 1 = N N/21 X k=N/2 N/21 X ei(k+N/2)(mj)2/N ei(N/2)(mj)2/N k=N/2 i(mj) =e eik(mj)2/N N 1 1 X ik(mj)2/N e N k=0 and, as we proved in the introduction ( N 1 X 0 if j m is not divisible by N 1 (34) eik(jm)2/N = N k=0 1 otherwise. Then (33) follows and PN (xm ) = f (xm ), m = 0, 1, . . . N 1. 8 Using the relations (16)-(18) between the ck and the ak and bk coefficients we find that \u0012 \u0013 N/21 X 1 N 1 (ak cos kx + bk sin kx) + aN/2 cos x PN (x) = a0 + 2 2 2 k=1 interpolates f at the equidistributed nodes xj = j(2/N ), j = 0, 1, . . . , N if and only if (35) N 2 X f (xj ) cos kxj , ak = N j=1 (36) N 2 X bk = f (xj ) sin kxj . N j=1 A smooth periodic function f can be approximated very accurately by its Fourier interpolant PN . Note that derivatives of PN can be easily computed (p) PN (x) (37) N/2 X 0 = (ik)p ck eikx k=N/2 The discrete Fourier coefficients of the p-th derivative of PN are (ik)p ck . Thus, once these Fourier coefficients have been computed a very accurate (p) approximation of the derivatives of f is obtained, f (p) (x) PN (x). Figure 1 shows the approximation of f (x) = sin xecos x on [0, 2] by P8 . The graph of f and of the Fourier interpolant are almost indistinguishable. Let us go back to the complex Fourier interpolant (30). Its coefficients ck are periodic of period N , (38) ck+N N 1 N 1 1 X 1 X i(k+N )xj fj e = fj eikxj eij2 = ck = N j=0 N j=0 and in particular cN/2 = cN/2 . Using the interpolation property and setting fj = f (xj ), we have (39) fj = N/2 X 0 N/21 ikxj ck e k=N/2 = X k=N/2 9 ck eikxj 1.5 1 0.5 0 0.5 1 1.5 0 1 2 3 4 5 6 Figure 1: S8 (x) for f (x) = sin xecos x on [0, 2] and N/21 X (40) ck e ikxj = k=N/2 N/21 1 X ikxj ck e + = ck eikxj k=0 k=N/2 N 1 X X N/21 ck e ikxj + X k=0 k=N/2 ikxj ck e = N 1 X ck eikxj , k=0 where we have used that ck+N = ck . Combining this with the formula for the ck 's we get Discrete Fourier Transform (DFT) pair (41) (42) N 1 1 X ck = fj eikxj , N j=0 fj = N 1 X ck eikxj , k = 0, . . . , N 1, j = 0, . . . , N 1. k=0 The set of discrete coefficients (41) is known as DFT of the periodic array f0 , f1 , . . . , fN 1 and (42) is referred to as the Inverse DFT. The direct evaluation of the DFT is computationally expensive, it requires order N 2 operations. However, there is a remarkable algorithm which 10 achieves this in merely order N log N operations. This is known as the Fast Fourier Transform. 11 The Fast Fourier Transform Hector D. Ceniceros 1 The Fast Fourier Transform The DFT was defined as (1) N 1 1 X ck = fj eikxj , N j=0 (2) fj = N 1 X ck eikxj , k = 0, . . . , N 1, j = 0, . . . , N 1. k=0 The direct computation of either (1) or (2) is requires order N 2 operations. As N increasing the cost quickly becomes prohibitive. In many applications N could easily be on the order of thousands, millions, etc. One of the top algorithms of all times is the Fast Fourier Transform (FFT). It is usually attributed to Cooley and Tukey (1965) but its origin can be tracked back to C. F. Gauss (1777-1855). We now look at the main ideas of this famous and widely used algorithm. Let us define dk = N ck for k = 0, 1, . . . , N 1. Then we can rewrite (1) as (3) dk = N 1 X kj fj N , j=0 These are lecture notes for Math 104 A. These notes and all course materials are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of the instructor. 1 where N = ei2/N . In matrix 0 N d0 0 d1 N (4) .. = .. . . 0 dN 1 N form 0 N 1 N .. . .. . N 1 N 0 N N 1 N .. . f0 f1 .. . (N 1)2 fN 1 N Let us call the matrix on the right FN . Then FN is N times the matrix of the DFT and the matrix of the Inverse DFT is simply the complex conjugate of FN . This follows from the identities: N 1 2 1 + N + N + + N = 0, 2(N 1) 2 4 1 + N + N + + N = 0, .. . 2(N 1) N 1 1 + N + N (N 1)2 + + N N (N 1) N 2N 1 + N + N + + N = 0, = N. We already proved the first of these identities. This geometric sum is equal N N to (1 N )/(1 N ), which in turn is equal to zero because N = 1. The other are proved similarly. We can summarize these identities as ( N 1 X 0 j 6= 0 (mod N ) ik (5) N = N j = 0 (mod N ), k=0 where j = k (mod N ) means j k is an integer multiple of N . Then ( N 1 N 1 X X 0 j= 6 l (mod N ) 1 1 1 (jl)k jk lk (6) (FN FN )jl = N = N N = N N k=0 N k=0 1 j = l (mod N ). which shows that FN is the inverse of N1 FN . Let N = 2n. Going back to (3), if we split the even-numbered and the odd-numbered points we have (7) dk = n1 X f2j k2jk + n1 X j=0 j=0 2 (2j+1)k f2j+1 N But (9) ijk 2 N 2 2jk N = ei2jk N = e (8) (2j+1)k N 2 2 2 = eijk n = nkj , 2 2 k kj n . = ei(2j+1)k N = eik N ei2jk N = N So denoting fje = f2j and fj0 = f2j+1 , we get (10) dk = n1 X fje njk + j=0 k N n1 X fj0 njk j=0 We have reduced the problem to two DFT of size n = N2 plus N multiplicak , k = 0, 1, . . . , N 1 only depend on tions (and N sums). The numbers N N so they can be precomputed once and store for other DFT of the same size N . If N = 2p , for p positive integer, we can repeat the process to reduce each of the DFT's of size n to a pair of DFT's of size n/2 plus n multiplications (and n additions), etc. We can do this p times so that we end up with 1-point DFT's, which require no multiplications! Let us count the number of operations in the FFT algorithm. For simplicity, let is count only the number of multiplications (the numbers of additions is of the same order). Let mN be the number of multiplications to compute the DFT for a periodic array of size N and assume that N = 2p . Then mN = 2m N + N 2 = 2m2p1 + 2p = 2(2m2p2 + 2p1 ) + 2p = 22 m2p2 + 2 2p = = 2p m20 + p 2p = p 2p = N log2 N, where we have used that m20 = m1 = 0 ( no multiplication is needed for DFT of 1 point). To illustrate the savings, if N = 220 , with the FFT we can obtain the DFT (or the Inverse DFT) in order 20 220 operations, whereas 1 20 the direct methods requires order 240 , i.e. a factor of 20 2 52429 more operations. 3 The FFT can also be implemented efficiently when N is the product of small primes. A very efficient implementation of the FFT is the FFTW (\"the Fastest Fourier Transform in the West\"), which employs a variety of code generation and runtime optimization techniques and is a free software. 4