Unformatted text preview: CHAPTER 2 TREFETHEN 1994 80 Chapter 2. Fourier Analysis
2.1. The Fourier transform 2.2. The semidiscrete Fourier transform 2.3. Interpolation and sinc functions 2.4. The discrete Fourier transform 2.5. Vectors and multiple space dimensions 2.6. Notes and references ::: Untwisting all the chains that tie the hidden soul of harmony.  JOHN MILTON, L'Allegro (1631) CHAPTER 2 TREFETHEN 1994 81 The last chapter dealt with time dependence, and this one is motivated by space dependence. Later chapters will combine the two. Fourier analysis touches almost every aspect of partial di erential equations and their numerical solution. Sometimes Fourier ideas enter into the analysis of a numerical algorithm derived from other principlesespecially in the stability analysis of nitedi erence formulas. Sometimes they underlie the design of the algorithm itselfspectral methods. And sometimes the situation is a mixture of both, as with iterative and multigrid methods for elliptic equations. For one reason or another, Fourier analysis will appear in all of the remaining chapters of this book. The impact of Fourier analysis is also felt in many elds besides di erential equations and their numerical solution, such as quantum mechanics, crystallography, signal processing, statistics, and information theory. There are four varieties of Fourier transform, depending on whether the spatial domain is unbounded or bounded, continuous or discrete: Name Space variable Transform variable Fourier transform unbounded, continuous continuous, unbounded Fourier series bounded, continuous discrete, unbounded semidiscrete Fourier transform unbounded, discrete continuous, bounded or ztransform discrete Fourier transform bounded, discrete discrete, bounded (DFT) (The second and third varieties are mathematically equivalent.) This chapter will describe the essentials of these operations, emphasizing the parallels between them. In discrete methods for partial di erential equations, one looks for a representation that will converge to a solution of the continuous problem as the mesh is re ned. Our de nitions are chosen so that the same kind of convergence holds also for the transforms. Rigorous Fourier analysis is a highly technical and highly developed area of mathematics, which depends heavily on the theory of Lebesgue measure and integration. We shall make use of L2 and `2 spaces, but for the most part this chapter avoids the technicalities. In particular, a number of statements made in this chapter hold not at every point of a domain, but \almost everywhere"  everywhere but on a set of measure zero. 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 82 2.1. The Fourier transform
kuk =
hZ 1 If u(x) is a (Lebesguemeasurable) function of x 2 R , the L2 norm of u is the nonnegative or in nite real number The symbol L2 (\Ltwo") denotes the set of all functions for which this integral is nite: L2 = fu : kuk < 1g: (2:1:2) Similarly, L1 and L1 are the sets of functions having nite L1  and L1 norms, de ned by
kuk1 =
Z1 ;1 ju(x)j2 dx i1=2 : (2:1:1) ;1 ju(x)jdx kuk1 = ;1sup 1 ju(x)j: <x< (2:1:3) Note that since the L2 norm is the norm used in most applications, because of its many desirable properties, we have reserved the symbol k k without a subscript for it. The convolution of two functions u v is the function u v de ned by (u v)(x) = (u v)(x) =
Z1 assuming these integrals exist. One way to think of u v is as a weighted moving average of values u(x) with weights de ned by v(x), or vice versa. For any u 2 L2, the Fourier transform of u is the function u( ) de ned ^ by
Z1 ;1 u(x ; y)v(y)dy = Z1 ;1 u(y)v(x ; y)dy (2:1:4) *If u 2 L1 , then u( ) exists for every and is continuous with respect to . According to the ^ RiemannLebesgue Lemma, it also satis es ju( )j! 0 as !1. ^ ;i x u( ) = (F u)( ) = ^ 2 R: ;1 e u(x)dx The quantity is known as the wave number, the spatial analog of frequency. For many functions u 2 L2 , this integral converges in the usual sense for all 2 R , but there are situations where this is not true, and in these cases one 2 must R M interpret the integral as a limit in a certain L norm sense of integrals ;M as M ! 1. The reader interested in such details should consult the various books listed in the references.* 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 83 x Figure 2.1.1. Space and wave number domains for the Fourier transform (compare Figures 2.2.1 and 2.4.1).
The following theorem summarizes some of the fundamental properties of Fourier transforms.
THE FOURIER TRANSFORM Theorem 2.1.* If u 2 L2, then the Fourier transform Z1 ;i x u( ) = (F u)( ) = ^ 2R (2:1:5) ;1 e u(x)dx belongs to L2 also, and u can be recovered from u by the inverse Fourier ^ transform ;1u)(x) = 1 Z 1 ei xu( )d u(x) = (F ^ x 2 R: (2:1:6) 2 ;1 ^ The L2 norms of u and u are related by Parseval's equality, ^
p kuk = 2 kuk: ^ (2:1:7) (2:1:8) If u 2 L2 and v 2 L1 (or vice versa), then u v 2 L2, and udv satis es udv( ) = u( )^( ): ^ v These four equations are of such fundamental importance that they are worth commenting on individually, although it is assumed the reader has already been exposed to Fourier analysis.
*As mentioned in the introduction to this chapter, some of these propertiesnamely equations (2.1.6) and (2.1.8)hold merely for \almost every" value of x or . In fact even if f (z) is a continuous function in L2 , its Fourier transform may fail to converge at certain points x. To ensure pointwise convergence one needs additional assumptions such as that f is of bounded variation (de ned below before Theorem 2.4) and belongs to L1 . These assumptions also ensure that at any point x where f has a jump discontinuity, its Fourier transform converges to the average value (f (x; )+ f (x+))=2. 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 84 First of all, (2.1.5) indicates that u( ) is a measure of the correlation of ^ u(x) with the function ei x. The idea behind Fourier analysis is to interpret u(x) as a superposition of monochromatic waves ei x with various wave numbers , and u( ) represents the complex amplitude (more precisely: amplitude ^ density with respect to ) of the component of u at wave number . Conversely, (2.1.6) expresses the synthesis of u(x) as a superposition of its components ei x, each multiplied by the appropriate factor u( ). The factor ^ 2 is a nuisance that could have been put in various places in our formulas, but is hard to eliminate entirely. Equation (2.1.7), Parseval's equality, is a statement of energy conservation: the L2 energy of any signal u(x) is equal p the sum of the energies of to its component vibrations (except for the factor 2 ). By \energy" we mean the square of the L2 norm. Finally, the convolution equation (2.1.8) is perhaps the most subtle of the four. The left side, udv( ), represents the strength of the wave number component that results when u is convolved with vin other words, the degree to which u and v beat in and out of phase with each other at wave number when multiplied together in reverse order with a varying o set. Such beating is caused by a quadratic interaction of the wave number component in u with the same component of vhence the righthand side u( )^( ). ^ v All of the assertions of Theorem 2.1 can be veri ed in the following example, which the reader should study carefully.
EXAMPLE 2.1.1. Bsplines. Suppose u is the function 1 u(x) = 2 for 1 x 1,
; p 0 otherwise (Figure 2.1.2). Then by (2.1.1) we have u = 1= 2, and (2.1.5) gives (2:1:9) u( ) = ^ 1 2 Z1 k k (This function u( ) is called a sinc function more on these in 2.3.) From (2.1.1) and the ^ indispensable identity* Z 1 sin2 s ds = (2:1:11)
x ;i x 1 = sin : e;i x dx = e 2i ;1 ;1
; (2:1:10) which can be derived by complex contour integration, we calculate u = , which con rms ^ (2.1.7). From the de nition (2.1.4) it is readily veri ed that in this example 1 (1 x =2) for 2 x 2, (u u)(x) = 2 (2:1:12) 0 otherwise
k k p ;j j ; ;1 s2 *worth memorizing! 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 85 u u ^
0 3 ;3
u u x ;15
d u u 0 15 ;3
u u u 0 3 x ;15
udu u 0 15 ;3
forms. 0 3 x ;15 0 15 Figure 2.1.2. The rst three Bsplines of Example 2.1.1 and their Fourier transand 0 otherwise, and by (2.1.8) and (2.1.10), the corresponding Fourier transforms must be
2 d u u ( ) = sin2 3 u d u ( ) = sin3 : u 8 3 1 x2 for 1 x 1, >4 4 < 1 (9 6 x + x2 ) for 1 (u u u)(x) = > 8 x 3, :
; ; ; j j j j (2:1:13) (2:1:14) See Figure 2.1.2. In general, a convolution u(p) of p copies of u has the Fourier transform u(p) ( ) = d Ff u u p u ( ) = sin :
g (2:1:15) Note that whenever u(p) or any other function is convolved with the function u of (2.1.9), it becomes smoother, since the convolution amounts to a local moving average. In particular, u itself is piecewise continuous, u u is continuous and has a piecewise continuous rst derivative, u u u has a continuous derivative and a piecewise continuous second derivative, and so on. In general u(p) is a piecewise polynomial of degree p 1 with a continuous (p 2)nd derivative and a piecewise continuous (p 1)st derivative, and is known as a Bspline. (See, for example, C. de Boor, A Practical Guide to Splines, Springer, 1978.) Thus convolution with u makes a function smoother, while the e ect on the Fourier transform is to multiply it by sin = and thereby make it decay more rapidly . This relationship is evident in Figure 2.1.2.
; ; ; !1 For applications to numerical methods for partial di erential equations, there are two properties of the Fourier transform that are most important. 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 86 One is equation (2.1.8): the Fourier transform converts convolution into multiplication. The second can be derived by integration by parts: ;i xux(x)dx = ; Z 1 (;i )e;i xu(x)dx = i u( ) (2:1:16) c ux( ) = ^ ;1 e ;1 Z1 assuming u(x) is smooth and decays at 1. That is, the Fourier transform converts di erentiation into multiplication by i . This result is rigorously valid for any absolutely continuous function u 2 L2 whose derivative belongs to L2 . Note that di erentiation makes a function less smooth, so the fact that it makes the Fourier transform decay less rapidly ts the pattern mentioned above for convolution. ; 2 0 2 Figure 2.1.3.
EXAMPLE 2.1.2. The function 8 > < u(x) = > : 1 4 1 ; 4 for 2 x < 0, for 0 < x 2, 0 otherwise,
; (2:1:17) illustrated in Figure 2.1.3, has Fourier transform
1 u( ) = 4 ^ Z0 1 = 4i (1 e2i
; ; ;2 e;i x dx ; 1 4 Z2
0 e;i x dx
4i
2 ;i )2 = i sin ;e ;2i +1) = 1 (ei ;e (2:1:18) which is i times the Fourier transform (2.1.14) of the triangular hat function (2.1.12). In keeping with (2.1.16), (2.1.17) is the derivative of (2.1.12). The following theorem collects (2.1.16) together with a number of additional properties of the Fourier transform: 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 87 PROPERTIES OF THE FOURIER TRANSFORM Theorem 2.2. Let u v 2 L2 have Fourier transforms u = F u, v = F v. ^ ^
Then: (a) Linearity. Ffu + vg( ) = u( )+^( ) Ffcug( ) = cu( ): ^ v ^ (b) Translation. If x0 2 R then Ffu(x + x0)g( ) = ei x0 u( ): ^ (c) Modulation. If 0 2 R then Ffei 0xu(x)g( ) = u( ; 0 ): ^ (d) Dilation. If c 2 R with c 6= 0 then Ffu(cx)g( ) = u( =c)=jcj: ^ ^ (e) Conjugation. Ff u g( ) = u(; ): (f) Di erentiation. If ux 2 L2 then Ffuxg( ) = i u( ): ^ (g) Inversion. F ;1fug( ) = 21 u(; ): ^ Proof. See Exercise 2.1.2. In particular, taking c = ;1 in part (d) above gives Ffu(;x)g = u(; ). ^ Combining this result with part (e) leads to the following elementary but useful results. De nitions: u(x) is even, odd, real, or imaginary if u(x) = u(;x), u(x) = ;u(;x), u(x) = u(x), or u(x) = ;u(x), respectively u(x) is hermitian or skewhermitian if u(x) = u(;x) or u(x) = ;u(;x), respectively.
SYMMETRIES OF THE FOURIER TRANSFORM Theorem 2.3. Let u 2 L2 have Fourier transform u = F u. Then ^ (a) u(x) is even (odd) () u( ) is even (odd) ^ (b) u(x) is real (imaginary) () u( ) is hermitian (skewhermitian) ^ and therefore (c) u(x) is real and even () u( ) is real and even ^ (d) u(x) is real and odd () u( ) is imaginary and odd ^ (e) u(x) is imaginary and even () u( ) is imaginary and even ^ (f) u(x) is imaginary and odd () u( ) is real and odd. ^ Proof. See Exercise 2.1.3. 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 88 In the discussion above we have twice observed the following relationships between the smoothness of a function and the decay of its Fourier transform: u( ) ^ u(x) smooth $ decays rapidly as j j ! 1 decays rapidly as jxj ! 1 $ smooth (Of course since the Fourier transform is essentially the same as the inverse Fourier transform, by Theorem 2.2g, the two rows of this summary are equivalent.) The intuitive explanation is that if a function is smooth, then it can be accurately represented as a superposition of slowlyvarying waves, so one does not need much energy in the high wave number components. Conversely, a nonsmooth function requires a considerable amplitude of high wave number components to be represented accurately. These relationships are the bedrock of analog and digital signal processing, where all kinds of smoothing operations are e ected by multiplying the Fourier transform by a \windowing function" that decays suitably rapidly. The following theorem makes these connections between u and u precise. ^ This theorem may seem forbidding at rst, but it is worth studying carefully. Each of the four parts of the theorem makes a stronger smoothness assumption on u than the last, and reaches a correspondingly stronger conclusion about the rate of decay of u( ) as j j ! 1. Parts (c) and (d) are known as the ^ PaleyWiener theorems. First, a standard de nition. A function u de ned on R is said to have bounded variation if there isPamconstant M such that for any nite m and any points x0 < x1 < < xm, j=1 ju(xj ) ; u(xj;1 )j M . 2.1. THE FOURIER TRANSFORM TREFETHEN 1994 89 Theorem 2.4. Let u be a function in L2 .
u( ) = O(j j;p;1) ^ SMOOTHNESS OF u AND DECAY OF u ^ (a) If u has p ; 1 continuous derivatives in L2 for some p 0, and a pth derivative in L2 that has bounded variation, then as j j ! 1: (2:1:19) (2:1:20) (b) If u has in nitely many continuous derivatives in L2 , then u( ) = O(j j;M ) ^ as j j ! 1 for all M and conversely. (c) If u can be extended to an analytic function of z = x + iy in the complex strip j Im zj < a for some a > 0, with ku(x + iy)k const uniformly for each constant ;a < y < a, then eaj j u( ) 2 L2 ^ (2:1:21) and conversely. (d) If u can be extended to an entire function* of z = x + iy with ju(z)j = O(eajzj) as jzj ! 1 (z 2 C ) for some a > 0, then u has compact support ^ contained in ;a a], i.e. u( ) = 0 ^
and conversely. for all j j > a (2:1:22) Proof. See, for example, xVI.7 of Y. Katznelson, An Introduction to Harmonic Analysis, Dover, 1976. Also see Rudin (p. 407), Paley & Wiener, Reed & Simon v. 2, Hormander v. 1 (p. 181), entire functions books::: ] A function of the kind described in (d) is said to be bandlimited, since only a nite band of wave numbers are represented in it. Since the Fourier transform and its inverse are essentially the same, by Theorem 2.2g, Theorem 2.4 also applies if the roles of u(x) and u( ) are inter^ changed.
EXAMPLE 2.1.1, CONTINUED. The square wave u of Example 2.1.1 (Figure 2.1.2) satis es condition (a) of Theorem 2.4 with p = 0, so its Fourier transform should satisfy
*An entire function is a function that is analytic throughout the complex plane C . 2.1. THE FOURIER TRANSFORM
j j j j TREFETHEN 1994 90 u( ) = O( ;1 ), as is veri ed by (2.1.10). On the other hand, suppose we interchange ^ the roles of u and u and apply the theorem again. The function u( ) = sin = is entire, ^ and since sin( ) = (ei e;i )=2i, it satis es u( ) = O(ej j ) as (with now taking complex values). By part (d) of Theorem 2.4, it follows that u(x) must have compact
; j j!1 support contained in 1 1], as indeed it does. Repeating the example for u u, condition (a) now applies with p = 1, and the Fourier transform (2.1.14) is indeed of magnitude O( ;2 ), as required. Interchanging u and u, ^ we note that sin2 = 2 is an entire function of magnitude O(e2j j ) as , and u u has support contained in 2 2].
; j j j j!1 ; . 2.1.1. Show that the two integrals in the de nition (2.1.4) of u v are equivalent. . 2.1.2. Derive conditions (a){(g) of Theorem 2.2. (Do not worry about justifying the usual EXERCISES operations on integrals.) . 2.1.3. Prove Theorem 2.3. . 2.1.4. (a) Which functions u L2 L1 satisfy u u = 0? (b) How about u u = u? . 2.1.5. Integration. (a) What does partR(f) of Theorem 2.2 suggest should be the Fourier transform of the x function U (x) = ;1 u(s)ds? R1 (b) Obviously U (x) cannot belong to L2 unless ;1 u(x)dx = 0, so by Theorem 2.1, this is R1 ^ a necessary condition for U to be in L2 also. Explain how the condition ;1 u(x)dx = 0 ^ relates to your formula of (a) for U in terms of u. ^ . 2.1.6. (a) Calculate the Fourier transform of u(x) = (1 + x2 );1 . (Hint: use a complex contour integral if you know how. Otherwise, look the result up in a table of integrals.) (b) How does this example t into the framework of Theorem 2.4? Which parts of the theorem apply to u? (c) If the roles of u and u are interchanged, how does the example now t Theorem 2.4? ^ Which parts of the theorem apply to u? ^ . 2.1.7. The autocorrelation function of a function u L2 L1 may be de ned by
2 \ 1 Z 1 u(x)u(x + c)dx: (c) = u 2
k k 2 \ ;1 Find an expression for (c) as an inverse Fourier transform of a product of Fourier transforms involving u. This expression is the basis of some algorithms for computing (c). 2.1. THE FOURIER TRANSFORM Fourier transform of the following function: TREFETHEN 1994 91 . 2.1.8. Without evaluating any integrals, use Theorem 2.2 and (2.1.10) to determine the 4 2 0 2 4 . 2.1.9. The uncertainty principle.2 Show by using Theorem 2.4 that if u(x) and u( ) both ^ have compact support, with u L , then u(x) 0.
; ; 2 2.2. THE SEMIDISCRETE FOURIER TRANSFORM TREFETHEN 1994 92 2.2. The semidiscrete Fourier transform
The semidiscrete Fourier transform is the reverse of the more familiar Fourier series: instead of a bounded, continuous spatial domain and an unbounded, discrete transform domain, it involves the opposite. This is just what is needed for the analysis of numerical methods for partial di erential equations, where we are perpetually concerned with functions de ned on discrete grids. For many analytical purposes it is simplest to think of these grids as in nite in extent. Let h > 0 be a real number, the space step, and let ::: x;1 x0 x1 ::: be de ned by xj = jh. Thus xj = hZ, where Z is the set of integers. We are concerned now with spatial grid functions v = vj , which may or may not be approximations to a continuous function u, vj u(xj ): As in the last chapter, it will be convenient to write v(xj ) sometimes for vj .
f g f g o * o * o * o * o * h o * o * o * x ; =h 0 =h Figure 2.2.1. Space and wave number domains for the semidiscrete Fourier
transform. L by a lowercase script letter `. (Both symbols honor Henri Lebesgue, the mathematician
1 X
j =;1
1=2 For functions de ned on discrete domains it is standard to replace the uppercase letter who laid the foundations of modern functional analysis early in the twentieth century.) The `2 norm of a grid function v is the nonnegative or in nite real number h
k k v = h j vj 2
j : (2:2:1) Notice the h in front of the summation. One can think of (2.2.1) as a discrete approximation to the integral (2.1.1) by the trapezoid rule or the rectangle rule for quadrature. (On an unbounded domain these two are equivalent.) The symbol `2 (\little Ltwo sub h") denotes h the set of grid functions of nite norm, `2 = v : v < h
f k k 1g and similarly with `1 and `1 . In contrast to the situation with L1 , L2, and L1 , these h h spaces are nested: `1 `2 `1 : (2:2:2) h h h (See Exercise 2.1.1.) 2.2. THE SEMIDISCRETE FOURIER TRANSFORM TREFETHEN 1994 93 The convolution v w of two functions v w is the function v w de ned by (v w)m = h
1 X j =;1 vm;j wj = h 1 X
j =;1 vj wm;j (2:2:3) provided that these sums exist. This formula is a trapezoid or rectangle rule approximation to (2.1.4). For any v `2 , the semidiscrete Fourier transform of v is the function v( ) de ned ^ h by
2 v ( ) = ( h v)( ) = h ^
F j =;1 1 X ;i xj e vj 2 ; =h =h]
2 a discrete approximation to (2.1.5). A priori, this sum de nes a function v( ) for all R . ^ However, notice that for any integer m, the exponential e2 imxj =h = e2 imj is exactly 1 at all of the grid points xj . More generally, any wave number is indistinguishable on the grid from all other wave numbers +2 m=h, where m is any integera phenomenon called aliasing. This means that the function v( ) is 2 =hperiodic on R . To make sense of the ^ idea of analyzing v into component oscillations, we shall normally restrict attention to one period of v by looking only at wave numbers in the range =h =h], and it is in this sense ^ that the Fourier transform of a grid function is de ned on a bounded domain. But the reader should bear in mind that the restriction of to any particular interval is a matter of labeling, not mathematics in principle e0 and e100 ij are equally valid representations of the grid function vj 1. Thus for discretized functions v, the transform v( ) inhabits a bounded domain. On ^ the other hand the domain is still continuous. This re ects the fact that arbitrarily ne gradations of wave number are distinguishable on an unbounded grid. Since x and belong to di erent sets, it is necessary to de ne an additional vector space for functions v. The L2 norm of a function v is the number ^ ^ h
; k k v = ^ hZ =h ; =h j v( ) 2 d ^
j i1=2 : (2:2:4) One can think of this as an approximation to (2.1.1) in which the wave number components with > =h have been replaced by zero. The symbol L2 denotes the set of (Lebesgueh measurable) functions on =h =h] of nite norm,
j j ; L2 = v : v < ^ ^ h
f k k 1g : (2:2:5) Now we can state a theorem analogous to Theorem 2.1: 2.2. THE SEMIDISCRETE FOURIER TRANSFORM TREFETHEN 1994 94 THE SEMIDISCRETE FOURIER TRANSFORM Theorem 2.5. If v `2 , then the semidiscrete Fourier transform h
2 v( ) = ( h v)( ) = h ^
F j =;1 1 X ;i xj e vj 2 ; =h =h] (2:2:6) belongs to L2 , and v can be recovered from v by the inverse semidiscrete Fourier ^ h transform i xj ^ j Z: 2 ; =h e v( )d The `2 norm of v and the L2 norm of v are related by Parseval's equality, ^ h h
2 k k ;^ vj = (Fh 1 v)(x) = 1 Z =h (2:2:7) (2:2:8) (2:2:9) v = 2 v: ^
k k 2 p If u `2 and v `1 (or vice versa), then v w `2 , and vd satis es w h h h
2 2 vd ( ) = v( )w ( ): w ^ ^ As in the continuous case, the following properties of the semidiscrete Fourier transform will be useful. In (c), and throughout this book wherever convenient, we take advantage of the fact that v( ) can be treated as a periodic function de ned for all R . ^
2 2.2. THE SEMIDISCRETE FOURIER TRANSFORM TREFETHEN 1994 95 PROPERTIES OF THE SEMIDISCRETE FOURIER TRANSFORM Theorem 2.6. Let v w `2 have Fourier transforms v w. Then: ^ ^ h
2 (a) Linearity. h v + w ( ) = v( )+ w ( ) ^ ^
F f g F f h cv ( ) = cv( ): ^
g (b) Translation. If j0 (c) Modulation. If (d) Dilation. If m
2 2 Z then h vj+j0 ( ) = ei xj0 v( ): ^
F f g 02R ^ then h ei 0 xj vj ( ) = v(
F f g 6 F f g ; 0 ):
j j Z with
g m = 0 then h vmj ( ) = v( =m)= m : ^
; (e) Conjugation. h v ( ) = v( ): ^
F f The symmetry properties of the Fourier transform summarized in Theorem 2.3 apply to the semidiscrete Fourier transform too we shall not repeat the list here. We come now to a fundamental result that describes the relationship of the Fourier transform of a continuous function u to that of a discretization v of uor if x and are interchanged, the relationship of Fourier series to Fourier transforms. Recall that because of the phenomenon of aliasing, all wave numbers +2 j=h, j Z, are indistinguishable on the grid hZ. Suppose that u L2 is a su ciently smooth function de ned on R , and let v `2 be the discretization obtained by sampling u(x) at the points xj . The aliasing h principle implies that v( ) must consist of the sum of all of the values u( +2 j=h). This ^ ^ result is known as the Poisson summation formula or the aliasing formula:
2 2 2 ALIASING FORMULA Theorem 2.7. Let u L2 be su ciently smooth ?], with Fourier transform u, and let ^
v `2 be the restriction of u to the grid hZ. Then h
2 2 v( ) = ^ 1 X
j =;1 u( +2 j=h) ^ 2 ; =h =h]: (2:2:10) Proof. Not yet written. See P. Henrici, Applied and Computational Complex Analysis, v. 3, Wiley, 1986. In applications, we are very often concerned with functions v obtained by discretization, and it will be useful to know how much the Fourier transform is a ected in the process. Theorems 2.4 and 2.7 combine to give the following answers to this question: 2.2. THE SEMIDISCRETE FOURIER TRANSFORM TREFETHEN 1994 96 EFFECT OF DISCRETIZATION ON THE FOURIER TRANSFORM Theorem 2.8. Let v be the restriction to the grid hZ of a function u L2. The following
2 2 ; ; ; estimates hold uniformly for all =h =h], or a forteriori, for in any xed interval A A]. (a) If u has p 1 continuous derivatives in L2 for some p 1 ?], and a pth derivative in L2 that has bounded variation, then v( ) u( ) = O(hp+1 ) as h 0: ^ ^ (2:2:11)
j ; j ! (b) If u has in nitely many continuous derivatives in L2 , then ^ v( ) u( ) = O(hM ) as h 0 for all M: ^
j ; j ! j k k ; (2:2:12)
j (c) If u can be extended to an analytic function of z = x+iy in the complex strip Im z < a for some a > 0, with u( + iy) const uniformly for each a < y < a, then for any > 0, v( ) u( ) = O(e; (a; )=h ) as h 0: ^ ^ (2:2:13)
j ; j ! (d) If u can be extended to an entire function of z = x + iy with u(z ) = O(eajzj) as z (z C ) for some a > 0, then v( ) = u( ) provided h < =a: ^ ^ (2:2:14)
j j!1 2 In part (c), u( +iy) denotes a function of x, namely u(x+iy) with x interpreted as a variable and y as a xed parameter. Proof. In each part of the theorem, u(x) is smooth enough for Theorem 2.7 to apply, which gives the identity v ( ) u( ) = ^ ^
; 2 ; 1 X
j =1 u( +2 j=h)+ u( ^ ^ ; 2 j=h): (2:2:15) Note that since =h =h], the arguments of u in this series have magnitudes at least ^ =h, 2 =h, 3 =h ::: . ^ For part (a), Theorem 2.4(a) asserts that u( ) C1 ;p;1 for some constant C1 . With (2.2.15) this implies
j j j j j v ( ) u( ) ^ ^
; j C1 1 X
j =1 1 ;p;1 = C hp+1 X j ;p;1 : (j =h) 2
j =1 Since p 1 this sum converges to a constant, which implies (2.2.11) as required. Part (b) follows from part (a). For part (c), ::: ?] For part (d), note that if h < =a, then =h > a. Thus (2.2.15) reduces to 0 for all =h =h], as claimed.
2 ; Note that part (d) of Theorem 2.8 asserts that on a grid of size h, the semidiscrete Fourier transform is exact for bandlimited functions containing energy only at wave numbers 2.2. THE SEMIDISCRETE FOURIER TRANSFORM
j j TREFETHEN 1994 97 smaller than =hthe Nyquist wave number, corresponding to two grid points per wavelength. This twopointsperwavelength restriction is famous among engineers, and has practical consequences in everything from ghter planes to compact disc players. When we come to discretize solutions of partial di erential equations, two points per wavelength will be the coarsest resolution we can hope for under normal circumstances. (a) Prove (2.2.2): `1 `2 `1. h h h (b) Give examples to show that these inclusions are proper: `2 `1 and `1 `2 . h h h h (c) Give examples to show that neither inclusion in (a) holds for functions on continuous domains: L1 L2 and L2 L1. . 2.2.2. Let : `2 `2 be the discrete di erentiation and smoothing operators de ned by h h
6 6 6 6 ! . 2.2.1. EXERCISES ( v)j = 21 (vj+1 vj;1 ) h
; 1 ( v)j = 2 (vj;1 + vj+1 ): (2:2:16)
2 (a) Show that and are equivalent to convolutions with appropriate sequences d m `2 . h (Be careful with factors of h.) ^ ^ (b) Compute the Fourier transforms d and m. How does d compare to the transform of ^ the exact di erentiation operator for functions de ned on R (Theorem 2.2f)? Illustrate ^ this comparison with a sketch of d( ) against . ^ , m , and m , and verify Parseval's equality. (c) Compute d , d ^ (d) Compute the Fourier transforms of the convolution sequences corresponding to the iterated operators p and p (p 2). Discuss how these results relate to the rule of thumb discussed in the last section: the smoother the function, the more rapidly its Fourier transform decays as . What imperfection in does this analysis bring to light? . 2.2.3. Continuation 2of;Exercise 2.1.6. Let v be the discretization on the grid hZ of the function u(x) = (1+ x ) 1 . (a) Determine v( ). (Hint: calculating it from the de nition (2.2.6) is very di cult.) ^ (b) How fast does v( ) approach u( ) as h 0? Give a precise answer based on (a), then ^ ^ compare your answer with the prediction of Theorem 2.8. (c) What would the answer to (b) have been if the roles of u and u had been interchanged ^ that is, if v had been the discretization not of u(x) but of its Fourier transform? . 2.2.4. Integration by the trapezoid rule. A function u L2 L1 can be integrated approximately by the trapezoid rule:
k k k k k k k k j j!1 ! 2 \ I=
j j ! 1 Z1 ;1 u(x) dx Ih = h 1 X
j =;1 u(xj ): (2:2:17) This is an in nite sum, but in practice one might delete the tails if u decays su ciently rapidly as x . (This idea leads to excellent quadrature algorithms even for nite intervals, which are rst transformed to the real axis by a change of variables for a survey 2.2. THE SEMIDISCRETE FOURIER TRANSFORM TREFETHEN 1994 98 see M. Mori, \Quadrature formulas obtained by variable transformation and the DErule," J. Comp. Appl. Math. 12 & 13 (1985), 119{130.) As h 0, how good an approximation is Ih to the exact integral I ? Of course the answer will depend on the smoothness of u(x). (a) State how Ih is related to the semidiscrete Fourier transform. (b) Give a bound for Ih I based on the theorems of this section. (c) In particular, what can you say about Ih I for the function u(x) = e;x2 ? (d) Show that your bound can be improved in a certain sense by a factor of 2. . . 2.2.5. Draw a plot of sin n as a function of n, where n ranges over the integers (not the real numbers) from 1 to 5000. (That is, your plot should contain 5000 dots in Matlab this can be done in one line.) Explain why the plot looks the way it does to the human eye, and what this has to do with aliasing. Make your explanation precise and quantitative. (See G. Strang, Calculus, WellesleyCambridge Press, 1991.)
! j ; j j ; j 2.3. INTERPOLATION AND SINC FUNCTIONS TREFETHEN 1994 99 2.3. Interpolation and sinc functions
This section is not written yet, but here's a hint as to what will be in it.] If j is the Kronecker delta function
j = 0 if j 6= 0, 1 if j = 0, (2:3:1) then (2.2.6) gives the semidiscrete Fourier transform ^j ( ) = h
j=
6 (for all ): If we now apply the inverse transform formula (2.2.7), we nd after a little algebra sin( xj =h) xj =h (2:3:2)
6 at least for j = 0. Since xj =h is a nonzero integer for each j = 0, the sines are zero and this formula matches (2.3.1). Suppose, however, that we evaluate (2.3.2) not just for x = xj but for all values x R . Then we've got a sinc function again, one that can be called a grid sinc function:
2 x=h Sh (x) = sin(x=h ) : (2:3:3) The plot of Sh (x) is the same as the upperright plot of Figure 2.1.2, except scaled so that the zeros are on the grid (i.e. at integer multiples of h). Obviously Sh (x) is a continuous interpolant to the discrete delta function j . Which one? It is the unique bandlimited c interpolant, bandlimited in the sense that its Fourier transform Sh( ) is zero for = =h =h]. (Proof: by construction it's bandlimited in that way, and uniqueness can be proved via an argument by contradiction, making use of Parseval's equality (2.2.8).) More generally, suppose we have an arbitrary grid function vj (well, not quite arbitrary we'll need certain integrability assumptions, but let's forget that for now). Then the bandlimited interpolant to vj is the unique function v(x) de ned for x R with v(xj ) = vj and v ( ) = 0 for = =h =h]. It can be derived in two equivalent ways: ^ Method 1: Fourier transform. Given vj , compute the semidiscrete Fourier transform v( ). Then invert that transform, and evaluate the resulting formula for all x rather than ^ just on the grid. Method 2: linear combination of sinc functions. Write
2 ; 2 2 ; vj = 1 X
m=;1 vm m;j 2.3. INTERPOLATION AND SINC FUNCTIONS and then set TREFETHEN 1994 100 v(x) = 1 X
m=;1 vm Sh (x xm ):
; The equivalence of Methods 1 and 2 is trivial it follows from the linearity and translationinvariance of all the processes in question. The consideration of bandlimited interpolation is a good way to get insight into the Aliasing Formula presented as Theorem 2.7. (In fact, maybe that should go in this section.) The following schema summarizes everything. Study it! u(x)
# F.T.
; ! DISCRETIZE vj
# F.T.
; ! BANDLIMITED INTERPOLATION F.T. v (x )
; ! u( ) ^ ALIASING FORMULA v( ) ^ ZERO HIGH WAVE NOS. v( ) ^
l # The Gibbs phenomenon is a particular phenomenon of bandlimited interpolation that has received much attention. After an initial discovery by Wilbraham in 1848, it was made famous by Michelson in 1898 in a letter to Nature, and then by an analysis by Gibbs in Nature the next year. Gibbs showed that if the step function u(x) = +1 x < 0, 1 x>0
; The ringing is scaleinvariant it does not go away as h 0. In the nal text I will illustrate the Gibbs phenomenon and include a quick derivation of (2.3.4).
! is sampled on a grid and then interpolated in the bandlimited manner, then the resulting function v(x) exhibits a ringing e ect: it overshoots the limits 1 by about 9%, achieving a maximum amplitude Z 1 sin( y) (2:3:4) y dy 1:089490:
;1 2.4. THE DISCRETE FOURIER TRANSFORM TREFETHEN 1994 101 Note: although the results of the last two sections will be used throughout the remainder of the book, the material of the present section will not be needed until Chapters 8 and 9. For the discrete Fourier transform, both x and inhabit discrete, bounded domains or if your prefer, they are periodic functions de ned on discrete, unbounded domains. Thus there is a pleasing symmetry here, as with the Fourier transform, that was missing in the semidiscrete case.
N ;N = ; 2 2
o * 2.4. The discrete Fourier transform x; N = 2
o * ; x0 = 0
o * o * o * o * h o * o * xN = 2
o * o * x 0
o * o * =0
o * o * 1 o * o * o * N 2 =N 2
o * Figure 2.4.1. Space and wave number domains for the discrete Fourier transform. For the fundamental spatial domain we shall take Let N be a positive even integer, set h= 2 (N even)
; ), as illustrated in Figure 2.4.1. (2:4:1)
; N= : (2:4:2) 2 h Let `2 denote the set of functions on xj that are N periodic with respect to j , i.e, N 2 periodic with respect to x, with the norm h N=2;1 2 i1=2 X v = h vj : (2:4:3)
f g k k and de ne xj = jh for any j . The grid points in the fundamental domain are x;N=2 = ::: x0 = 0 ::: xN=2;1 = h: An invaluable identity to keep in mind is this:
; N j =;N=2 j j (Since the sum is nite, the norm is nite, so every function of the required type is guaranteed to belong to `2 and to `1 and `1 .) The discrete Fourier transform (DFT) of a N N N function v `2 is de ned by N
2 v( ) = ( N v)( ) = h ^
F N=2;1 X j =;N=2 e;i xj vj 2 Z: 2.4. THE DISCRETE FOURIER TRANSFORM TREFETHEN 1994 102 Since the spatial domain is periodic, the set of wave numbers is discrete, and in fact ranges precisely over the set of integers Z. Thus it is natural to use as a subscript, v = ( N v) = h ^
F N=2;1 X j =;N=2 e;i jh vj 2 Z and since h = 2 =N , v is N periodic as a function of . We shall take N=2 N=2] as the ^ fundamental domain of wave numbers, and let L2 denote the set of N periodic functions N on the grid Z, with norm
; k k v = ^ h N=2;1 X This is nonstandard notation, for an upper case L is normally reserved for a family of functions de ned on a continuum. We use it here to highlight the relationship of the discrete Fourier transform with the semidiscrete Fourier transform. The convolution of two functions in `2 is de ned by N (v w)m = h
N=2;1 X j =;N=2 =;N=2 j v ^ j 2 i1=2 : (2:4:4) vm;j wj = h N=2;1 X j =;N=2 vj wm;j : (2:4:5) Again, since the sum is nite, there is no question of convergence. Here is a summary of the discrete Fourier transform: 2.4. THE DISCRETE FOURIER TRANSFORM TREFETHEN 1994 103 THE DISCRETE FOURIER TRANSFORM Theorem 2.9. If v `2 , then the discrete Fourier transform N
2 v = ( N v) = h ^
F N=2;1 X j =;N=2 e;i jh vj N
; N 2 2 ; 1 (2:4:6) belongs to L2 , and v can be recovered from v by the inverse discrete Fourier trans^ N form The `2 norm of v and the L2 norm of v are related by Parseval's equality, ^ N N
k k N=2;1 ;1 v) = 1 X ei jh v : ^ vj = (FN ^ j 2 =;N=2 (2:4:7) v = 2 v: ^
k k p (2:4:8) (2:4:9) If v w `2 , then vd satis es w N
2 (vd ) = v w : w ^ ^ As with the other Fourier transforms we have considered, the following properties of the discrete Fourier transform will be useful. Once again we take advantage of the fact that v( ) can be treated as a periodic function de ned for all Z. ^
2 PROPERTIES OF THE DISCRETE FOURIER TRANSFORM Theorem 2.10. Let v w `2 have discrete Fourier transforms v w. Then: ^ ^ N
2 (a) Linearity. N v + w ( ) = v( )+ w ( ) ^ ^
F f g F ^ N fcv g( ) = cv( ): ): (b) Translation. If j0 (c) Modulation. If
F 2 Z then N vj+j0 ( ) = ei xj0 v( ): ^
F f g 02Z
f g ^ then N ei 0 xj vj ( ) = v(
F f g ; ; 0 (e) Conjugation. N v ( ) = v( ): ^
; (g) Inversion. N 1 v ( ) =
F f g 1 2
h v( ): ^
; An enormously important fact about discrete Fourier transforms is that they can be computed rapidly by the recursive algorithm known as the fast Fourier transform (FFT).* A direct implementation of (2.4.6) or (2.4.7) requires (N 2 ) arithmetic operations, but the
*The fast Fourier transform was discovered by Gauss in 1805 at the age of 28, but although he wrote a paper on the subject, he did not publish it, and the idea was more or less lost until its celebrated 2.4. THE DISCRETE FOURIER TRANSFORM TREFETHEN 1994 104 FFT is based upon a recursion that reduces this gure to (N log N ). We shall not describe the details of the FFT here, but refer the reader to various books in numerical analysis, signal processing, or other elds. However, to illustrate how simple an implementation of this idea may be, Figure 2.4.2 reproduces the original Fortran program that appeared in a 1965 paper ; by Cooley, Lewis, and Welch. Assuming that N is a power of 2, it computes 2 N 1 , in our notation: the vector A(1 : N ) represents v0 ::: vN ;1 on input and 2 v0 ::: 2 vN ;1 on ^ ^ output.
y F subroutine fft(a,m) complex a(1),u,w,t n = 2**m nv2 = n/2 nm1 = n1 j=1 do 7 i = 1,nm1 if (i.ge.j) goto 5 t = a(j) a(j) = a(i) a(i) = t 5 k = nv2 6 if (k.ge.j) goto 7 j = jk k = k/2 goto 6 7 j = j+k do 20 l = 1,m le = 2**l le1 = le/2 u = 1. ang = 3.14159265358979/le1 w = cmplx(cos(ang),sin(ang)) do 20 j = 1,le1 do 10 i = j,n,le ip = i+le1 t = a(ip)*u a(ip) = a(i)t 10 a(i) = a(i)+t 20 u = u*w return end (1965). Figure 2.4.2. Complex inverse FFT program of Cooley, Lewis, and Welch As mentioned above, this program computes the inverse Fourier transform according to our de nitions, times 2 . The same program can be used for the forward transform by making use of the following identity: rediscovery by Cooley and Tukey in 1965. (See M. T. Heideman, et al., \Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine, October 1984.) Since then, fast Fourier transforms have changed prevailing computational practices in many areas. y Before publication, permission to print this program will be secured. 2.4. THE DISCRETE FOURIER TRANSFORM
; v = N v ( ) = 2 h N 1 v ( ): ^
F f g ; F f g TREFETHEN 1994 105 (2:4:10) These equalities follow from parts (e) and (g) of Theorem 2.10, respectively. 2.5. VECTORS AND MULTIPLE SPACE DIMENSIONS TREFETHEN 1994 106 2.5. Vectors and multiple space dimensions
Fourier analysis generalizes with surprising ease to situations where the independent variable x and/or the dependent variable u are vectors. We shall only sketch the essentials, which are based on the following two ideas: If x is a dvector, then the dual variable is a dvector too, and the Fourier integral is a multiple integral involving the inner product x If u is an N vector, then its Fourier transform u is an N vector too, and is de ned ^ componentwise. As these statements suggest, our notation will be as follows: d = number of space dimensions: x = (x1 ::: N = number of dependent variables: u = (u1 :::
Both and u become vectors of the same dimensions, ^ = ( 1 ::: d )T and x becomes the dot product transform and its inverse read
F xd )T uN )T : u = (^1 ::: uN )T ^ u ^ x = 1 x1 + + d xd . The formulas for the Fourier
(2:5:1) Z u( ) = ( u)( ) = e;i x u(x) dx ^ Z1 Z1 = e;i x u(x) dx1 dxd
;1 ;1 for 2 R d, and Z u(x) = ( ;1 u)(x) = (2 );d ei x u( )d ^ ^ Z1 Z1 ;d = (2 ) ei xu( ) d 1 d d ^
F (2:5:2) ;1 ;1 for x 2 R d. In other words, u and u are related componentwise: ^ u( ) = (^(1) ( ) ::: u(N ) ( ))T : ^ u ^
If the vector L2 norm is de ned by
k (2:5:3) u =
k 2 Z k u(x) dx =
k 2 Z1 Z1
;1 ;1
k u(x) 2 dx1 dxd
k (2:5:4) 2.5. VECTORS AND MULTIPLE SPACE DIMENSIONS where the symbol
k k TREFETHEN 1994 107 Parseval's equality for vector Fourier transforms takes the form
k in the integrand denotes the 2norm on vectors of length N , then u = (2 )d=2 u : ^
k k k (2:5:5) The set of vector functions with bounded vector 2norms can be written simply as (L2 )N . Before speaking of convolutions, we have to go a step further and allow u(x) and u( ) ^ to be M N matrices rather than just N vectors. The de nitions above extend to this further case unchanged, if the symbol in the integrand of (2.5.4) now represents the 2norm (largest singular value) of a matrix. If u(x) is an M P matrix and v(x) is a P N matrix, then the convolution u v is de ned by
k k (u v)(x) = = Z Z u(x y)v(y) dy
; ; u(y)v(x y) dy Z1 Z1 = u(x y)v(y) dy1 dyd ;1 ;1
; (2:5:6) (2:5:7) Since matrices do not commute in general, it is no longer possible to exchange u and v as in (2.1.4). This generalization of Fourier transforms and convolutions to matrix functions is far from idle, for we shall need it for the Fourier analysis of multistep nite di erence approximations such as the leap frog formula. Similar generalizations of our scalar results hold for semidiscrete and discrete Fourier transforms. and it satis es d u v( ) = u( )^( ): ^ v . 2.5.1. What is the Fourier transform of the vector function
T u(x) = sin x sin2x x 2x EXERCISES de ned for x R ? . 2.5.2. What is the Fourier transform of the scalar function
2 1 1 2 u(x) = e; 2 (x2 +x2 ) de ned for x = (x1 x2 )T 2 R2? ...
View
Full Document
 Spring '10
 RunyiYu
 Fourier Series

Click to edit the document details