13 Chapter Fourier Analysis
In addition to their inestimable importance in mathematics and its applications, Fourier series also serve as the entry point into the wonderful world of Fourier analysis and its wide-ranging extensions and generalizations. An entire industry is devoted to further developing the theory and enlarging the scope of applications of Fourierinspired methods. New directions in Fourier analysis continue to be discovered and exploited in a broad range of physical, mathematical, engineering, chemical, biological, financial, and other systems. In this chapter, we will concentrate on four of the most important variants: discrete Fourier sums leading to the Fast Fourier Transform (FFT); the modern theory of wavelets; the Fourier transform; and, finally, its cousin, the Laplace transform. In addition, more general types of eigenfunction expansions associated with partial differential equations in higher dimensions will appear in the following chapters. Modern digital media, such as CD's, DVD's and MP3's, are based on discrete data, not continuous functions. One typically samples an analog signal at equally spaced time intervals, and then works exclusively with the resulting discrete (digital) data. The associated discrete Fourier representation re-expresses the data in terms of sampled complex exponentials; it can, in fact, be handled by finite-dimensional vector space methods, and so, technically, belongs back in the linear algebra portion of this text. However, the insight gained from the classical continuous Fourier theory proves to be essential in understanding and analyzing its discrete digital counterpart. An important application of discrete Fourier sums is in signal and image processing. Basic data compression and noise removal algorithms are applied to the sample's discrete Fourier coefficients, acting on the observation that noise tends to accumulate in the high frequency Fourier modes, while most important features are concentrated at low frequencies. The first Section 13.1 develops the basic Fourier theory in this discrete setting, culminating in the Fast Fourier Transform (FFT), which produces an efficient numerical algorithm for passing between a signal and its discrete Fourier coefficients. One of the inherent limitations of classical Fourier methods, both continuous and discrete, is that they are not well adapted to localized data. (In physics, this lack of localization is the basis of the Heisenberg Uncertainty Principle.) As a result, Fourier-based signal processing algorithms tend to be inaccurate and/or inefficient when confronting highly localized signals or images. In the second section, we introduce the modern theory of wavelets, which is a recent extension of Fourier analysis that more naturally incorporates multiple scales and localization. Wavelets are playing an increasingly dominant role in many modern applications; for instance, the new JPEG digital image compression format is based on wavelets, as are the computerized FBI fingerprint data used in law enforcement 2/25/07 691
c 2006 Peter J. Olver
in the United States. Spectral analysis of non-periodic functions defined on the entire real line requires replacing the Fourier series by a limiting Fourier integral. The resulting Fourier transform plays an essential role in functional analysis, in the analysis of ordinary and partial differential equations, and in quantum mechanics, data analysis, signal processing, and many other applied fields. In Section 13.3, we introduce the most important applied features of the Fourier transform. The closely related Laplace transform is a basic tool in engineering applications. To mathematicians, the Fourier transform is the more fundamental of the two, while the Laplace transform is viewed as a certain real specialization. Both transforms change differentiation into multiplication, thereby converting linear differential equations into algebraic equations. The Fourier transform is primarily used for solving boundary value problems on the real line, while initial value problems, particularly those involving discontinuous forcing terms, are effectively handled by the Laplace transform.
13.1. Discrete Fourier Analysis and the Fast Fourier Transform.
In modern digital media -- audio, still images or video -- continuous signals are sampled at discrete time intervals before being processed. Fourier analysis decomposes the sampled signal into its fundamental periodic constituents -- sines and cosines, or, more conveniently, complex exponentials. The crucial fact, upon which all of modern signal processing is based, is that the sampled complex exponentials form an orthogonal basis. The section introduces the Discrete Fourier Transform, and concludes with an introduction to the Fast Fourier Transform, an efficient algorithm for computing the discrete Fourier representation and reconstructing the signal from its Fourier coefficients. We will concentrate on the one-dimensional version here. Let f (x) be a function representing the signal, defined on an interval a x b. Our computer can only store its measured values at a finite number of sample points a x0 < x1 < < xn b. In the simplest and, by far, the most common case, the sample points are equally spaced, and so b-a n indicates the sample rate. In signal processing applications, x represents time instead of space, and the xj are the times at which we sample the signal f (x). Sample rates can be very high, e.g., every 1020 milliseconds in current speech recognition systems. For simplicity, we adopt the "standard" interval of 0 x 2 , and the n equally spaced sample points xj = a + j h, j = 0, . . . , n, where h= 4 2j 2(n - 1) 2 , x2 = , . . . xj = , . . . xn-1 = . (13.1) n n n n (Signals defined on other intervals can be handled by simply rescaling the interval to have length 2 .) Sampling a (complex-valued) signal or function f (x) produces the sample x0 = 0, x1 =
We will find it convenient to omit the final sample point xn = 2 from consideration.
2/25/07
692
c 2006
Peter J. Olver
1 0.5
1 0.5
1 0.5
1 -0.5 -1
2
3
4
5
6 -0.5 -1
1
2
3
4
5
6 -0.5 -1
1
2
3
4
5
6
1 0.5
1 0.5
1 0.5
1 -0.5 -1
2
3
4
5
6 -0.5 -1
1
2
3
4
5
6 -0.5 -1
1
2
3
4
5
6
Figure 13.1. vector
Sampling e- i x and e7 i x on n = 8 sample points.
f = f0 , f1 , . . . , fn-1 where
T
= f (x0 ), f (x1 ), . . . , f (xn-1 ) 2j n .
T
,
fj = f (xj ) = f
(13.2)
Sampling cannot distinguish between functions that have the same values at all of the sample points -- from the sampler's point of view they are identical. For example, the periodic complex exponential function f (x) = e i n x = cos n x + i sin n x has sampled values fj = f 2j n = exp in 2j n = e2 j i = 1 for all j = 0, . . . , n - 1,
and hence is indistinguishable from the constant function c(x) 1 -- both lead to the T same sample vector ( 1, 1, . . . , 1 ) . This has the important implication that sampling at n equally spaced sample points cannot detect periodic signals of frequency n. More generally, the two complex exponential signals e i (k+n)x and eikx
are also indistinguishable when sampled. This has the important consequence that we need only use the first n periodic complex exponential functions f0 (x) = 1, f1 (x) = e i x , f2 (x) = e2 i x , ... fn-1 (x) = e(n-1) i x , (13.3)
in order to represent any 2 periodic sampled signal. In particular, exponentials e- i k x of "negative" frequency can all be converted into positive versions, namely e i (n-k)x , by the same sampling argument. For example, e- i x = cos x - i sin x 2/25/07 and e(n-1) i x = cos(n - 1) x + i sin(n - 1) x 693
c 2006 Peter J. Olver
have identical values on the sample points (13.1). However, off of the sample points, they are quite different; the former is slowly varying, while the latter represents a high frequency oscillation. In Figure 13.1, we compare e- i x and e7 i x when there are n = 8 sample values, indicated by the dots on the graphs. The top row compares the real parts, cos x and cos 7 x, while the bottom row compares the imaginary parts, sin x and - sin 7 x. Note that both functions have the same pattern of sample values, even though their overall behavior is strikingly different. This effect is commonly referred to as aliasing . If you view a moving particle under a stroboscopic light that flashes only eight times, you would be unable to determine which of the two graphs the particle was following. Aliasing is the cause of a well-known artifact in movies: spoked wheels can appear to be rotating backwards when our brain interprets the discretization of the high frequency forward motion imposed by the frames of the film as an equivalently discretized low frequency motion in reverse. Aliasing also has important implications for the design of music CD's. We must sample an audio signal at a sufficiently high rate that all audible frequencies can be adequately represented. In fact, human appreciation of music also relies on inaudible high frequency tones, and so a much higher sample rate is actually used in commercial CD design. But the sample rate that was selected remains controversial; hi fi aficionados complain that it was not set high enough to fully reproduce the musical quality of an analog LP record! The discrete Fourier representation decomposes a sampled function f (x) into a linear combination of complex exponentials. Since we cannot distinguish sampled exponentials of frequency higher than n, we only need consider a finite linear combination
n-1
f (x) p(x) = c0 + c1 e i x + c2 e2 i x + + cn-1 e(n-1) i x =
ck e i k x
k=0
(13.4)
of the first n exponentials (13.3). The symbol in (13.4) means that the function f (x) and the sum p(x) agree on the sample points: f (xj ) = p(xj ), j = 0, . . . , n - 1. (13.5)
Therefore, p(x) can be viewed as a (complex-valued) interpolating trigonometric polynomial of degree n - 1 for the sample data fj = f (xj ). Remark : If f (x) is real, then p(x) is also real on the sample points, but may very well be complex-valued in between. To avoid this unsatisfying state of affairs, we will usually discard its imaginary component, and regard the real part of p(x) as "the" interpolating trigonometric polynomial. On the other hand, sticking with a purely real construction unnecessarily complicates the analysis, and so we will retain the complex exponential form (13.4) of the discrete Fourier sum.
In computer graphics, the term "aliasing" is used in a much broader sense that covers a variety of artifacts introduced by discretization -- particularly, the jagged appearance of lines and smooth curves on a digital monitor.
2/25/07
694
c 2006
Peter J. Olver
Since we are working in the finite-dimensional vector space C n throughout, we may reformulate the discrete Fourier series in vectorial form. Sampling the basic exponentials (13.3) produces the complex vectors k = e i k x0 , e i k x1 , e i k x2 , . . . , e i k xn =
T T
1, e2 k i /n , e4 k i /n , . . . , e2 (n-1) k i /n
,
k = 0, . . . , n - 1.
(13.6)
The interpolation conditions (13.5) can be recast in the equivalent vector form f = c0 0 + c1 1 + + cn-1 n-1 . (13.7)
In other words, to compute the discrete Fourier coefficients c0 , . . . , cn-1 of f , all we need to do is rewrite its sample vector f as a linear combination of the sampled exponential vectors 0 , . . . , n-1 . Now, as with continuous Fourier series, the absolutely crucial property is the orthonormality of the basis elements 0 , . . . , n-1 . Were it not for the power of orthogonality, Fourier analysis might have remained a mere mathematical curiosity, rather than today's indispensable tool. Proposition 13.1. The sampled exponential vectors 0 , . . . , n-1 form an orthonormal basis of C n with respect to the inner product f ,g 1 = n
n-1
fj gj
j =0
1 = n
n-1
f (xj ) g(xj ) ,
j =0
f , g Cn.
(13.8)
Remark : The inner product (13.8) is a rescaled version of the standard Hermitian dot product (3.90) between complex vectors. We can interpret the inner product between the sample vectors f , g as the average of the sampled values of the product signal f (x) g(x). Remark : As usual, orthogonality is no accident. Just as the complex exponentials are eigenfunctions for a self-adjoint boundary value problem, so their discrete sampled counterparts are eigenvectors for a self-adjoint matrix eigenvalue problem; details can be found in Exercise 8.4.12. Here, though, to keep the discussion on track, we shall outline a direct proof. Proof : The crux of the matter relies on properties of the remarkable complex numbers n = e2 i /n = cos Particular cases include 2 = -1,
1 3 = - 2 + 3 2
2 2 + i sin , n n
where
n = 1, 2, 3, . . . .
(13.9)
i,
4 = i ,
n
and
8 =
2 2
+
2 2
i.
(13.10)
The nth power of n is
n n = e2 i /n
= e2 i = 1,
c 2006 Peter J. Olver
2/25/07
695
5
2 5 0 5 = 1 3 5 4 5
Figure 13.2.
The Fifth Roots of Unity.
and hence n is one of the complex nth roots of unity: n = n 1. There are, in fact, n different complex nth roots of 1, including 1 itself, namely the powers of n :
k n = e2 k i /n = cos
2k 2k + i sin , n n
k = 0, . . . , n - 1.
(13.11)
Since it generates all the others, n is known as the primitive nth root of unity. Geometrically, the nth roots (13.11) lie on the vertices of a regular unit ngon in the complex plane; see Figure 13.2. The primitive root n is the first vertex we encounter as we go around the ngon in a counterclockwise direction, starting at 1. Continuing around, the other roots 2 3 n-1 n appear in their natural order n , n , . . . , n , and finishing back at n = 1. The complex conjugate of n is the "last" nth root e-2 i /n = n = 1 n-1 = n = e2(n-1) i /n . n (13.12)
The complex numbers (13.11) are a complete set of roots of the polynomial z n - 1, which can therefore be factored:
2 n-1 z n - 1 = (z - 1)(z - n )(z - n ) (z - n ).
On the other hand, elementary algebra provides us with the real factorization z n - 1 = (z - 1)(1 + z + z 2 + + z n-1 ). Comparing the two, we conclude that
2 n-1 1 + z + z 2 + + z n-1 = (z - n )(z - n ) (z - n ). k Substituting z = n into both sides of this identity, we deduce the useful formula (n-1) k k 2 = 1 + n + n k + + n
n, 0,
k = 0, 0 < k < n.
(13.13)
n+k k Since n = n , this formula can easily be extended to general integers k; the sum is equal to n if n evenly divides k and is 0 otherwise.
2/25/07
696
c 2006
Peter J. Olver
Now, let us apply what we've learned to prove Proposition 13.1. First, in view of (13.11), the sampled exponential vectors (13.6) can all be written in terms of the nth roots of unity:
k 2 3 (n-1) k k = 1, n , n k , n k , . . . , n T
,
k = 0, . . . , n - 1. 1, 0, k = l, k = l,
(13.14)
Therefore, applying (13.12, 13), we conclude that k , l = 1 n
n-1 jl j n k n = j =0
1 n
n-1 j(k-l) n = j =0
0 k, l < n, Q.E.D.
which establishes orthonormality of the sampled exponential vectors.
Orthonormality of the basis vectors implies that we can immediately compute the Fourier coefficients in the discrete Fourier sum (13.4) by taking inner products: ck = f , k 1 = n
n-1
fj
j =0
e i k xj
1 = n
n-1
fj e
j =0
- i k xj
1 = n
n-1 - n j k fj . j =0
(13.15)
In other words, the discrete Fourier coefficient ck is obtained by averaging the sampled values of the product function f (x) e- i k x . The passage from a signal to its Fourier coefficients is known as the Discrete Fourier Transform or DFT for short. The reverse procedure of reconstructing a signal from its discrete Fourier coefficients via the sum (13.4) (or (13.7)) is known as the Inverse Discrete Fourier Transform or IDFT. The Discrete Fourier Transform and its inverse define mutually inverse linear transformations on the space C n , whose matrix representations can be found in Exercise 5.7.9. Example 13.2. If n = 4, then 4 = i . The corresponding sampled exponential vectors 1 1 1 1 1 i -1 -i 0 = , 1 = 2 = 3 = , , , 1 -1 1 -1 1 -i -1 i form an orthonormal basis of C 4 with respect to the averaged Hermitian dot product w0 v0 v1 w where v = , w = 1 . v , w = 1 v0 w0 + v1 w1 + v2 w2 + v3 w3 , 4 w2 v2 w3 v3 Given the sampled function values f0 = f (0), f1 = f
1 2
,
f2 = f (),
f3 = f
3 2
,
we construct the discrete Fourier representation f = c0 0 + c1 1 + c2 2 + c3 3 , where
1 c0 = f , 0 = 4 (f0 + f1 + f2 + f3 ), c2 = f , 2 = 1 (f0 - f1 + f2 - f3 ), 4
(13.16)
c1 = f , 1 = 1 (f0 - i f1 - f2 + i f3 ), 4 1 c3 = f , 3 = 4 (f0 + i f1 - f2 - i f3 ). 697
c 2006 Peter J. Olver
2/25/07
10 8 6 4 2 1 10 8 6 4 2 1 2 3 4 5 6 2 3 4 5 6
10 8 6 4 2 1 10 8 6 4 2 1 2 3 4 5 6 2 3 4 5 6
10 8 6 4 2 1 10 8 6 4 2 1 2 3 4 5 6 2 3 4 5 6
Figure 13.3.
The Discrete Fourier Representation of x2 - 2 x.
We interpret this decomposition as the complex exponential interpolant f (x) p(x) = c0 + c1 e i x + c2 e2 i x + c3 e3 i x that agrees with f (x) on the sample points. For instance, if f (x) = 2 x - x2 , then f0 = 0., and hence c0 = 6.1685, c1 = -2.4674, c2 = -1.2337, c3 = -2.4674. (13.17) f1 = 7.4022, f2 = 9.8696, f3 = 7.4022,
Therefore, the interpolating trigonometric polynomial is given by the real part of p(x) = 6.1685 - 2.4674 e i x - 1.2337 e2 i x - 2.4674 e3 i x , namely, Re p(x) = 6.1685 - 2.4674 cos x - 1.2337 cos 2 x - 2.4674 cos 3 x. (13.18)
In Figure 13.3 we compare the function, with the interpolation points indicated, and discrete Fourier representations (13.18) for both n = 4 and n = 16 points. The resulting graphs point out a significant difficulty with the Discrete Fourier Transform as developed so far. While the trigonometric polynomials do indeed correctly match the sampled function values, their pronounced oscillatory behavior makes them completely unsuitable for interpolation away from the sample points. However, this difficulty can be rectified by being a little more clever. The problem is that we have not been paying sufficient attention to the frequencies that are represented in the Fourier sum. Indeed, the graphs in Figure 13.3 might remind you of our earlier observation that, due to aliasing, low and high frequency exponentials can have the same sample data, but differ wildly in between the sample points. While the first half of the 2/25/07 698
c 2006 Peter J. Olver
10 8 6 4 2 1 10 8 6 4 2 1 2 3 4 5 6 2 3 4 5 6
10 8 6 4 2 1 10 8 6 4 2 1 2 3 4 5 6 2 3 4 5 6
10 8 6 4 2 1 10 8 6 4 2 1 2 3 4 5 6 2 3 4 5 6
Figure 13.4.
The Low Frequency Discrete Fourier Representation of x2 - 2 x.
summands in (13.4) represent relatively low frequencies, the second half do not, and can be replaced by equivalent lower frequency, and hence less oscillatory exponentials. Namely, if 0 < k 1 n, then e- i k x and e i (n-k) x have the same sample values, but the former is of 2 lower frequency than the latter. Thus, for interpolatory purposes, we should replace the second half of the summands in the Fourier sum (13.4) by their low frequency alternatives. If n = 2 m + 1 is odd, then we take
m
p(x) = c- m e
- i mx
+ +c-1 e
-ix
+c0 +c1 e
ix
+ +cm e
i mx
=
k=-m
ck e i k x (13.19)
as the equivalent low frequency interpolant. If n = 2 m is even -- which is the most common case occurring in applications -- then
m-1
p(x) = c- m e
- i mx
+ +c-1 e
-ix
+c0 +c1 e
ix
+ +cm-1 e
i (m-1) x
=
k=-m - i mx
ck e i k x
(13.20) will be our choice. (It is a matter of personal taste whether to use e or e i m x to represent the highest frequency term.) In both cases, the Fourier coefficients with negative indices are the same as their high frequency alternatives: c- k = cn-k = f , n-k = f , - k , (13.21)
where - k = n-k is the sample vector for e- i k x e i (n-k) x . Returning to the previous example, for interpolating purposes, we should replace (13.17) by the equivalent low frequency interpolant p (x) = - 1.2337 e- 2 i x - 2.4674 e- i x + 6.1685 - 2.4674 e i x , with real part Re p (x) = 6.1685 - 4.9348 cos x - 1.2337 cos 2 x. Graphs of the n = 4 and 16 low frequency trigonometric interpolants can be seen in Figure 13.4. Thus, by utilizing only the lowest frequency exponentials, we successfully 2/25/07 699
c 2006 Peter J. Olver
(13.22)
suppress the aliasing artifacts, resulting in a quite reasonable trigonometric interpolant to the given function. Remark : The low frequency version also serves to unravel the reality of the Fourier representation of a real function f (x). Since - k = k , formula (13.21) implies that c- k = ck , and so the common frequency terms c- k e- i k x + ck e i k x = ak cos k x + bk sin k x add up to a real trigonometric function. Therefore, the odd n interpolant (13.19) is a real trigonometric polynomial, whereas in the even version (13.20) only the highest frequency term c- m e- i m x produces a complex term -- which is, in fact, 0 on the sample points. Compression and Noise Removal In a typical experimental signal, noise primarily affects the high frequency modes, while the authentic features tend to appear in the low frequencies. Think of the hiss and static you hear on an AM radio station or a low quality audio recording. Thus, a very simple, but effective, method for denoising a corrupted signal is to decompose it into its Fourier modes, as in (13.4), and then discard the high frequency constituents. A similar idea underlies the Dolby T recording system used on most movie soundtracks: during the recording process, the high frequency modes are artificially boosted, so that scaling them back when showing the movie in the theater has the effect of eliminating much of the extraneous noise. The one design issue is the specification of a cut-off between low and high frequency, that is, between signal and noise. This choice will depend upon the properties of the measured signal, and is left to the discretion of the signal processor. A correct implementation of the denoising procedure is facilitated by using the unaliased forms (13.19, 20) of the trigonometric interpolant, in which the low frequency summands only appear when | k | is small. In this version, to eliminate high frequency components, we replace the full summation by
l
ql (x) =
1 2 k=-l
ck e i k x ,
(13.23)
where l < (n + 1) specifies the selected cut-off frequency between signal and noise. The 2 l + 1 n low frequency Fourier modes retained in (13.23) will, in favorable situations, capture the essential features of the original signal while simultaneously eliminating the high frequency noise. In Figure 13.5 we display a sample signal followed by the same signal corrupted by adding in random noise. We use n = 28 = 256 sample points in the discrete Fourier representation, and to remove the noise, we retain only the 2 l + 1 = 11 lowest frequency modes. In other words, instead of all n = 512 Fourier coefficients c-256 , . . . , c-1 , c0 , c1 , . . . , c255 , we only compute the 11 lowest order ones c-5 , . . . , c5 . Summing up just those 11 exponentials produces the denoised signal q(x) = c-5 e- 5 i x + + c5 e5 i x . To compare, we plot both the original signal and the denoised version on the same graph. In this case, the maximal deviation is less than .15 over the entire interval [ 0, 2 ]. 2/25/07 700
c 2006 Peter J. Olver
7 6 5 4 3 2 1 1 2 3 4 5 6
7 6 5 4 3 2 1 1 2 3 4 5 6
The Original Signal
7 6 5 4 3 2 1 1 2 3 4 5 6 7 6 5 4 3 2 1
The Noisy Signal
1
2
3
4
5
6
The Denoised Signal Figure 13.5.
7 6 5 4 3 2 1 1 2 3 4 5 6 7 6 5 4 3 2 1 1 2 3 4
Comparison of the Two Denoising a Signal.
7 6 5 4 3 2 1 5 6 1 2 3 4 5 6
The Original Signal
Moderate Compression Compressing a Signal.
High Compression
Figure 13.6.
The same idea underlies many data compression algorithms for audio recordings, digital images and, particularly, video. The goal is efficient storage and/or transmission of the signal. As before, we expect all the important features to be contained in the low frequency constituents, and so discarding the high frequency terms will, in favorable situations, not lead to any noticeable degradation of the signal or image. Thus, to compress a signal (and, simultaneously, remove high frequency noise), we retain only its low frequency discrete Fourier coefficients. The signal is reconstructed by summing the associated truncated discrete Fourier series (13.23). A mathematical justification of Fourier-based compression algorithms relies on the fact that the Fourier coefficients of smooth functions tend rapidly to zero -- the smoother the function, the faster the decay rate. Thus, the small high frequency Fourier coefficients will be of negligible importance. In Figure 13.6, the same signal is compressed by retaining, respectively, 2 l + 1 = 21 and 2 l + 1 = 7 Fourier coefficients only instead of all n = 512 that would be required for complete accuracy. For the case of moderate compression, the maximal deviation between the signal and the compressed version is less than 1.5 10-4 over the entire interval, 2/25/07 701
c 2006 Peter J. Olver
while even the highly compressed version deviates at most .05 from the original signal. Of course, the lack of any fine scale features in this particular signal means that a very high compression can be achieved -- the more complicated or detailed the original signal, the more Fourier modes need to be retained for accurate reproduction. The Fast Fourier Transform While one may admire an algorithm for its intrinsic beauty, in the real world, the bottom line is always efficiency of implementation: the less total computation, the faster the processing, and hence the more extensive the range of applications. Orthogonality is the first and most important feature of many practical linear algebra algorithms, and is the critical feature of Fourier analysis. Still, even the power of orthogonality reaches its limits when it comes to dealing with truly large scale problems such as three-dimensional medical imaging or video processing. In the early 1960's, James Cooley and John Tukey, [44], discovered a much more efficient approach to the Discrete Fourier Transform, exploiting the rather special structure of the sampled exponential vectors. The resulting algorithm is known as the Fast Fourier Transform, often abbreviated FFT, and its discovery launched the modern revolution in digital signal and data processing, [30, 31]. In general, computing all the discrete Fourier coefficients (13.15) of an n times sampled signal requires a total of n2 complex multiplications and n2 - n complex additions. Note also that each complex addition z + w = (x + i y) + (u + i v) = (x + u) + i (y + v) generally requires two real additions, while each complex multiplication z w = (x + i y) (u + i v) = (x u - y v) + i (x v + y u) x v + y u = (x + y) (u + v) - x u - y v (13.25) (13.24)
requires 4 real multiplications and 2 real additions, or, by employing the alternative formula (13.26)
for the imaginary part, 3 real multiplications and 5 real additions. (The choice of formula (13.25) or (13.26) will depend upon the processor's relative speeds of multiplication and addition.) Similarly, given the Fourier coefficients c0 , . . . , cn-1 , reconstruction of the sampled signal via (13.4) requires n2 -n complex multiplications and n2 -n complex additions. As a result, both computations become quite labor intensive for large n. Extending these ideas to multi-dimensional data only exacerbates the problem. In order to explain the method without undue complication, we return to the original, aliased form of the discrete Fourier representation (13.4). (Once one understands how the FFT works, one can easily adapt the algorithm to the low frequency version (13.20).) The seminal observation is that if the number of sample points n = 2m
In fact, the key ideas can be found in Gauss' hand computations in the early 1800's, but his insight was not fully appreciated until modern computers arrived on the scene.
2/25/07
702
c 2006
Peter J. Olver
is even, then the primitive mth root of unity m = nth root: 2 m = n .
m
1 equals the square of the primitive
We use this fact to split the summation (13.15) for the order n discrete Fourier coefficients k into two parts, collecting together the even and the odd powers of n : ck = 1 -k - - f + f1 n + f2 n 2 k + + fn-1 n (n-1) k n 0 1 - - - = f0 + f2 n 2 k + f4 n 4 k + + f2m-2 n (2 m-2) k + n -k 1 - - - + n f + f3 n 2 k + f5 n 4 k + + f2m-1 n (2 m-2) k n 1 1 1 - - - + f + f2 m k + f4 m 2 k + + f2m-2 m (m-1) k = 2 m 0 +
- n k 2
(13.27)
1 - - - f + f3 m k + f5 m 2 k + + f2m-1 m (m-1) k m 1
.
Now, observe that the expressions in braces are the order m Fourier coefficients for the sample data f e = f0 , f2 , f4 , . . . , f2m-2 f o = f1 , f3 , f5 , . . . , f2m-1
T T
= f (x0 ), f (x2 ), f (x4 ), . . . , f (x2m-2 ) = f (x1 ), f (x3 ), f (x5 ), . . . , f (x2m-1 )
T T
, .
(13.28)
Note that f e is obtained by sampling f (x) on the even sample points x2j , while f o is obtained by sampling the same function f (x), but now at the odd sample points x2j+1 . In other words, we are splitting the original sampled signal into two "half-sampled" signals obtained by sampling on every other point. The even and odd Fourier coefficients are 1 - - - ce = f0 + f2 m k + f4 m 2 k + + f2m-2 m (m-1) k , k m k = 0, . . . , m - 1. 1 -k -2k - (m-1) k o f + f3 m + f5 m + + f2m-1 m , ck = (13.29) m 1 Since they contain just m data values, both the even and odd samples require only m distinct Fourier coefficients, and we adopt the identification
e ce k+m = ck , o co k+m = ck ,
k = 0, . . . , m - 1.
(13.30)
Therefore, the order n = 2 m discrete Fourier coefficients (13.27) can be constructed from a pair of order m discrete Fourier coefficients via ck =
1 2 -k ce + n co , k k
k = 0, . . . , n - 1.
(13.31)
Now if m = 2 l is also even, then we can play the same game on the order m Fourier coefficients (13.29), reconstructing each of them from a pair of order l discrete Fourier coefficients -- obtained by sampling the signal at every fourth point. If n = 2r is a power of 2, then this game can be played all the way back to the start, beginning with the trivial order 1 discrete Fourier representation, which just samples the function at a single point. 2/25/07 703
c 2006 Peter J. Olver
The result is the desired algorithm. After some rearrangement of the basic steps, we arrive at the Fast Fourier Transform, which we now present in its final form. We begin with a sampled signal on n = 2r sample points. To efficiently program the Fast Fourier Transform, it helps to write out each index 0 j < 2r in its binary (as opposed to decimal) representation j = jr-1 jr-2 . . . j2 j1 j0 , where j = 0 or 1; (13.32)
the notation is shorthand for its r digit binary expansion j = j0 + 2 j1 + 4 j2 + 8 j3 + + 2r-1 jr-1 . We then define the bit reversal map (jr-1 jr-2 . . . j2 j1 j0 ) = j0 j1 j2 . . . jr-2 jr-1 . (13.33)
For instance, if r = 5, and j = 13, with 5 digit binary representation 01101, then (j) = 22 has the reversed binary representation 10110. Note especially that the bit reversal map = r depends upon the original choice of r = log2 n. Secondly, for each 0 k < r, define the maps k (j) = jr-1 . . . jk+1 0 jk-1 . . . j0 , k (j) = jr-1 . . . jk+1 1 jk-1 . . . j0 = k (j) + 2k , k th for j = jr-1 jr-2 . . . j1 j0 .
(13.34) In other words, k (j) sets the binary digit of j to 0, while k (j) sets it to 1. In the preceding example, 2 (13) = 9, with binary form 01001, while 2 (13) = 13 with binary form 01101. The bit operations (13.33, 34) are especially easy to implement on modern binary computers. Given a sampled signal f0 , . . . , fn-1 , its discrete Fourier coefficients c0 , . . . , cn-1 are computed by the following iterative algorithm: cj
(0)
= f(j) ,
cj
(k+1)
=
1 2
-j ck (j) + 2k+1 ck (j) ,
(k)
(k)
j = 0, . . . , n - 1, k = 0, . . . , r - 1,
(13.35)
in which 2k+1 is the primitive 2k+1 root of unity. The final output of the iterative procedure, namely (r) (13.36) j = 0, . . . , n - 1, cj = cj , are the discrete Fourier coefficients of our signal. The preprocessing step of the algorithm, (0) where we define cj , produces a more convenient rearrangement of the sample values. The subsequent steps successively combine the Fourier coefficients of the appropriate even and odd sampled subsignals together, reproducing (13.27) in a different notation. The following example should help make the overall process clearer. Example 13.3. Consider the case r = 3, and so our signal has n = 23 = 8 sampled values f0 , f1 , . . . , f7 . We begin the process by rearranging the sample values c0 = f0 , c1 = f4 , c2 = f2 , c3 = f6 , c4 = f1 , c5 = f5 , c6 = f3 , c7 = f7 , 2/25/07 704
c 2006 Peter J. Olver (0) (0) (0) (0) (0) (0) (0) (0)
in the order specified by the bit reversal map . For instance (3) = 6, or, in binary notation, (011) = 110. The first stage of the iteration is based on 2 = -1. Equation (13.35) gives
1 c0 = 1 (c0 + c1 ), c1 = 1 (c0 - c1 ), c2 = 2 (c2 + c3 ), c3 = 1 (c2 - c3 ), 2 2 2 1 1 c4 = 2 (c4 + c5 ), c5 = 1 (c4 - c5 ), c6 = 2 (c6 + c7 ), c7 = 1 (c6 - c7 ), 2 2 (1) (0) (0) (1) (0) (0) (1) (0) (0) (1) (0) (0) (1) (0) (0) (1) (0) (0) (1) (0) (0) (1) (0) (0)
where we combine successive pairs of the rearranged sample values. The second stage of the iteration has k = 1 with 4 = i . We find
1 c0 = 1 (c0 + c2 ), c1 = 2 (c1 - i c3 ), c2 = 1 (c0 - c2 ), c3 = 1 (c1 + i c3 ), 2 2 2 (2) (1) (1) (2) (1) (1) (2) (1) (1) (2) (1) (1) (2) (1) (1) (2) (1) (1) (2) (1) (1) (2) (1) (1)
1 c4 = 1 (c4 + c6 ), c5 = 2 (c5 - i c7 ), c6 = 1 (c4 - c6 ), c7 = 1 (c5 + i c7 ). 2 2 2
Note that the indices of the combined pairs of coefficients differ by 2. In the last step, where k = 2 and 8 = 22 (1 + i ), we combine coefficients whose indices differ by 4 = 22 ; the final output
1 c0 = c0 = 2 (c0 + c4 ), (3) (2) (2)
c1 = c1 = c2 = c2 = c3 = c3 =
(3) (3)
(3)
1 2 1 2 1 2
c1 + c2 - c3 -
(2) (2)
(2)
2 2 (1 (2) i c6 2 2 (1
- i ) c5 , + i ) c7
(2)
c4 = c4 = 1 (c0 - c4 ), 2 , c5 = c5 = c6 = c6 =
(3) (3) 1 2 1 2 1 2
(3)
(2)
(2)
c1 - c2 + c3 +
(2) (2)
(2)
(2)
,
c7 = c7 =
(3)
2 2 (1 (2) i c6 2 2 (1
- i ) c5 , + i ) c7
(2)
,
(2)
,
is the complete set of discrete Fourier coefficients. Let us count the number of arithmetic operations required in the Fast Fourier Transform algorithm. At each stage in the computation, we must perform n = 2r complex additions/subtractions and the same number of complex multiplications. (Actually, the number of multiplications is slightly smaller since multiplications by 1 and i are extremely simple. However, this does not significantly alter the final operations count.) There are r = log2 n stages, and so we require a total of r n = n log2 n complex additions/subtractions and the same number of multiplications. Now, when n is large, n log2 n is significantly smaller than n2 , which is the number of operations required for the direct algorithm. For instance, if n = 210 = 1, 024, then n2 = 1, 048, 576, while n log2 n = 10, 240 -- a net savings of 99%. As a result, many large scale computations that would be intractable using the direct approach are immediately brought into the realm of feasibility. This is the reason why all modern implementations of the Discrete Fourier Transform are based on the FFT algorithm and its variants. The reconstruction of the signal from the discrete Fourier coefficients c0 , . . . , cn-1 is -1 speeded up in exactly the same manner. The only differences are that we replace n = n 1 by n , and drop the factors of 2 since there is no need to divide by n in the final result (13.4). Therefore, we apply the slightly modified iterative procedure fj
(0)
= c(j) ,
fj
(k+1)
j = fk (j) + 2k+1 fk (j) ,
(k)
(k)
j = 0, . . . , n - 1, k = 0, . . . , r - 1,
c 2006
(13.37)
Peter J. Olver
2/25/07
705
and finish with
f (xj ) = fj = fj ,
(r)
j = 0, . . . , n - 1.
(13.38)
Example 13.4. The reconstruction formulae in the case of n = 8 = 23 Fourier coefficients c0 , . . . , c7 , which were computed in Example 13.3, can be implemented as follows. First, we rearrange the Fourier coefficients in bit reversed order: f0
(0)
= c0 , f1 f0
(1) (1)
(0)
= c4 , f2 + f1 , + f5 ,
(1) (0) (0)
(0)
= c2 , f3
(1) (1)
(0)
= c6 , f4
(0)
(0)
= c1 , f5 = f2
(0) (0)
(0)
= c5 , f6
(0)
(0)
= c3 , f7
(0) (0)
(0)
= c7 ,
Then we begin combining them in successive pairs: = f0
(0) (0)
f1
= f0
(0) (0)
- f1 , f2
(0)
(1) (1)
+ f3 , + f7 , - f2 , - f6 ,
(3) (1) (1) (0)
f3
(1) (1)
= f2
- f3 , - f7 . - i f3 , - i f7 .
(1) (1) (0)
(0)
f4 Next, f0
= f4
f5
= f4
- f5 , f6
(1)
= f6
f7
= f6
(2) (2)
= f0
(1) (1)
+ f2 , f1
(1)
(2) (2)
= f1
(1) (1)
+ i f3 , f2
(1)
(2) (2)
= f0
(1) (1)
f3
(2) (2)
= f1
(1) (1)
f4
= f4
+ f6 , f5 = f0 = f2
(2) (2)
= f5
+ i f7 , f6
= f4
f7
= f5
Finally, the sampled signal values are f (x0 ) = f0 f (x2 ) = f2
(3) (3)
+ f4 , +
2 2
(2)
f (x4 ) = f4
(2)
= f0 = f1 = f2 = f3
(2)
f (x1 ) = f1 f (x3 ) = f3
= f1 = f3
(1 + i ) f5 ,
(2)
f (x5 ) = f5 f (x7 ) = f7
(3) (3)
(2) (2)
- f4 , - +
2 2
(2)
(1 + i ) f5 ,
(2)
(2)
(3) (3)
(2) (2)
+ i f6 , -
2 2
f (x6 ) = f6
(2)
(1 - i ) f7 ,
(3)
(2)
- i f6 ,
2 2
(1 - i ) f7 .
(2)
13.2. Wavelets.
Trigonometric Fourier series, both continuous and discrete, are amazingly powerful, but they do suffer from one potentially serious defect. The basis functions e i k x = cos k x + i sin k x are spread out over the entire interval [ - , ], and so are not well-suited to processing localized signals -- meaning data that are concentrated in a relatively small regions. Indeed, the most concentrated data of all -- a single delta function -- has every Fourier component of equal magnitude in its Fourier series (12.61) and its high degree of localization is completely obscured. Ideally, one would like to construct a system of functions that is orthogonal, and so has all the advantages of the Fourier trigonometric functions, but, in addition, adapts to localized structures in signals. This dream was the inspiration for the development of the modern theory of wavelets. The Haar Wavelets Although the modern era of wavelets started in the mid 1980's, the simplest example of a wavelet basis was discovered by the Hungarian mathematician Alfrd Haar in 1910, e [88]. We consider the space of functions (signals) defined the interval [ 0, 1 ], equipped with the standard L2 inner product
1
f ,g =
0
f (x) g(x) dx.
c 2006
(13.39)
Peter J. Olver
2/25/07
706
1
1
0.5
0.5
-0.2
0.2
0.4
0.6
0.8
1
1.2
-0.2
0.2
0.4
0.6
0.8
1
1.2
-0.5
-0.5
-1
-1
1 (x)
1 1
2 (x)
0.5
0.5
-0.2
0.2
0.4
0.6
0.8
1
1.2
-0.2
0.2
0.4
0.6
0.8
1
1.2
-0.5
-0.5
-1
-1
3 (x) Figure 13.7.
4 (x) The First Four Haar Wavelets.
This choice is merely for convenience, being slightly better suited to our construction than [ - , ] or [ 0, 2 ]. Moreover, the usual scaling arguments can be used to adapt the wavelet formulas to any other interval. The Haar wavelets are certain piecewise constant functions. The first four are graphed in Figure 13.7 The first is the box function 1 (x) = (x) = 1, 0, 0 < x 1, otherwise, (13.40)
is known as the mother wavelet. The third and fourth Haar functions are compressed 2/25/07 707
c 2006 Peter J. Olver
known as the scaling function, for reasons that shall appear shortly. Although we are only interested in the value of (x) on the interval [ 0, 1 ], it will be convenient to extend it, and all the other wavelets, to be zero outside the basic interval. Its values at the points of discontinuity, i.e., 0, 1, is not critical, but, unlike the Fourier series midpoint value, it will be more convenient to consistently choose the left hand limiting value. The second Haar function 1 1, 0 < x 2, 1 (13.41) 2 (x) = w(x) = -1, 2 < x 1, 0, otherwise,
called daughter wavelets. One can easily check, by direct evaluation of the integrals, that the four Haar wavelet functions are orthogonal with respect to the L2 inner product (13.39). The scaling transformation x 2 x serves to compress the wavelet function, while the translation 2 x 2 x - 1 moves the compressed version to the right by a half a unit. Furthermore, we can represent the mother wavelet by compressing and translating the scaling function,: w(x) = (2 x) - (2 x - 1). (13.42) It is these two operations of scaling and compression -- coupled with the all-important orthogonality -- that underlies the power of wavelets. The Haar wavelets have an evident discretization. If we decompose the interval ( 0, 1 ] into the four subintervals
1 0, 4 , 1 4
versions of the mother wavelet: 0 < x 1, 1, 4 1 1 3 (x) = w(2 x) = -1, 4 < x 2, 0, otherwise,
4 (x) = w(2 x - 1) =
1,
1 2 3 4
< x 3, 4 < x 1,
-1, 0,
otherwise,
,1 , 2
1 2
3 ,4 ,
3 4
,1 ,
(13.43)
that form the orthogonal wavelet basis of R 4 we encountered in Examples 2.35 and 5.10. Orthogonality of the vectors (13.44) with respect to the standard Euclidean dot product is equivalent to orthogonality of the Haar wavelet functions with respect to the inner product (13.39). Indeed, if f (x) f = (f1 , f2 , f3 , f4 ) and g(x) g = (g1 , g2 , g3 , g4 )
on which the four wavelet functions are constant, then we can represent each of them by a vector in R 4 whose entries are the values of each wavelet function sampled at the left endpoint of each subinterval. In this manner, we obtain the wavelet sample vectors 1 1 1 0 1 1 -1 0 v1 = , v2 = v3 = v4 = (13.44) , , , 1 -1 0 1 1 -1 0 -1
are piecewise constant real functions that achieve the indicated values on the four subintervals (13.43), then their L2 inner product
1
f ,g =
0
f (x) g(x) dx =
1 4
f1 g1 + f2 g2 + f3 g3 + f4 g4 =
1 4
f g,
is equal to the averaged dot product of their sample values -- the real form of the inner product (13.8) that was used in the discrete Fourier transform. Since the vectors (13.44) form an orthogonal basis of R 4 , we can uniquely decompose any such piecewise constant function as a linear combination of wavelets f (x) = c1 1 (x) + c2 2 (x) + c3 3 (x) + c4 4 (x), 2/25/07 708
c 2006 Peter J. Olver
or, equivalently, in terms of the sample vectors, f = c1 v1 + c2 v2 + c3 v3 + c4 v4 . The required coefficients ck = f , k f vk = k 2 vk 2
1 2 1 2
are fixed by our usual orthogonality formula (5.7). Explicitly, c1 = c2 =
1 4 1 4
(f1 + f2 + f3 + f4 ), (f1 + f2 - f3 - f4 ),
c3 = c4 =
(f1 - f2 ), (f3 - f4 ).
Before proceeding to the more general case, let us introduce an important analytical definition that quantifies precisely how localized a function is. Definition 13.5. The support of a function f (x), written supp f , is the closure of the set where f (x) = 0. Thus, a point will belong to the support of f (x), provided f is not zero there, or at least is not zero at nearby points. More precisely: Lemma 13.6. If f (a) = 0, then a supp f . More generally, a point a supp f if and only if there exist a convergent sequence xn a such that f (xn ) = 0. Conversely, a supp f if and only if f (x) 0 on an interval a - < x < a + for some > 0. Intuitively, the smaller the support of a function, the more localized it is. For example, the support of the Haar mother wavelet (13.41) is supp w = [ 0, 1 ] -- the point x = 0 is included, even though w(0) = 0, because w(x) = 0 at nearby points. The two daughter wavelets have smaller support:
1 supp 3 = 0 , 2 ,
supp 4 =
1 2
,1 ,
and so are twice as localized. An extreme case is the delta function, whose support is a single point. In contrast, the support of the Fourier trigonometric basis functions is all of R, since they only vanish at isolated points. The effect of scalings and translations on the support of a function is easily discerned. Lemma 13.7. If supp f = [ a, b ], and g(x) = f (r x - ), then supp g = a+ b+ , . r r
In other words, scaling x by a factor r compresses the support of the function by a factor 1/r, while translating x translates the support of the function. The key requirement for a wavelet basis is that it contains functions with arbitrarily small support. To this end, the full Haar wavelet basis is obtained from the mother wavelet by iterating the scaling and translation processes. We begin with the scaling function (x), 2/25/07 709
c 2006
(13.45)
Peter J. Olver
from which we construct the mother wavelet via (13.42). For any "generation" j 0, we form the wavelet offspring by first compressing the mother wavelet so that its support fits into an interval of length 2-j , wj,0 (x) = w(2j x), so that supp wj,0 = [ 0, 2-j ], (13.46)
and then translating wj,0 so as to fill up the entire interval [ 0, 1 ] by 2j subintervals, each of length 2-j , defining wj,k (x) = wj,0 (x - k) = w(2j x - k), where k = 0, 1, . . . 2j - 1.
2j -1 [ k=0
(13.47)
Lemma 13.7 implies that supp wj,k = [ 2-j k, 2-j (k + 1) ], and so the combined supports of all the j th generation of wavelets is the entire interval: supp wj,k = [ 0, 1 ] . The
primal generation, j = 0, just consists of the mother wavelet w0,0 (x) = w(x).
The first generation, j = 1, consists of the two daughter wavelets already introduced as 3 and 4 , namely w1,0 (x) = w(2 x), w1,1 (x) = w(2 x - 1).
The second generation, j = 2, appends four additional granddaughter wavelets to our basis: w2,0 (x) = w(4 x), w2,1 (x) = w(4 x - 1), w2,2 (x) = w(4 x - 2), w2,3 (x) = w(4 x - 3).
The 8 Haar wavelets , w0,0 , w1,0 , w1,1 , w2,0 , w2,1 , w2,2 , w2,3 are constant on the 8 subintervals of length 1 , taking the successive sample values indicated by the columns of the 8 matrix 1 1 1 0 1 0 0 0 1 1 0 -1 0 0 0 1 1 -1 0 0 1 0 0 1 1 -1 0 0 -1 0 0 1 W8 = (13.48) . 0 1 0 0 1 0 1 -1 0 1 0 0 -1 0 1 -1 1 -1 0 -1 0 0 0 1 1 -1 0 -1 0 0 0 -1
Orthogonality of the wavelets is manifested in the orthogonality of the columns of W8 . (Unfortunately, the usual terminological constraints of Definition 5.18 prevent us from calling W8 an orthogonal matrix because its columns are not orthonormal!) The nth stage consists of 2n+1 different wavelet functions comprising the scaling functions and all the generations up to the nth : w0 (x) = (x) and wj,k (x) for 0 j n and 0 k < 2j . They are all constant on each subinterval of length 2-n-1 .
Theorem 13.8. The wavelet functions (x), wj,k (x) form an orthogonal system with respect to the inner product (13.39). 2/25/07 710
c 2006 Peter J. Olver
2
-j-1
Proof : First, note that each wavelet wj,k (x) is equal to + 1 on an interval of length and to - 1 on an adjacent interval of the same length. Therefore,
1
wj,k , =
0
wj,k (x) dx = 0,
(13.49)
since the + 1 and - 1 contributions cancel each other. If two different wavelets wj,k and wl,m with, say j l, have supports which are either disjoint, or just overlap at a single point, then their product wj,k (x) wl,m (x) 0, and so their inner product is clearly zero:
1
wj,k , wl,m =
0
wj,k (x) wl,m (x) dx = 0.
Otherwise, except in the case when the two wavelets are identical, the support of wl,m is entirely contained in an interval where wj,k is constant and so wj,k (x) wl,m (x) = wl,m (x). Therefore, by (13.49),
1 1
wj,k , wl,m = Finally, we compute
1 2
0
wj,k (x) wl,m (x) dx =
1
0
wl,m (x) dx = 0.
=
0
dx = 1,
wj,k
2
=
0
wj,k (x)2 dx = 2-j .
(13.50)
The second formula follows from the fact that | wj,k (x) | = 1 on an interval of length 2-j and is 0 elsewhere. Q.E.D. In direct analogy with the trigonometric Fourier series, the wavelet series of a signal f (x) is given by 2j -1 cj,k wj,k (x). (13.51) f (x) c0 (x) +
j =0 k=0
Orthogonality implies that the wavelet coefficients c0 , cj,k can be immediately computed using the standard inner product formula coupled with (13.50): c0 = cj,k f , = 2
1
f (x) dx,
0 2-j k+2-j-1 2-j k 2-j (k+1)
f , wj,k = 2j = wj,k 2
(13.52) f (x) dx.
f (x) dx - 2
j 2-j k+2-j-1
The convergence properties of the wavelet series (13.51) are similar to those of Fourier series; details can be found [52]. Example 13.9. In Figure 13.8, we plot the Haar expansions of the signal in the first plot. The next plots show the partial sums over j = 0, . . . , r with r = 2, 3, 4, 5, 6. We have used a discontinuous signal to demonstrate that there is no nonuniform Gibbs phenomenon in a Haar wavelet expansion. Indeed, since the wavelets are themselves discontinuous, they do not have any difficult uniformly converging to a discontinuous function. On the other hand, it takes quite a few wavelets to begin to accurately reproduce the signal. In the last plot, we combine a total of 26 = 64 Haar wavelets, which is considerably more than would be required in a comparably accurate Fourier expansion (excluding points very close to the discontinuity). 2/25/07 711
c 2006 Peter J. Olver
7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1
7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1
7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1
7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1
7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1
7 6 5 4 3 2 1 0.2 0.4 0.6 0.8 1
Figure 13.8.
Haar Wavelet Expansion.
Remark : To the novice, there may appear to be many more wavelets than trigonometric functions. But this is just another illusion of the magic show of infinite dimensional space. The point is that they both form a countably infinite set of functions, and so could, if necessary, but less conveniently, be numbered in order 1, 2, 3, . . . . On the other hand, accurate reproduction of functions usually does require many more Haar wavelets. This handicap makes the Haar system of less use in practical situations, and served to motivate the search for a more sophisticated choice of wavelet basis. Just as the discrete Fourier representation arises from a sampled version of the full Fourier series, so there is a discrete wavelet transformation for suitably sampled signals. To take full advantage of the wavelet basis, we sample the signal f (x) at n = 2r equally spaced sample points xk = k/2r , for k = 0, . . . , n - 1, on the interval [ 0, 1 ]. As before, we can identify the sampled signal with the vector f = f (x0 ), f (x1 ), . . . , f (xn-1 )
T
= f0 , f1 , . . . , fn-1
T
Rn.
(13.53)
Since we are only sampling on intervals of length 2-r , the discrete wavelet transform of our sampled signal will only use the first n = 2r wavelets (x) and wj,k (x) for j = 0, . . . , r - 1 and k = 0, . . . , 2j - 1. Let w0 (x) and wj,k wj,k (x) denote the corresponding sampled wavelet vectors; all of their entries are either +1, -1, or 0. (When r = 3, so n = 8, these are the columns of the wavelet matrix (13.48).) Moreover, orthogonality of the wavelets immediately implies that the wavelet vectors form an orthogonal basis of R n with n = 2r . We can decompose our sample vector (13.53) as a linear combination of the sampled wavelets, f = c0 w0 + 2/25/07
r-1 2j -1 j =0 k=0
cj,k wj,k ,
c 2006
(13.54)
712
Peter J. Olver
where, by our usual orthogonality formulae, f , w0 1 c0 = = r 2 w0 2 cj,k
2r -1 i=0
fi ,
k+2r-j-1 -1 i=k k+2r-j -1 i = k+2r-j-1
These are the basic formulae connecting the functions f (x), or, rather, its sample vector f , and its discrete wavelet transform consisting of the 2r coefficients c0 , cj,k . The reconstructed function f (x) = c0 (x) +
r-1 2j -1 j =0 k=0
f , wj,k = 2j-r = wj,k 2
fi -
fi .
(13.55)
cj,k wj,k (x)
(13.56)
is constant on each subinterval of length 2-r , and has the same value f (x) = f (xi ) = f (xi ) = fi , xi = 2-r i x < xi+1 = 2-r (i + 1),
as our signal at the left hand endpoint of the interval. In other words, we are interpolating the sample points by a piecewise constant (and thus discontinuous) function. Modern Wavelets The main defect of the Haar wavelets is that they do not provide a very efficient means of representing even very simple functions -- it takes quite a large number of wavelets to reproduce signals with any degree of precision. The reason for this is that the Haar wavelets are piecewise constant, and so even an affine function y = x + requires many sample values, and hence a relatively extensive collection of Haar wavelets, to be accurately reproduced. In particular, compression and denoising algorithms based on Haar wavelets are either insufficiently precise or hopelessly inefficient, and hence of minor practical value. For a long time it was thought that it was impossible to simultaneously achieve the requirements of localization, orthogonality and accurate reproduction of simple functions. The breakthrough came in 1988, when, in her Ph.D. thesis, the Dutch mathematician Ingrid Daubechies produced the first examples of wavelet bases that realized all three basic criteria. Since then, wavelets have developed into a sophisticated and burgeoning industry with major impact on modern technology. Significant applications include compression, storage and recognition of fingerprints in the FBI's data base, and the JPEG2000 image format, which, unlike earlier Fourier-based JPEG standards, incorporates wavelet technology in its image compression and reconstruction algorithms. In this section, we will present a brief outline of the basic ideas underlying Daubechies' remarkable construction. The recipe for any wavelet system involves two basic ingredients -- a scaling function and a mother wavelet. The latter can be constructed from the scaling function by a prescription similar to that in (13.42), and therefore we first concentrate on the properties of the scaling function. The key requirement is that the scaling function must solve a 2/25/07 713
c 2006 Peter J. Olver
1.2 1 0.8 0.6 0.4 0.2 -2 -1 -0.2 1 2 3 4
Figure 13.9. dilation equation of the form
p
The Hat Function.
(x) =
k=0
ck (2 x - k) = c0 (2 x) + c1 (2 x - 1) + + cp (2 x - p)
(13.57)
for some collection of constants c0 , . . . , cp . The dilation equation relates the function (x) to a finite linear combination of its compressed translates. The coefficients c0 , . . . , cp are not arbitrary, since the properties of orthogonality and localization will impose certain rather stringent requirements. Example 13.10. The Haar or box scaling function (13.40) satisfies the dilation equation (13.57) with c0 = c1 = 1, namely (x) = (2 x) + (2 x - 1). (13.58)
We recommend that you convince yourself of the validity of this identity before continuing. Example 13.11. Another example of a scaling function is the hat function 0 x 1, x, 2 - x, 1 x 2, (13.59) (x) = 0, otherwise, (x) = which is (13.57) with c0 = this identity by hand.
1 2
graphed in Figure 13.9, whose variants play a starring role in the finite element method, cf. (11.167). The hat function satisfies the dilation equation (2 x) + (2 x - 1) + = 1, c2 =
1 2. 1 2
(2 x - 2),
(13.60)
1 2 , c1
Again, the reader should be able to check
The dilation equation (13.57) is a kind of functional equation, and, as such, is not so easy to solve. Indeed, the mathematics of functional equations remains much less well developed than that of differential equations or integral equations. Even to prove that (nonzero) solutions exist is a nontrivial analytical problem. Since we already know two 2/25/07 714
c 2006 Peter J. Olver
explicit examples, let us defer the discussion of solution techniques until we understand how the dilation equation can be used to construct a wavelet basis. Given a solution to the dilation equation, we define the mother wavelet to be
p
w(x) =
k=0
(-1)k cp-k (2 x - k)
(13.61)
= cp (2 x) - cp-1 (2 x - 1) + cp-2 (2 x - 2) + c0 (2 x - p), This formula directly generalizes the Haar wavelet relation (13.42), in light of its dilation equation (13.58). The daughter wavelets are then all found, as in the Haar basis, by iteratively compressing and translating the mother wavelet: wj,k (x) = w(2j x - k). (13.62)
In the general framework, we do not necessarily restrict our attention to the interval [ 0, 1 ] and so j and k can, in principle, be arbitrary integers. Let us investigate what sort of conditions should be imposed on the dilation coefficients c0 , . . . , cp in order that we obtain a viable wavelet basis by this construction. First, localization of the wavelets requires that the scaling function has bounded support, and so (x) 0 when x lies outside some bounded interval [ a, b ]. If we integrate both sides of (13.57), we find
b
(x) dx =
a
-
p
(x) dx =
k=0
ck
-
(2 x - k) dx. dy, we find
b
(13.63)
Now using the change of variables y = 2 x - k, with dx =
-
1 2
(2 x - k) dx =
b
1 2
-
(y) dy =
1 2
(x) dx,
a
(13.64)
where we revert to x as our (dummy) integration variable. We substitute this result back into (13.63). Assuming that satisfy
a
(x) dx = 0, we discover that the dilation coefficients must c0 + + cp = 2. (13.65)
Example 13.12. Once we impose the constraint (13.65), the very simplest version of the dilation equation is (x) = 2 (2 x) (13.66) where c0 = 2 is the only (nonzero) coefficient. Up to constant multiple, the only "solutions" of the functional equation (13.66) with bounded support are scalar multiples of the delta function (x), which follows from the identity in Exercise 11.2.5. Other solutions, such as (x) = 1/x, are not localized, and thus not useful for constructing a wavelet basis. 2/25/07 715
c 2006 Peter J. Olver
The second condition we require is orthogonality of the wavelets. For simplicity, we only consider the standard L2 inner product f ,g =
-
f (x) g(x) dx.
It turns out that the orthogonality of the complete wavelet system is guaranteed once we know that the scaling function (x) is orthogonal to all its integer translates: (x) , (x - m) = We first note the formula (2 x - k) , (2 x - l) = =
-
(x) (x - m) dx = 0
for all
m = 0.
(13.67)
1 2
-
(2 x - k) (2 x - l) dx (x) (x + k - l) dx = 1 (x) , (x + k - l) 2
(13.68)
-
follows from the same change of variables y = 2 x - k used in (13.64). Therefore, since satisfies the dilation equation (13.57),
p p
(x) , (x - m) =
p
j =0
cj (2 x - j) ,
k=0
ck (2 x - 2 m - k)
p
(13.69)
=
j,k = 0
cj ck
1 (2 x - j) , (2 x - 2 m - k) = 2
j,k = 0
cj ck (x) , (x + j - 2 m - k) .
If we require orthogonality (13.67) of all the integer translates of , then the left hand side of this identity will be 0 unless m = 0, while only the summands with j = 2 m + k will be nonzero on the right. Therefore, orthogonality requires that c2m+k ck =
0 k p-2m
2, 0,
m = 0, m = 0.
(13.70)
The algebraic equations (13.65, 70) for the dilation coefficients are the key requirements for the construction of an orthogonal wavelet basis. For example, if we have just two nonzero coefficients c0 , c1 , then (13.65, 70) reduce to c0 + c1 = 2, c2 + c2 = 2, 0 1
and so c0 = c1 = 1 is the only solution, resulting in the Haar dilation equation (13.58). If we have three coefficients c0 , c1 , c2 , then (13.65), (13.70) require c0 + c1 + c2 = 2,
c2 + c2 + c2 = 2, 0 1 2
c0 c2 = 0.
In all instances, the functions have bounded support, and so the inner product integral can be reduced to an integral over a finite interval where both f and g are nonzero.
2/25/07
716
c 2006
Peter J. Olver
Thus either c2 = 0, c0 = c1 = 1, and we are back to the Haar case, or c0 = 0, c1 = c2 = 1, and the resulting dilation equation is a simple reformulation of the Haar case; see Exercise . In particular, the hat function (13.59) does not give rise to orthogonal wavelets. The remarkable fact, discovered by Daubechies, is that there is a nontrivial solution for four (and, indeed, any even number) of nonzero coefficients c0 , c1 , c2 , c3 . The basic equations (13.65), (13.70) require c0 + c1 + c2 + c3 = 2, The particular values c0 =
1+ 3 4 ,
c2 + c2 + c2 + c2 = 2, 0 1 2 3
3+ 3 4 , 3- 3 4 ,
c0 c2 + c1 c3 = 0.
1- 3 4 ,
(13.71)
c1 =
3+ 3 4
c2 =
3- 3 4
c3 =
1- 3 4
(13.72)
solve (13.71). These coefficients correspond to the Daubechies dilation equation (x) =
1+ 3 4
(2 x) +
(2 x - 1) +
(2 x - 2) +
(2 x - 3).
(13.73)
Any nonzero solution of bounded support to this remarkable functional equation will give rise to a scaling function (x), a mother wavelet w(x) =
1- 3 4
(2 x) -
3- 3 4
(2 x - 1) +
3+ 3 4
(2 x - 2) -
1+ 3 4
(2 x - 3),
(13.74)
and then, by compression and translation (13.62), the complete system of orthogonal wavelets wj,k (x). Before explaining how to solve the Daubechies dilation equation, let us complete the proof of orthogonality. It is easy to see that, by translation invariance, since (x) and (x - m) are orthogonal for any m = 0, so are (x - k) and (x - l) for any k = l. Next we prove orthogonality of (x - m) and w(x):
p p
w(x) , (x - m) = =
(-1)
j =0 p
j+1
cj (2 x - 1 + j) ,
k=0
ck (2 x - 2 m - k)
j,k = 0 p
(-1)j+1 cj ck (2 x - 1 + j) , (2 x - 2 m - k) (-1)j+1 cj ck (x) , (x - 1 + j - 2 m - k) ,
=
1 2
j,k = 0
using (13.68). By orthogonality (13.67) of the translates of , the only summands that are nonzero are when j = 2 m + k + 1; the resulting coefficient of (x) 2 is (-1)k c1-2 m-k ck = 0,
k
where the sum is over all 0 k p such that 0 1 - 2 m - k p. Each term in the sum appears twice, with opposite signs, and hence the result is always zero -- no matter what the coefficients c0 , . . . , cp are! The proof of orthogonality of the translates w(x - m) of the mother wavelet, along with all her wavelet descendants w(2j x - k), relies on a similar argument, and the details are left as an exercise for the reader. 2/25/07 717
c 2006 Peter J. Olver
Figure 13.10. Solving the Dilation Equation
Approximating the Daubechies Wavelet.
Let us next discuss how to solve the dilation equation (13.57). The solution we are after does not have an elementary formula, and we require a slightly sophisticated approach to recover it. The key observation is that (13.57) has the form of a fixed point equation = F [ ], not in ordinary Euclidean space, but in an infinite-dimensional function space. With luck, the fixed point (or, more correctly, fixed function) will be stable, and so starting with a suitable initial guess 0 (x), the successive iterates n+1 = F [ n ] will converge to the desired solution: n (x) - (x). In detail, the iterative version of the dilation equation (13.57) reads
p
n+1 (x) =
k=0
ck n (2 x - k),
n = 0, 1, 2, . . . .
(13.75)
Before attempting to prove convergence of this iterative procedure to the Daubechies scaling function, let us experimentally investigate what happens. A reasonable choice for the initial guess might be the Haar scaling or box function 0 (x) = 1, 0, 0 < t 1. otherwise.
In Figure 13.10 we graph the next 5 iterates 1 (x), . . . , 5 (x). There clearly appears to be converging to some function (x), although the final result does look a little bizarre. Bolstered by this preliminary experimental evidence, we can now try to prove convergence of the iterative scheme. This turns out to be true; a fully rigorous proof relies on the Fourier transform, [52], but is a little too advanced for this text and will be omitted. Theorem 13.13. The functions converge n (x) defined by the iterative functional equation (13.75) converge uniformly to a continuous function (x), called the Daubechies scaling function. 2/25/07 718
c 2006 Peter J. Olver
Once we have established convergence, we are now able to verify that the scaling function and consequential system of wavelets form an orthogonal system of functions. Proposition 13.14. All integer translates (x - k), for k Z of the Daubechies scaling function, and all wavelets wj,k (x) = w(2j x - k), j 0, are mutually orthogonal functions with respect to the L2 inner product. Moreover, 2 = 1, while wj,k 2 = 2- j . Proof : As noted earlier, the orthogonality of the entire wavelet system will follow once we know the orthogonality (13.67) of the scaling function and its integer translates. We use induction to prove that this holds for all the iterates n (x), and so, in view of uniform convergence, the limiting scaling function also satisfies this property. Details are relegated to Exercise . Q.E.D. In practical computations, the limiting procedure for constructing the scaling function is not so convenient, and an alternative means of computing its values is employed. The starting point is to determine its values at integer points. First, the initial box function has values 0 (m) = 0 for all integers m Z except 0 (1) = 1. The iterative functional equation (13.75) will then produce the values of the iterates n (m) at integer points m Z. A simple induction will convince you that n (m) = 0 except for m = 1 and m = 2, and, therefore, by (13.75), n+1 (1) =
3+ 3 4
n (1) +
1+ 3 4
n (2),
n+1 (2) =
1- 3 4
n (1) +
3- 3 4
n (2),
since all other terms are 0. This has the form of a linear iterative system v(n+1) = A v(n) with coefficient matrix A=
3+ 3 4 1- 3 4 1+ 3 4 3- 3 4
(13.76)
and where
v(n) =
n (1) . n (2)
Referring back to Chapter the 10, solution to such an iterative system is specified by the eigenvalues and eigenvectors of the coefficient matrix, which are 1 = 1, v1 =
1+ 3 4 1- 3 4
,
2 = 1 , 2
v2 =
-1 . 1
We write the initial condition as a linear combination of the eigenvectors 1- 3 1 0 (1) = 2 v1 - = v(0) = v2 . 0 0 (2) 2 The solution is v(n) = An v(0) 2/25/07 1 1- 3 1- 3 n A v2 = 2 v1 - n v2 . = 2 An v1 - 2 2 2 719
c 2006 Peter J. Olver
2 1.5 1 0.5 -0.5 -0.5 -1 -1.5 0.5 1 1.5 2 2.5 3 3.5
2 1.5 1 0.5 -0.5 -0.5 -1 -1.5 0.5 1 1.5 2 2.5 3 3.5
Figure 13.11. The limiting vector
The Daubechies Scaling Function and Mother Wavelet.
(1) (2)
= lim v(n) = 2 v1 =
n
1+ 3 2 1- 3 2
gives the desired values of the scaling function: 1- 3 = 1.366025 . . . , (1) = 2 (m) = 0, -1 + 3 (2) = = - .366025 . . . , 2 for all m = 1, 2.
(13.77)
With this in hand, the Daubechies dilation equation (13.73) then prescribes the func1 tion values 2 m at all half integers, because when x = 1 m then 2 x - k = m - k is an 2 integer. Once we know its values at the half integers, we can re-use equation (13.73) to give 1 its values at quarter integers 4 m. Continuing onwards, we determine the values of (x) at all dyadic points, meaning rational numbers of the form x = m/2j for m, j Z. Continuity will then prescribe its value at any other x R since x can be written as the limit of dyadic numbers xn -- namely those obtained by truncating its binary (base 2) expansion at the nth digit beyond the decimal (or, rather "binary") point. But, in practice, this latter step is unnecessary, since all computers are ultimately based on the binary number system, and so only dyadic numbers actually reside in a computer's memory. Thus, there is no real need to determine the value of at non-dyadic points. The preceding scheme was used to produce the graphs of the Daubechies scaling function in Figure 13.11. It is continuous, but non-differentiable function -- and its graph has a very jagged, fractal-like appearance when viewed at close range. The Daubechies scaling function is, in fact, a close relative of the famous example of a continuous, nowhere differentiable function originally due to Weierstrass, [127, 154], whose construction also relies on a similar scaling argument. With the values of the Daubechies scaling function on a sufficiently dense set of dyadic points in hand, the consequential values of the mother wavelet are given by formula (13.74). Note that supp = supp w = [ 0, 3 ]. The daughter wavelets are then found by the usual compression and translation procedure (13.62). 2/25/07 720
c 2006 Peter J. Olver
8 6 4 2
8 6 4 2
8 6 4 2
0.2 8 6 4 2
0.4
0.6
0.8
1 8 6 4 2
0.2
0.4
0.6
0.8
1 8 6 4 2
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
0.2
0.4
0.6
0.8
1
Figure 13.12.
Daubechies Wavelet Expansion.
The Daubechies wavelet expansion of a function whose support is contained in [ 0, 1 ] is then given by f (x) c0 (x) +
2j -1
cj,k wj,k (x).
(13.78)
j = 0 k = -2
The inner summation begins at k = -2 so as to include all the wavelet offspring wj,k whose support has a nontrivial intersection with the interval [ 0, 1 ]. The wavelet coefficients c0 , cj,k are computed by the usual orthogonality formula
3
c0 = f , =
f (x) (x) dx,
0 2-j (k+3) 3
(13.79) f 2- j (x + k) w(x) dx,
cj,k = f , wj,k = 2j
2-j
k
f (x) wj,k (x) dx =
0
where we agree that f (x) = 0 whenever x < 0 or x > 1. In practice, one employs a numerical integration procedure, e.g., the trapezoid rule, [9], based on dyadic nodes to speedily evaluate the integrals (13.79). A proof of completeness of the resulting wavelet basis functions can be found in [52]. Compression and denoising algorithms based on retaining only low frequency modes proceed as before, and are left as exercises for the reader to implement. Example 13.15. In Figure 13.12, we plot the Daubechies wavelet expansions of the same signal for Example 13.9. The first plot is the original signal, and the following show the partial sums of (13.78) over j = 0, . . . , r with r = 2, 3, 4, 5, 6. Unlike the Haar expansion, the Daubechies wavelets do exhibit the nonuniform Gibbs phenomenon at the interior discontinuity as well as the endpoints, since the function is set to 0 outside the interval [ 0, 1 ]. Indeed, the Daubechies wavelets are continuous, and so cannot converge uniformly to a discontinuous function.
For functions with larger support, one should include additional terms in the expansion corresponding to further translates of the wavelets so as to cover the entire support of the function. Alternatively, one can translate and rescale x to fit the function's support inside [ 0, 1 ].
2/25/07
721
c 2006
Peter J. Olver
13.3. The Fourier Transform.
Fourier series and their ilk are designed to solve boundary value problems on bounded intervals. The extension of Fourier methods to unbounded intervals -- the entire real line -- leads naturally to the Fourier transform, which is a powerful mathematical tool for the analysis of non-periodic functions. The Fourier transform is of fundamental importance in a broad range of applications, including both ordinary and partial differential equations, quantum mechanics, signal processing, control theory, and probability, to name but a few. We begin by motivating the Fourier transform as a limiting case of Fourier series. Although the rigorous details are rather exacting, the underlying idea is not so difficult. Let f (x) be a reasonably nice function defined for all - < x < . The goal is to construct a Fourier expansion for f (x) in terms of basic trigonometric functions. One evident approach is to construct its Fourier series on progressively larger and larger intervals, and then take the limit as the intervals' length becomes infinite. The limiting process converts the Fourier sums into integrals, and the resulting representation of a function is renamed the Fourier transform. Since we are dealing with an infinite interval, there are no longer any periodicity requirements on the function f (x). Moreover, the frequencies represented in the Fourier transform are no longer constrained by the length of the interval, and so we are effectively decomposing a quite general, non-periodic function into a continuous superposition of trigonometric functions of all possible frequencies. Let us present the details of this construction in a more concrete form. The computations will be significantly simpler if we work with the complex version of the Fourier series from the outset. Our starting point is the rescaled Fourier series (12.85) on a symmetric interval [ - , ] of length 2 , which we rewrite in the adapted form f (x)
= -
f (k ) i k x e . 2
(13.80)
The sum is over the discrete collection of frequencies , = 0, 1, 2, . . . , k = (13.81) corresponding to those trigonometric functions that have period 2 . For reasons that will soon become apparent, the Fourier coefficients of f are now denoted as c = f , e i k x = so that 1 2
f (x) e- i k x dx =
-
f (k ) , 2
1 f (x) e- i k x dx. (13.82) f (k ) = 2 - This reformulation of the basic Fourier series formula allows us to smoothly pass to the limit when the intervals' length . On an interval of length 2 , the frequencies (13.81) required to represent a function in Fourier series form are equally distributed, with interfrequency spacing k = k+1 - k = .
2/25/07
722
c 2006
Peter J. Olver
As , the spacing k 0, and so the relevant frequencies become more and more densely packed in the space of all possible frequencies: - < k < . In the limit, we anticipate that all possible frequencies will be represented. Indeed, letting k = k be arbitrary in (13.82), and sending , results in the infinite integral 1 f (k) = 2
f (x) e- i k x dx
(13.83)
-
known as the Fourier transform of the function f (x). If f (x) is a reasonably nice function, e.g., piecewise continuous and decayin to 0 reasonably quickly as | x | , its Fourier transform f (k) is defined for all possible frequencies - < k < . This formula will sometimes conveniently be abbreviated as f (k) = F[ f (x) ], (13.84)
where F is the Fourier transform operator . The extra 2 factor in front of the integral, which is sometimes omitted, is for later convenience. To reconstruct the function from its Fourier transform, we employ a similar limiting procedure on the Fourier series (13.80), which we first rewrite in a more suggestive form: 1 f (x) 2
= -
f (k ) e i k x k.
(13.85)
For each fixed value of x, the right hand side has the form of a Riemann sum, [9, 155], over the entire frequency space - < k < , for the function g (k) = f (k) e i k x . Under reasonable hypotheses, as , the functions f (k) f (k) converge to the Fourier transform; moreover, the interfrequency spacing k 0 and so one expects the Riemann sums to converge to the integral
1 f (x) f (k) e i k x dk, (13.86) 2 - known as the inverse Fourier transform, that serves to recover the original signal from its Fourier transform. In abbreviated form
is the inverse of the Fourier transform operator (13.84). In this manner, the Fourier series (13.85) becomes a Fourier integral that reconstructs the function f (x) as a (continuous) superposition of complex exponentials e i k x of all possible frequencies. The contribution of each such exponential is the value of the Fourier transform f (k) at the exponential's frequency. It is also worth pointing out that both the Fourier transform (13.84) and its inverse (13.87) define linear maps on function space. This means that the Fourier transform of the sum of two functions is the sum of their individual transforms, while multiplying a function by a constant multiplies its Fourier transform by the same factor: F[ f (x) + g(x) ] = F[ f (x) ] + F[ g(x) ] = f (k) + g(k), F[ c f (x) ] = c F[ f (x) ] = c f (k). 723 2/25/07
c 2006
f (x) = F -1 [ f (k) ]
(13.87)
(13.88)
Peter J. Olver
1.2 1
2
0.8 0.6 0.4 0.2 -2 -1 -0.2 1 2
-3 -2 -1 -0.5 1.2 1 0.8 0.6 0.4 0.2 1 -0.2 2 -2 -1 -0.2 1 2 1.5 1 0.5 1 2 3
1.2 1 0.8 0.6 0.4 0.2 -2 -1
Figure 13.13.
Fourier Transform of a Rectangular Pulse.
Recapitulating, by letting the length of the interval go to , the discrete Fourier series has become a continuous Fourier integral, while the Fourier coefficients, which were defined only at a discrete collection of possible frequencies, have become an entire function f (k) defined on all of frequency space k R. The reconstruction of f (x) from its Fourier transform f (k) via (13.86) can be rigorously justified under suitable hypotheses. For example, if f (x) is piecewise C1 on all of R and f (x) 0 decays reasonably rapidly as | x | ensuring that its Fourier integral (13.83) converges, then it can be proved that the inverse Fourier integral (13.86) will converge to f (x) at all points of continuity, and 1 to the midpoint 2 (f (x- ) + f (x+ )) at jump discontinuities -- just like a Fourier series. In particular, its Fourier transform f (k) 0 must also decay as | k | 0, implying that (as with Fourier series) the very high frequency modes make negligible contributions to the reconstruction of the signal. A more precise, general result will be formulated in Theorem 13.28 below. Example 13.16. The Fourier transform of a rectangular pulse or box function f (x) = (x + a) - (x - a) = of width 2 a is easily computed: 1 f (k) = 2
a
A similar statement hold for the inverse Fourier transform F -1 .
1, 0,
| x | > a.
- a < x < a,
(13.89)
e- i k x dx =
-a
e i ka - e- i ka = 2 i k
2 sin a k . k
(13.90)
On the other hand, the reconstruction of the pulse via the inverse transform (13.86) tells 2/25/07 724
c 2006 Peter J. Olver
us that 1
-
e
ikx
Note the convergence to the middle of the jump discontinuities at x = a. Splitting this complex integral into its real and imaginary parts, we deduce a pair of striking trigonometric integral identities 1, -a < x < a, cos k x sin a k sin k x sin a k 1 1 1 , x = a, dk = dk = 0. 2 - k - k 0, | x | > a, (13.92) Just as many Fourier series yield nontrivial summation formulae, the reconstruction of a function from its Fourier transform often leads to nontrivial integral identities. One cannot compute the integral (13.91) by the Fundamental Theorem of Calculus, since there is no elementary function whose derivative equals the integrand. Moreover, it is not even clear that the integral converges; indeed, the amplitude of the oscillatory integrand decays like 1/| k |, but the latter function does not have a convergent integral, and so the usual comparison test for infinite integrals, [9, 46, 167], fails to apply. Thus, the convergence of the integral is marginal at best: the trigonometric oscillations somehow overcome the slow rate of decay of 1/k and thereby induce the (conditional) convergence of the integral! In Figure 13.13 we display the box function with a = 1, its Fourier transform, along with a reconstruction obtained by numerically integrating (13.92). Since we are dealing with an infinite integral, we must break off the numerical integrator by restricting it to a finite interval. The first graph is obtained by integrating from -5 k 5 while the second is from -10 k 10. The non-uniform convergence of the integral leads to the appearance of a Gibbs phenomenon at the two discontinuities, similar to that of a Fourier series. Example 13.17. Consider an exponentially decaying right-handed pulse fr (x) = e- a x , 0, x > 0, x < 0, (13.93)
1, sin a k 1 , dk = f (x) = 2 k 0,
x = a, | x | > a.
-a < x < a,
(13.91)
where a > 0. We compute its Fourier transform directly from the definition: 1 f (k) = 2
0
e
-ax - i kx
e
1 e-(a+ i k)x dx = - 2 a + i k
x=0
=
1 . 2 (a + i k)
One can use Euler's formula (3.84) to reduce Zthe integrand to one of the form e k /k, but it can be proved, [ int ], that there is no formula for (e k /k) dk in terms of elementary functions.
Note that we can't Fourier transform the entire exponential function e- a x because it does not go to zero at both , which is required for the integral (13.83) to converge.
2/25/07
725
c 2006
Peter J. Olver
1 0.8 0.6 0.4 0.2 -3 -2 -1 -0.2 1 2 3 -3 -2 -1
1 0.8 0.6 0.4 0.2 1 -0.2 2 3
Right Pulse fr (x)
1 0.8 0.6 0.4 0.2 -3 -2
Left Pulse fl (x)
1 0.5
-1 -0.5
1
2
3
-3
-2
-1 -0.2
1
2
3 -1
Even Pulse fe (x) Figure 13.14.
Odd Pulse fo (x) Exponential Pulses.
As in the preceding example, the inverse Fourier transform for this function produces a nontrivial complex integral identity: e- a x , x > 0, ikx e 1 1 (13.94) dk = , x = 0, 2 2 - a + i k 0, x < 0. Similarly, a pulse that decays to the left, fl (x) = ea x , x < 0, 0, x > 0, where a > 0 is still positive, has Fourier transform f (k) = 1 . 2 (a - i k) (13.95)
(13.96)
This also follows from the general fact that the Fourier transform of f (- x) is f (- k); see Exercise . The even exponentially decaying pulse fe (x) = e- a | x | is merely the sum of left and right pulses: fe = fr + fl . Thus, by linearity, f e (k) = 1 1 + = 2 (a + i k) 2 (a - i k) a 2 , 2 + a2 k (13.98) (13.97)
The result is real and even because fe (x) is a real even function; see Exercise . The inverse Fourier transform (13.86) produces another nontrivial integral identity: e- a | x | = 2/25/07 1
-
a eikx a dk = 2 + a2 k 726
-
cos k x dk. k 2 + a2
c 2006
(13.99)
Peter J. Olver
(The imaginary part of the integral vanishes because the integrand is odd.) On the other hand, the odd exponentially decaying pulse, fo (x) = (sign x) e- a | x | = e- a x , - ea x , x > 0, x < 0, (13.100)
is the difference of the right and left pulses, fo = fr - fl , and has purely imaginary and odd Fourier transform f o (k) = 1 1 - = -i 2 (a + i k) 2 (a - i k) i
-
k 2 . 2 + a2 k
(13.101)
The inverse transform is (sign x) e- a | x | = - 1 k eikx dk = 2 + a2 k
-
k sin k x dk. k 2 + a2
(13.102)
As a final example, consider the rational function f (x) = x2 1 , + c2 where c > 0. (13.103)
Its Fourier transform requires integrating 1 f (k) = 2
-
e- i k x dx. x2 + a2
(13.104)
The indefinite integral (anti-derivative) does not appear in basic integration tables, and, in fact, cannot be done in terms of elementary functions. However, we have just managed to evaluate this particular integral! Look at (13.99). If we change x to k and k to - x, then we exactly recover the integral (13.104) up to a factor of a 2/. Therefore, in accordance with Theorem 13.18, we conclude that the Fourier transform of (13.103) is f (k) = e- a | k | . 2 a (13.105)
This last example is indicative of an important general fact. The reader has no doubt already noted the remarkable similarity between the Fourier transform (13.83) and its inverse (13.86). Indeed, the only difference is that the former has a minus sign in the exponential. This implies the following Symmetry Principle relating the direct and inverse Fourier transforms. Theorem 13.18. If the Fourier transform of the function f (x) is f (k), then the Fourier transform of f (x) is f (- k). Symmetry allows us to reduce the tabulation of Fourier transforms by half. For instance, referring back to Example 13.16, we deduce that the Fourier transform of the function 2 sin a x f (x) = x 2/25/07 727
c 2006 Peter J. Olver
is
Note that, by linearity, we can divide both f (x) and f (k) by 2/ to deduce the Fourier sin a x . transform of x Warning: Some authors leave out the 2 factor in the definition (13.83) of the Fourier transform f (k). This alternative convention does have a slight advantage of elim inating many 2 factors in the Fourier transforms of elementary functions. However, this requires an extra such factor in the reconstruction formula (13.86), which is achieved by replacing 2 by 2 . A significant disadvantage is that the resulting formulae for the Fourier transform and its inverse are not as close, and so the Symmetry Principle of Theorem 13.18 requires some modification. (On the other hand, convolution -- to be discussed below -- is a little easier without the extra factor.) When consulting any particular reference book, the reader always needs to check which convention is being used. See also Exercise . All of the functions in Example 13.17 required a > 0 for the Fourier integrals to converge. The functions that emerge in the limit as a goes to 0 are of fundamental importance. Let us start with the odd exponential pulse (13.100). When a 0, the function fo (x) converges to the sign function f (x) = sign x = (x) - (- x) = +1, -1, x > 0, x < 0. (13.107)
1, 1 , f (k) = (- k + a) - ( -k - a) = (k + a) - (k - a) = 2 0,
- a < k < a, | k | > a. k = a,
(13.106)
Taking the limit of the Fourier transform (13.101) leads to f (k) = - i 2 1 . k (13.108)
The nonintegrable singularity of f (k) at k = 0 is indicative of the fact that the sign function does not decay as | x | . In this case, neither the Fourier transform integral nor its inverse are well-defined as standard (Riemann, or even Lebesgue, [155]) integrals. Nevertheless, it is possible to rigorously justify these results within the framework of generalized functions. More interesting are the even pulse functions fe (x), which, in the limit a 0, become the constant function f (x) 1. (13.109) The limit of the Fourier transform (13.98) is lim 2 2a = 2 + a2 k 728 0, , k = 0, k = 0.
c 2006
a0
(13.110)
Peter J. Olver
2/25/07
This limiting behavior should remind the reader of our construction (11.31) of the delta function as the limit of the functions (x) = lim
n
n a = lim , 2 x2 ) 2 + x2 ) a 0 (a (1 + n
Comparing with (13.110), we conclude that the Fourier transform of the constant function (13.109) is a multiple of the delta function in the frequency variable: f (k) = 2 (k). (13.111) The direct transform integral (k) = 1 2
-
e- i k x dx
is, strictly speaking, not defined because the infinite integral of the oscillatory sine and cosine functions doesn't converge! However, this identity can be validly interpreted within the framework of weak convergence and generalized functions. On the other hand, the inverse transform formula (13.86) yields
-
(k) e i k x dk = e i k 0 = 1,
which is in accord with the basic definition (11.36) of the delta function. As in the previous case, the delta function singularity at k = 0 manifests the lack of decay of the constant function. Conversely, the delta function (x) has constant Fourier transform 1 (k) = 2 1 e- i k 0 . (x) e- i k x dx = 2 2 -
(13.112)
a result that also follows from the symmerty principle of Theorem 13.18. To determine the Fourier transform of a delta spike y (x) = (x - y) concentrated at position x = y, we compute 1 e- i k y y (k) = . (13.113) (x - y) e- i k x dx = 2 - 2 The result is a pure exponential in frequency space. Applying the inverse Fourier transform (13.86) leads, formally, to the remarkable identity 1 y (x) = (x - y) = 2
-
e- i k(x-y) dk =
1 eiky , eikx , 2
(13.114)
where , denotes the usual L2 inner product on R. Since the delta function vanishes for x = y, this identity is telling us that complex exponentials of differing frequencies are mutually orthogonal. However, this statement must be taken with a grain of salt, since the integral does not converge in the normal (Riemann or Lebesgue) sense. But 2/25/07 729
c 2006 Peter J. Olver
it is possible to make sense of this identity within the language of generalized functions. Indeed, multiplying both sides by f (x), and then integrating with respect to x, we find f (y) = 1 2
- -
f (x) e- i k(x-y) dx dk.
(13.115)
This is a perfectly valid formula, being a restatement (or, rather, combination) of the basic formulae (13.83, 86) connecting the direct and inverse Fourier transforms of the function f (x). Conversely, the Symmetry Principle tells that the Fourier transform of a pure us i lx exponential e will be a shifted delta spike 2 (k - l), concentrated in frequency space. Both results are particular cases of the general Shift Theorem, whose proof is left as an exercise for the reader. Theorem 13.19. If f (x) has Fourier transform f (k), then the Fourier transform of the shifted function f (x-y) is e- i k y f (k). Similarly, the transform of the product function e i l x f (x) for l real is the shifted transform f (k - l). Since the Fourier transform uniquely associates a function f (k) on frequency space with each (reasonable) function f (x) on physical space, one can characterize functions by their transforms. Practical applications rely on tables (or, even better, computer algebra systems) that recognize a wide variety of transforms of basic functions of importance in applications. The accompanying table lists some of the most important examples of functions and their Fourier transforms. Note that, according to the Symmetry Principle of Theorem 13.18, each tabular entry can be used to deduce two different Fourier transforms. A more extensive collection of Fourier transforms can be found in [140]. Derivatives and Integrals One of the most remarkable and important properties of the Fourier transform is that it converts calculus into algebra! More specifically, the two basic operations in calculus -- differentiation and integration of functions -- are realized as algebraic operations on their Fourier transforms. (The downside is that algebraic operations become more complicated in the transform domain.) Let us begin with derivatives. If we differentiate the basic inverse Fourier transform formula 1 f (x) f (k) e i k x dk. 2 - with respect to x, we obtain
1 f (x) i k f (k) e i k x dk. (13.116) 2 - The resulting integral is itself in the form of an inverse Fourier transform, namely of 2 i k f (k) which immediately implies the following basic result.
We are assuming the integrand is sufficiently nice so that we can bring the derivative under the integral sign; see [ 64, 195 ] for a fully rigorous derivation.
2/25/07
730
c 2006
Peter J. Olver
Short Table of Fourier Transforms
f (x) 1 (x) (x) sign x (x + a) - (x - a) e- a x (x) ea x (1 - (x)) e- a | x | e
- a x2
f (k) 2 (k) 1 2
i (k) - 2 2 k -i 2 1 k
2 sin a k k 1 2 (a + i k) 1 2 (a - i k) 2 a 2 + a2 k e- k /4 a 2a -i e- | k | 3/2 + (k) 2 k 2 e i k d/c f |c| f (- k) f (- k) i k f (k) i f (k) 2 f (k) g(k) k c
2
tan-1 x f (c x + d) f (x) f (x) f (x) x f (x) f g(x)
Note: The parameter a > 0 is always real and positive, while c = 0 is any nonzero complex number. 2/25/07 731
c 2006 Peter J. Olver
Proposition 13.20. The Fourier transform of the derivative f (x) of a function is obtained by multiplication of its Fourier transform by i k: F[ f (x) ] = i k f (k). (13.117)
Similarly, the Fourier transform of x f (x) is obtained by differentiating the Fourier transform of f (x): df F[ x f (x) ] = i k . (13.118) dk The second statement follows from the first by use of the Symmetry Principle of Theorem 13.18. Example 13.21. The derivative of the even exponential pulse fe (x) = e- a | x | is a multiple of the odd exponential pulse fo (x) = (sign x) e- a | x | : fe (x) = - a (sign x) e- a | x | = - afo (x). Proposition 13.20 says that their Fourier transforms are related by f o (k) = - ik f (k) = - i a e ik k 2 =- f (k), 2 + a2 k a o
as previously noted in (13.98, 101). On the other hand, since the odd exponential pulse has a jump discontinuity of magnitude 2 at x = 0, its derivative contains a delta function, and is equal to fo (x) = - a e- a | x | + 2 (x) = - a fe (x) + 2 (x). This is reflected in the relation between their Fourier transforms. If we multiply (13.101) by i k we obtain i k f o (k) = 2 k2 = k 2 + a2 2 - 2 a2 = 2 (k) - a f e (k). k 2 + a2
The Fourier transform, just like Fourier series, is completely compatible with the calculus of generalized functions. Higher order derivatives are handled by iterating formula (13.117), and so: Corollary 13.22. The Fourier transform of f (n) (x) is ( i k)n f (k). This result has an important consequence: the smoothness of f (x) is manifested in the rate of decay of its Fourier transform f (k). We already noted that the Fourier transform of a (nice) function must decay to zero at large frequencies: f (k) 0 as | k | . If the nth derivative f (n) (x) is also a reasonable function, then its Fourier transform f (n) (k) = ( i k)n f (k) must go to zero as | k | . This requires that f (k) go to zero more rapidly than | k |- n . Thus, the smoother f (x), the more rapid the decay of its Fourier transform. As a general rule of thumb, local features of f (x), such as smoothness, are manifested by global features of f (k), such as decay for large | k |. The Symmetry 2/25/07 732
c 2006 Peter J. Olver
Principle implies that reverse is also true: global features of f (x) correspond to local features of f (k). This local-global duality is one of the major themes of Fourier theory. Integration of functions is the inverse operation to differentiation, and so should correspond to division by 2 i k in frequency space. As with Fourier series, this is not completely correct; there is an extra constant involved, and this contributes an extra delta function in frequency space. Proposition 13.23. If f (x) has Fourier transform f (k), then the Fourier transform
x
of its integral g(x) =
-
f (y) dy is g(k) = - i f (k) + f (0) (k). k
-
(13.119)
Proof : First notice that
x -
lim
g(x) = 0,
x +
lim
g(x) =
f (x) dx =
2 f (0).
Therefore, by subtracting a suitable multiple of the step function from the integral, the resulting function h(x) = g(x) - 2 f (0) (x) decays to 0 at both . Consulting our table of Fourier transforms, we find h(k) = g(k) - f (0) (k) + On the other hand, h (x) = f (x) - i f (0) . k (13.120)
2 f (0) (x).
Since h(x) 0 as | x | , we can apply our differentiation rule (13.117), and conclude that i k h(k) = f (k) - f (0). (13.121) Combining (13.120) and (13.121) establishes the desired formula (13.119). Example 13.24. The Fourier transform of the inverse tangent function
x
Q.E.D.
f (x) = tan-1 x =
0
dy = 1 + y2
x -
dy - 2 1+y 2
can be computed by combining Proposition 13.23 with (13.105): f (k) = - i Applications to Differential Equations The fact that the Fourier transform changes differentiation in the physical domain into multiplication in the frequency domain is one of its most compelling features. A 2/25/07 733
c 2006 Peter J. Olver
e- | k | 3/2 + (k). 2 k 2
particularly important consequence is that the Fourier transform effectively converts differential equations into algebraic equations, and thereby opens the door to their solution by elementary algebra! One begins by applying the Fourier transform to both sides of the differential equation under consideration. Solving the resulting algebraic equation will produce a formula for the Fourier transform of the desired solution, which can then be immediately reconstructed via the inverse Fourier transform. The Fourier transform is particularly well adapted to boundary value problems on the entire real line. In place of the boundary conditions used on finite intervals, we look for solutions that decay to zero sufficiently rapidly as | x | -- in order that their Fourier transform be well-defined (in the context of ordinary functions). In quantum mechanics, [129], these solutions are known as the bound states of the system, and correspond to subatomic particles that are trapped or localized in a region of space by some sort of force field. For example, the bound electrons in an atom are localized by the electrostatic attraction of the nucleus. As a specific example, consider the boundary value problem - d2 u + 2 u = h(x), dx2 - < x < , (13.122)
where > 0 is a positive constant. In lieu of boundary conditions, we require that the solution u(x) 0 as | x | . We will solve this problem by applying the Fourier transform to both sides of the differential equation. Taking Corollary 13.22 into account, the result is a linear algebraic equation k 2 u(k) + 2 u(k) = h(k) relating the Fourier transforms of u and h. Unlike the differential equation, the transformed equation can be immediately solved for u(k) = h(k) + 2 (13.123)
k2
Therefore, we can reconstruct the solution by applying the inverse Fourier transform formula (13.86): 1 u(x) = 2
-
h(k) e i k x dk. k2 + 2
(13.124)
For example, if the forcing function is an even exponential pulse, h(x) = e- | x | with h(k) = 1 2 , k2 + 1
then (13.124) writes the solution as a Fourier integral: u(x) = 2/25/07 1
-
eikx 1 dk = 2 + 2 )(k 2 + 1) (k 734
-
(k 2
cos k x dk . + 2 )(k 2 + 1)
c 2006 Peter J. Olver
noting that the imaginary part of the complex integral vanishes because the integrand is an odd function. (Indeed, if the forcing function is real, the solution must also be real.) The Fourier integral can be explicitly evaluated by using partial fractions to rewrite u(k) = 1 2 = 2 + 2 )(k 2 + 1) (k 1 2 2-1 k2 1 1 - 2 + 1 k + 2 , 2 = 1,
Thus, according to our Fourier Transform Table, the solution to this boundary value problem is 1 e- | x | - e- | x | u(x) = when 2 = 1. (13.125) 2 - 1 The reader may wish to verify that this function is indeed a solution, meaning that it is twice continuously differentiable (which is not so immediately apparent from the formula), decays to 0 as | x | , and satisfies the differential equation everywhere. The "resonant" case 2 = 1 is left as an exercise. Remark : The method of partial fractions that you learned in first year calculus is often an effective tool for evaluating (inverse) Fourier transforms of rational functions. A particularly important case is when the forcing function h(x) = y (x) = (x - y) represents a unit impulse concentrated at x = y. The resulting square-integrable solution is the Green's function G(x, y) for the boundary value problem. According to (13.123), its Fourier transform with respect to x is e- i k y 1 G(k, y) = . 2 k2 + 2 which is the product of an exponential factor e- i k y , representing the Fourier transform of y (x), times a multiple of the Fourier transform of the even exponential pulse e- | x | . We apply Theorem 13.19, and conclude that the Green's function for this boundary value problem is an exponential pulse centered at y, namely 1 - | x-y | e . 2 Observe that, as with other self-adjoint boundary value problems, the Green's function is symmetric under interchange of x and y. As a function of x, it satisfies the homogeneous differential equation - u + 2 u = 0, except at the point x = y when its derivative has a jump discontinuity of unit magnitude. It also decays as | x | , as required by the boundary conditions. The Green's function superposition principle tells us that the solution to the inhomogeneous boundary value problem (13.122) under a general forcing can be represented in the integral form G(x, y) = 1 u(x) = G(x, y) h(y) dy = 2 -
-
e- | x-y | h(y) dy.
(13.126)
The reader may enjoy recovering the particular exponential solution (13.125) from this integral formula. 2/25/07 735
c 2006 Peter J. Olver
Convolution The final Green's function formula (13.126) is indicative of a general property of Fourier transforms. The right hand side has the form of a convolution product between functions. Definition 13.25. The convolution of scalar functions f (x) and g(x) is the scalar function h = f g defined by the formula h(x) = f (x) g(x) =
-
f (x - y) g(y) dy.
(13.127)
We record the basic properties of the convolution product, leaving their verification as exercises for the reader. All of these assume that the implied convolution integrals converge. (a) Symmetry: f g = g f, f (a g + b h) = a (f g) + b (f h), (b) Bilinearity: a, b C, (a f + b g) h = a (f h) + b (g h), (c) Associativity: f (g h) = (f g) h, (d) Zero function: f 0 = 0, (e) Delta function: f = f. One tricky feature is that the constant function 1 is not a unit for the convolution product; indeed, f 1=
f (y) dy
-
is a constant function -- the total integral of f -- not the original function f (x). In fact, according to the last property, the delta function plays the role of the "convolution unit": f (x) (x) =
-
f (x - y) (y) dy = f (x),
which follows from the basic property (11.36) of the delta function. In particular, our solution (13.126) has the form of a convolution product between an 1 - | x | e and the forcing function: even exponential pulse g(x) = 2 u(x) = g(x) h(x). On the other hand, its Fourier transform (13.123) is the ordinary multiplicative product u(k) = g(k) h(k) of the Fourier transforms of g and h. In fact, this is a general property of the Fourier transform: convolution in the physical domain corresponds to multiplication in the frequency domain, and conversely. 2/25/07 736
c 2006 Peter J. Olver
Theorem 13.26. The Fourier transform of the convolution u(x) = f (x) g(x) of two functions is a multiple of the product of their Fourier transforms: u(k) = 2 f (k) g(k). Vice versa, the Fourier transform of their product h(x) = f (x) g(x) is, up to multiple, the convolution of their Fourier transforms: 1 h(k) = f (k) g(k). 2 Proof : Combining the definition of the Fourier transform with the convolution formula (13.127), we find 1 u(k) = 2 1 h(x) e- i k x dx = 2 -
- -
f (x - y) g(y) e- i k x dx dy.
where we are assuming that the integrands are sufficiently nice to allow us to interchange the order of integration, [9]. Applying the change of variables z = x - y in the inner integral produces 1 u(k) = 2 = 2
-
f (z) g(y) e- i k(y+z) dy dz f (z) e
-ikz
1 2
- -
dz
1 2
-
g(y) e- i k y dy
=
2 f (k) g(k).
The second statement follows directly from the Symmetry Principle of Theorem 13.18. Q.E.D. Example 13.27. We already know, (13.106), that the Fourier transform of f (x) = is the box function f (k) = sin x x , 2 0,
We also know that the Fourier transform of g(x) = 1 x is
(k + 1) - (k - 1) = 2
-1 < k < 1, | k | > 1, sign k. 2
g(k) = - i
Therefore, the Fourier transform of their product h(x) = f (x) g(x) = 2/25/07 737 sin x x2
c 2006 Peter J. Olver
1 0.5 -3 -2 -1 -0.5 -1 1 2 3
Figure 13.15.
The Fourier transform of
sin x . x2
can be obtained by convolution: 1 1 f g(k) = h(k) = 2 2
-
= -i
8
A graph of h(k) appears in Figure 13.15.
i 2 1 sign(k - l) dl = k, -i 2 -1 -i 2
f (l) g(k - l) dl
k < -1, -1 < k < 1, k > 1.
The Fourier Transform on Hilbert Space While we do not have the space to embark on a fully rigorous treatment of the theory underlying the Fourier transform, it is worth outlining a few of the most basic features. We have already noted that the Fourier transform, when defined, is a linear map, taking functions f (x) on physical space to functions f (k) on frequency space. A critical question is precisely which function space should the theory be applied to. Not every function admits a Fourier transform in the classical sense -- the Fourier integral (13.83) is required to converge, and this places restrictions on the function and its asymptotics at large distances. It turns out the proper setting for the rigorous theory is the Hilbert space of complexvalued square-integrable functions -- the same infinite-dimensional vector space that lies at the heart of modern quantum mechanics. In Section 12.5, we already introduced the Hilbert space L2[ - , ] on a finite interval; here we adapt Definition 12.33 to the entire real line. Thus, the Hilbert space L2 = L2 (R) is the infinite-dimensional vector space
We leave aside the more advanced issues involving generalized functions in this subsection.
2/25/07
738
c 2006
Peter J. Olver
consisting of all complex-valued functions f (x) which are defined for all x R and have finite L2 norm: f
2
=
For example, any piecewise continuous function that satisfies the decay criterion | f (x) | M , | x |1/2+ for all sufficiently large | x | 0, (13.129)
-
| f (x) |2 dx < .
(13.128)
for some M > 0 and > 0, belongs to L2 . However, as in Section 12.5, Hilbert space contains many more functions, and the precise definitions and identification of its elements is quite subtle. The inner product on the Hilbert space L2 is prescribed in the usual manner, f ,g = The CauchySchwarz inequality
f (x) g(x) dx.
-
ensures that the inner product integral is finite whenever f, g L2 . Let us state the fundamental theorem governing the effect of the Fourier transform on functions in Hilbert space. It can be regarded as a direct analog of the pointwise convergence Theorem 12.7 for Fourier series. A fully rigorous proof can be found in [175, 195]. Theorem 13.28. If f (x) L2 is square-integrable, then its Fourier transform f (k) L is a well-defined, square-integrable function of the frequency variable k. If f (x) is continuously differentiable at a point x, then its inverse Fourier transform (13.86) equals its value f (x). More generally, if the right and left hand limits f (x- ), f (x- ), and f (x+ ), f (x+ ) exist, then the inverse Fourier transform integral converges to the average value 1 - + 2 f (x ) + f (x ) .
2
| f ,g | f
g
Thus, the Fourier transform f = F[ f ] defines a linear transformation from L2 functions of x to L2 functions of k. In fact, the Fourier transform preserves inner products. This important result is known as Parseval's formula, whose Fourier series counterpart appeared in (12.117). Theorem 13.29. If f (k) = F[ f (x) ] and g(k) = F[ g(x) ], then f , g = f , g , i.e.,
f (x) g(x) dx =
f (k) g(k) dk.
(13.130)
-
-
Proof : Let us sketch a formal proof that serves to motivate why this result is valid. We use the definition (13.83) of the Fourier transform to evaluate
-
f (k) g(k) dk = =
- -
1 2
-
-
f (x) e- i k x dx 1 2
-
1 2
-
g(y) e+ i k y dy
dk
f (x) g(y) 739
e- i k(x-y) dk
dx dy.
c 2006 Peter J. Olver
2/25/07
Now according to (13.114), the inner k integral can be replaced by a delta function (x-y), and hence
-
f (k) g(k) dk =
-
-
f (x) g(y) (x - y) dx dy =
-
f (x) g(x) dx. Q.E.D.
This completes our "proof"; see [175, 195] for a rigorous version.
In particular, orthogonal functions, f , g = 0, will have orthogonal Fourier transforms, f , g = 0. Choosing f = g in (13.130) results in the Plancherel formula f 2 = f 2 , or, explicitly,
-
| f (x |2 dx =
-
| f (k) |2 dk.
(13.131)
We conclude that the Fourier transform map F: L2 L2 defines a unitary or normpreserving linear transformation on Hilbert space, mapping L2 functions of the physical variable x to L2 functions of the frequency variable k. Quantum Mechanics and the Uncertainty Principle In quantum mechanics, the wave functions or probability densities of a quantum system are characterized as unit elements, = 1 of the underlying state space, which, in a one-dimensional model of a single particle, is the Hilbert space L2 = L2 (R). All measurable physical quantities are represented as linear operators A: L2 L2 acting on (a suitable dense subspace of) the state space. These are obtained by the tahter myusterious process of "quantizing" their classical counterparts, which are ordinary functions. In particular, the particle's position, usually denoted by Q, is represented by the operation of multiplication by x, so Q[ ] = x (x), whereas its momentum, denoted by P for historical reasons, is represented by the differentiation operator P = d/dx, so that P [ ] = (x). In its popularized form, Heisenberg's Uncertainty Principle is a by now familiar philosophical concept. The principle, first formulated by the twentieth century German physicist Werner Heisenberg, one of the founders of modern quantum mechanics, states that, in a physical system, certain quantities cannot be simultaneously measured with complete accuracy. For instance, the more precisely one measures the position of a particle, the less accuracy there will be in the measurement of its momentum; vice versa, the greater the accuracy in the momentum, the less certainty in its position. A similar uncertainty couples energy and time. Experimental verification of the uncertainty principle can be found even in fairly simple situations. Consider a light beam passing through a small hole. The position of the photons is constrained by the hole; the effect of their momenta is in the pattern of light diffused on a screen placed beyond the hole. The smaller the hole, the more constrained the position, and the wider the image on the screen, meaning the less certainty there is in the momentum. This is not the place to discuss the philosophical and experimental consequences of Heisenberg's principle. What we will show is that the Uncertainty Principle is, in fact, a rigorous theorem concerning the Fourier transform! In quantum theory, each of the paired quantities, e.g., position and momentum, are interrelated by the Fourier transform. Indeed, 2/25/07 740
c 2006 Peter J. Olver
Proposition 13.20 says that the Fourier transform of the differentiation operator representing monmentum is a multiplication operator representing position on the frequency space and vice versa. This Fourier transform-based duality between position and momentum, or multiplication and differentiation, lies at the heart of the Uncertainty Principle. In general, suppose A is a linear operator on Hilbert space representing a physical quantity. If the quantum system is in a state represented by a particular wave function , then the localization of the quantity A is measured by the norm A[ ] . The smaller this norm, the more accurate the measurement. Thus, Q[ ] = x (x) measures the localization of the position of the particle represented by ; the smaller Q[ ] , the more concentrated the probability of finding the particle near x = 0 and hence the smaller the error in the measurement of its position. Similarly, by Plancherel's formula (13.131), P [ ] = (x) = i k (k) = k (k) ,
measures the localization in the momentum of the particle, which is small if and only if its Fourier transform is concentrated near k = 0, and hence the smaller the error in its measured momentum. With this interpretation, the Uncertainty Principle states that these two quantities cannot simultaneously be arbitrarily small. Theorem 13.30. If (x) is a wave function, so = 1, then Q[ ] P[] 1. 2 (x) = Q[ ] P[] . (13.132)
Proof : The proof rests on the CauchySchwarz inequality x (x) , (x) x (x) (13.133)
On the other hand, writing out the inner product term x (x) , (x) = Let us integrate by parts, using the fact that d 1 2 2 (x) . dx Since (x) 0 as | x | , the boundary terms vanish, and hence (x) (x) = x (x) , (x) = since
- -
x (x) (x) dx.
x (x) (x) dx = -
-
1 2
(x)2 dx = - 1 , 2 Q.E.D.
= 1. Substituting back into (13.133) completes the proof.
The inequality (13.132) quantifies the statement that the more accurately we measure the momentum Q, the less accurately we are able to measure the position P , and vice versa. For more details and physical consequences, you should consult an introductory text on mathematical quantum mechanics, e.g., [124, 129].
At another position a, one replaces x by x - a.
2/25/07
741
c 2006
Peter J. Olver
13.4. The Laplace Transform.
In engineering applications, the Fourier transform is often overshadowed by a close relative. The Laplace transform plays an essential role in control theory, linear systems analysis, electronics, and many other fields of practical engineering and science. However, the Laplace transform is most properly interpreted as a particular real form of the more fundamental Fourier transform. When the Fourier transform is evaluated along the imaginary axis, the complex exponential factor becomes real, and the result is the Laplace transform, which maps real-valued functions to real-valued functions. Since it is so closely allied to the Fourier transform, the Laplace transform enjoys many of its featured properties, including linearity. Moreover, derivatives are transformed into algebraic operations, which underlies its applications to solving differential equations. The Laplace transform is one-sided; it only looks forward in time and prefers functions that decay -- transients. The Fourier transform looks in both directions and prefers oscillatory functions. For this reason, while the Fourier transform is used to solve boundary value problems on the real line, the Laplace transform is much better adapted to initial value problems. Since we will be applying the Laplace transform to initial value problems, we switch our variable from x to t to emphasize this fact. Suppose f (t) is a (reasonable) function which vanishes on the negative axis, so f (t) = 0 for all t < 0. The Fourier transform of f is 1 f (k) = f (t) e- i k t dt, 2 0 since, by our assumption, its negative t values make no contribution to the integral. The Laplace transform of such a function is obtained by replacing i k by a real variable s, leading to where, in accordance with the standard convention, the factor of 2 has been omitted. By allowing complex values of the Fourier frequency variable k, we may identify the Laplace transform with 2 times the evaluation of the Fourier transform for values of k = - i s on the imaginary axis: F (s) = 2 f (- i s). (13.135) Since the exponential factor in the integral has become real, the Laplace transform L takes real functions to real functions. Moreover, since the integral kernel e- s t is exponentially decaying for s > 0, we are no longer required to restrict our attention to functions that decay to zero as t . Example 13.31. Consider an exponential function f (t) = e t , where the exponent is allowed to be complex. Its Laplace transform is F (s) =
0
F (s) = L[ f (t) ] =
f (t) e- s t dt,
(13.134)
0
e(-s) t dt =
1 . s-
(13.136)
One can also define the Laplace transform at complex values of s, but this will not be required in the applications discussed here.
2/25/07
742
c 2006
Peter J. Olver
Note that the integrand is exponentially decaying, and hence the integral converges, if and only if Re ( - s) < 0. Therefore, the Laplace transform (13.136) is, strictly speaking, only defined at sufficiently large s > Re . In particular, for an oscillatory exponential, L[ e i t ] = 1 s+ i = 2 s- i s + 2 provided s > 0.
Taking real and imaginary parts of this identity, we discover the formulae for the Laplace transforms of the simplest trigonometric functions: L[ cos t ] =
0
s , s2 + 2
0
L[ sin t ] = 1 s
. s2 + 2
(13.137)
Two additional important transforms are L[ 1 ] = e- s t dt = 1 , s L[ t ] = t e- s t dt = e- s t dt = 1 . s2 (13.138)
0
The second computation relies on an integration by parts, making sure that the boundary terms at s = 0, vanish. Remark : In every case, we really mean the Laplace transform of the function whose values are given for t > 0 and is equal to 0 for all negative t. Therefore, the function 1 in reality signifies the step function (t) = 1, 0, t > 0, t < 0, (13.139)
and so the first formula in (13.138) should more properly be written 1 . (13.140) s However, in the traditional approach to the Laplace transform, one only considers the functions on the positive t axis, and so the step function and the constant function are, from this viewpoint, indistinguishable. However, once one moves beyond a purely mechanistic approach, any deeper understanding of the properties of the Laplace transform requires keeping this distinction firmly in mind. L[ (t) ] = Let us now pin down the precise class of functions to which the Laplace transform can be applied. Definition 13.32. A function f (t) is said to have exponential growth of order a if | f (t) | < M ea t , for some M > 0 and t0 > 0. Note that the exponential growth condition only depends upon the function's behavior for large values of t. If a < 0, then f is, in fact, exponentially decaying as x . Since ea t < eb t for a < b and all t > 0, if f (t) has exponential growth of order a, it automatically has exponential growth of any higher order b > a. All polynomial, trigonometric, 2/25/07 743
c 2006 Peter J. Olver
for all
t > t0 ,
(13.141)
and exponential functions (with linear argument) have exponential growth. The simplest 2 example of a function that does not satisfy any exponential growth bound is f (t) = et , since it grows faster than any simple exponential ea t . The following result guarantees the existence of the Laplace transform, at least for sufficiently large values of the transform variable s, for a rather broad class of functions that includes almost all of the functions that arise in applications. Theorem 13.33. If f (t) is piecewise continuous and has exponential growth of order a, then its Laplace transform F (s) = L[ f (t) ] is defined for all s > a. Proof : The exponential growth inequality (13.141) implies that we can bound the integrand in (13.134) by | f (t) e- s t | < M e(a-s) t . Therefore, as soon as s > a, the integrand is exponentially decaying as t , and this suffices to ensure the convergence of the Laplace transform integral. Q.E.D. Theorem 13.33 is an existential result, and of course, in practice, we may not be able to explicitly evaluate the Laplace transform integral. Nevertheless, the Laplace transforms of most common functions are not hard to find, and extensive lists have been tabulated, [141]. An abbreviated table of Laplace transforms can be found on the following page. Nowadays, the most convenient sources of transform formulas are computer algebra packages, including Mathematica and Maple. According to Theorem 13.28, when it exists, the Fourier transform uniquely specifies the function, except possibly at jump discontinuities where the limiting value must be half way in between. An analogous result can be established for the Laplace transform. Lemma 13.34. If f (t) and g(t) are piecewise continuous functions that are of exponential growth, and L[ f (t) ] = L[ g(t) ] for all s sufficiently large, then f (t) = g(t) at all points of continuity t > 0. In fact, there is an explicit formula for the inverse Laplace transform, which follows from its identification, (13.135), with the Fourier transform along the imaginary axis. Under suitable hypotheses, a given function F (s) is the Laplace transform of the function f (t) determined by the complex integral formula f (t) = 1 2 i
i -i
F (s) es t ds,
t > 0.
(13.142)
In practice, one hardly ever uses this complicated formula to compute the inverse Laplace transform. Rather, one simply relies on tables of known Laplace transforms, coupled with a few basic rules that will be covered in the following subsection.
See Section 16.5 for details on complex integration. The stated formula doesn't apply to all functions of exponential growth. A more universally valid inverse Laplace transform formula is obtained by shifting the complex contour to run from b - i to b + i for some b > a, the order of exponential growth of f .
2/25/07
744
c 2006
Peter J. Olver
Table of Laplace Transforms
f (t) 1 t tn (t - c) e t cos t sin t ec t f (t) (t - c) f (t - c) t f (t) f (t) f (n) (t) 1 s 1 s2 n! sn+1 e- s c 1 s- s2 s2
F (s)
s + 2 + 2
F (s - c) e- s c F (s) - F (s) s F (s) - f (0) sn F (s) - sn-1 f (0) - - sn-2 f (0) - - f (n-1) (0) F (s) G(s)
f (t) g(t)
In this table, n is a non-negative integer, is any real number, while c 0 is any non-negative real number.
2/25/07
745
c 2006
Peter J. Olver
The Laplace Transform Calculus The first observation is that the Laplace transform is a linear operator, and so L[ f + g ] = L[ f ] + L[ g ], L[ c f ] = c L[ f ], (13.143)
for any constant c. Moreover, just like its Fourier progenitor, the Laplace transform converts calculus into algebra. In particular, differentiation turns into multiplication by the transform variable, but with one additional term that depends upon the value of the function at t = 0. Theorem 13.35. Let f (t) be continuously differentiable for t > 0 and have exponential growth of order a. If L[ f (t) ] = F (s) then, for s > a, L[ f (t) ] = s F (s) - f (0). Proof : The proof relies on an integration by parts: L[ f (t) ] =
0 t -st
(13.144)
f (t) e
-st
dt = f (t) e
-st
t=0
+s
0
f (t) e- s t dt
= lim f (t) e
- f (0) + s F (s).
The exponential growth inequality (13.141) implies that first term vanishes for s > a, and the remaining terms agree with (13.144). Q.E.D. Example 13.36. According to Example 13.31, the Laplace transform of the function sin t is L[ sin t ] = 2 . s + 2 Its derivative is d sin t = cos t, dt and therefore L[ cos t ] = s L[ sin t ] = s2 s , + 2
since sin t vanishes at t = 0. The result agrees with (13.137). On the other hand, d cos t = - sin t, dt and so L[ - sin t ] = s L[ cos t ] - 1 = again in agreement with the known formula. 2/25/07 746
c 2006 Peter J. Olver
s2 2 -1=- 2 , s2 + 2 s + 2
Remark : The final term - f (0) in (13.144) is a manifestation of the discontinuity in f (t) at t = 0. Keep in mind that the Laplace transform only applies to functions with f (t) = 0 for all t < 0, and so if f (0) = 0, then f (t) has a jump discontinuity of magnitude f (0) at t = 0. Therefore, by the calculus of generalized functions, its derivative f (t) should include a delta function term, namely f (0) (0), which accounts for the additional constant term in its transform. In the practical approach to the Laplace transform calculus, one suppresses the delta function when computing the derivative f (t). However, its effect must reappear on the other side of the differentiation formula (13.144), and the upshot is the extra term - f (0). Laplace transforms of higher order derivatives are found by iterating the first order formula (13.144). For example, if f C2 , then In general, for an n times continuously differentiable function, L[ f (t) ] = s L[ f (t) ] - f (0) = s2 F (s) - s f (0) - f (0). (13.145) (13.146)
On the other hand, multiplying the function by t corresponds to differentiation of its Laplace transform (up to a change of sign): The proof of this formula is left as an exercise for the reader. Conversely, integration corresponds to dividing the Laplace transform by s, so F (s) . (13.148) s 0 Unlike the Fourier transform, there are no additional terms in the integration formula as long as we start the integral at t = 0. For instance, L f ( ) d = 1 2 n! L[ 2 t ] = 3 , and, more generally, L[ tn ] = n+1 . (13.149) s s s There is also a shift formula, analogous to Theorem 13.19 for Fourier transforms, but with one important caveat. Since all functions must vanish for t < 0, we are only allowed to shift them to the right, a shift to the left would produce nonzero function values for some t < 0. In general, the Laplace transform of the function f (t - c) shifted to the right by an amount c > 0, is L[ t2 ] = L[ f (t - c) ] =
0 0 t
L[ f (n) (t) ] = sn F (s) - sn-1 f (0) - sn-2 f (0) - - f (n-1) (0).
L[ t f (t) ] = - F (s).
(13.147)
f (t - c) e- s t dt = dt +
0
-c
f (t) e- s (t+c) dt dt = e
-sc 0
(13.150) f (t) e- s t dt = e- s c F (s).
=
-c
f (t) e
- s (t+c)
f (t) e
- s (t+c)
In this computation, we first used a change of variables in the integral, replacing t - c by t; then, the fact that f (t) 0 for t < 0 was used to eliminate the integral from - c to 0. When using the shift formula in practice, it is important to keep in mind that c > 0 and the function f (t - c) vanishes for all t < c. In the table, the factor (t - c) is used to remind the user of this fact. 2/25/07 747
c 2006 Peter J. Olver
1, b < t < c, for some 0, otherwise, 0 < b < c. To compute its Laplace transform, we write it as the difference Example 13.37. Consider the square wave pulse f (t) = f (t) = (t - b) - (t - c) of shifted versions of the step function (13.139). Combining the shift formula (13.150) and the formula (13.140) for the Laplace transform of the step function, we find L[ f (t) ] = L[ (t - b) ] - L[ (t - c) ] = e- s b - e - s c . s (13.151)
We already noted that the Fourier transform of the convolution product of two functions is realized as the ordinary product of their individual transforms. A similar result holds for the Laplace transform. Let f (t), g(t) be given functions. Since we are implicitly assuming that the functions vanish at all negative values of t, their convolution product (13.127) reduces to a finite integral
t
h(t) = f (t) g(t) =
0
f (t - ) g( ) d.
(13.152)
In particular h(t) = 0 for all t < 0 also. Further, it is not hard to show that the convolution of two functions of exponential growth also has exponential growth. Theorem 13.38. If L[ f (t) ] = F (s) and L[ g(t) ] = G(s), then the convolution h(t) = f (t) g(t) has Laplace transform given by the product H(s) = F (s) G(s). The proof of the convolution theorem for the Laplace transform proceeds along the same lines as its Fourier transform version Theorem 13.26, and is left as Exercise for the reader. Applications to Initial Value Problems The key application of the Laplace transform is to facilitate the solution of initial value problems for linear, constant coefficient ordinary differential equations. As a prototypical example, consider the second order initial value problem a du d2 u +b + c u = f (t), 2 dt dt u(0) = , du (0) = , dt (13.153)
in which a, b, c are constant. We will solve the initial value problem by applying the Laplace transform to both sides of the differential equation. In view of the differentiation formulae (13.144, 145), a s2 L[ u(t) ] - s u(0) - u(0) + b s L[ u(t) ] - u(0) + c L[ u(t) ] = L[ f (t) ]. Setting L[ u(t) ] = U (s) and L[ f (t) ] = F (s), and making use of the initial conditions, the preceding equation takes the form (a s2 + b s + c) U (s) = F (s) + (a s + b) + a . 2/25/07 748
c 2006
(13.154)
Peter J. Olver
Thus, by applying the Laplace transform, we have effectively reduced the entire initial value problem to a single elementary algebraic equation! Solving for U (s) = F (s) + (a s + b) + a , a s2 + b s + c (13.155)
we then recover solution u(t) to the initial value problem by finding the inverse Laplace transform of U (s). As noted earlier, in practice the inverse transform is found by suitably massaging the answer (13.155) to be a combination of known transforms. Example 13.39. Consider the initial value problem u + u = 10 e- 3 t , u(0) = 1, u(0) = 2.
Taking the Laplace transform of the differential equation as above, we find (s2 + 1) U (s) - s - 2 = 10 , s+3 and so U (s) = s+2 10 + . s2 + 1 (s + 3)(s2 + 1)
The second summand does not directly correspond to any of the entries in our table of Laplace transforms. However, we can use the method of partial fractions to write it as a sum s+2 1 3-s 1 5 U (s) = 2 + + 2 = + 2 s +1 s+3 s +1 s+3 s +1 of terms appearing in the table. We conclude that the solution to our initial value problem is u(t) = e- 3 t + 5 sin t. Of course, the last example is a problem that you can easily solve directly. The standard method learned in your first course on differential equations is just as effective in finding the final solution, and does not require all the extra Laplace transform machinery! The Laplace transform method is, however, particularly effective when dealing with complications that arise in cases of discontinuous forcing functions. Example 13.40. Consider a mass vibrating on a spring with fixed stiffness c = 4. Assume that the mass starts at rest, is then subjected to a unit force over time interval 1 2 < t < 2 , after which it left to vibrate on its own. The initial value problem is d2 u + 4u = dt2 1, 0,
1 2
< t < 2 ,
otherwise,
u(0) = u(0) = 0.
Taking the Laplace transform of the differential equation, and using (13.151), we find (s2 + 4) U (s) = e- s/2 - e-2 s , s and so U (s) = e- s/2 - e-2 s . s(s2 + 4)
Therefore, by the shift formula (13.150)
1 u(t) = h t - 2 - h(t - 2 ),
2/25/07
749
c 2006
Peter J. Olver
where h(t) is the function with Laplace transform L[ h(t) ] = H(s) = s(s2 1 1 = + 4) 4 s 1 - 2 s s +4 ,
which has been conveniently rewritten using partial fractions. Referring to our table of Laplace transforms, 1 1 t > 0, 4 - 4 cos 2t, h(t) = 0, t < 0. Therefore, our desired solution is 0, 1 1 u(t) = 4 + cos 2t, 1 4 2 cos 2t,
2 t.
1 0 t 2 , 1 2 t 2 ,
Note that the solution u(t) is only C1 at the points of discontinuity of the forcing function. Remark : A direct solution of this problem would proceed as follows. One solves the differential equation on each interval of continuity of the forcing function, leading to a solution on that interval depending upon two integration constants. The integration constants are then adjusted so that the solution satisfies the initial conditions and is continuously differentiable at each point of discontinuity of the forcing function. The details are straightforward, but messy. The Laplace transform method successfully bypasses these intervening manipulations. Example 13.41. Consider the initial value problem d2 u + 2 u = f (t), dt2 u(0) = u(0) = 0,
involving a general forcing function f (t). Applying the Laplace transform to the differential equation and using the initial conditions, (s2 + 2 ) U (s) = F (s), and hence U (s) = s2 F (s) . + 2
The right hand side is the product of the Laplace transform of the forcing function f (t) sin t . Therefore, Theorem 13.38 implies and that of the trigonometric function k(t) = that the solution can be written as their convolution
t t
u(t) = f k(t) = where
0
k(t - ) f ( ) d = sin t , 0, 750
0
sin (t - ) f ( ) d. t > 0, t < 0.
c 2006
(13.156)
k(t) = 2/25/07
Peter J. Olver
The integral kernel k(t - ) is known as the fundamental solution to the initial value problem, and prescribes the response of the system to a unit impulse force that is applied instantaneously at the time t = . Note particularly that (unlike boundary value problems) the impulse only affect the solutions at later times t > . For initial value problems, the fundamental solution plays a role similar to that of a Green's function in a boundary value problem. The convolution formula (13.156) can be viewed as a linear superposition of fundamental solution responses induced by expressing the forcing function f (t) =
0
f ( ) (t - ) d
as a superposition of concentrated impulses over time. This concludes our brief introduction to the Laplace transform and a few of its many applications to physical problems. More details can be found in almost all applied texts on mechanics, electronic circuits, signal processing, control theory, and many other areas, [Ltr].
2/25/07
751
c 2006
Peter J. Olver