fa

Chapter 13 FourierAnalysis In addition to their inestimable importance in mathematics and its applications, Fourier series also serve as the entry point into the wonderf ul world of Fourier analy- sis and its wide-ranging extensions and generalizations. A n entire industry is devoted to further developing the theory and enlarging the scope of app lications of Fourier–inspired methods. New directions in Fourier analysis continue to be d iscovered and exploited in a broad range of physical, mathematical, engineering, chem ical, biological, financial, and other systems. In this chapter, we will concentrate on four o f 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 ad- dition, more general types of eigenfunction expansions ass ociated 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 base d on discrete data, not continuous functions. One typically samples an analog s ignal at equally spaced time intervals, and then works exclusively with the resulting di screte (digital) data. The asso- ciated discrete Fourier representation re-expresses the d ata in terms of sampled complex exponentials; it can, in fact, be handled by finite-dimensio nal 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 understand- ing and analyzing its discrete digital counterpart. An impo rtant application of discrete Fourier sums is in signal and image processing. Basic data co mpression and noise removal algorithms are applied to the sample’s discrete Fourier coe fficients, acting on the obser- vation that noise tends to accumulate in the high frequency F ourier modes, while most important features are concentrated at low frequencies. Th e first Section 13.1 develops the basic Fourier theory in this discrete setting, culminat ing in the Fast Fourier Transform (FFT), which produces an efficient numerical algorithm for pa ssing between a signal and its discrete Fourier coefficients. One of the inherent limitations of classical Fourier method s, both continuous and discrete, is that they are not well adapted to localized data . (In physics, this lack of localizationisthebasisoftheHeisenbergUncertaintyPri nciple.) Asaresult, Fourier-based signal processing algorithms tend to be inaccurate and/or i nefficient when confronting highly localized signals or images. In the second section, w e introduce the modern theory ofwavelets, which isa recent extensionof Fourier analysis that more naturallyincorporates multiple scales and localization. Wavelets are playing an i ncreasingly dominant role in many modern applications; for instance, the new JPEG digita l image compression format is based on wavelets, as are the computerized FBI fingerprint data used in law enforcement 12/11/12 678 c/circleco√yrt2012 Peter J. Olver in the United States. Spectral analysis of non-periodic functions defined on the e ntire real line requires re- placing the Fourier series by a limiting Fourier integral. T he resulting Fourier transform plays an essential role in functional analysis, in the analy sis of ordinary and partial differ- ential equations, and in quantum mechanics, data analysis, signal processing, and many other applied fields. In Section 13.3, we introduce the most i mportant applied features of the Fourier transform. The closely related Laplace transform is a basic tool in engi neering applications. To mathematicians, the Fourier transform is the more fundamen tal of the two, while the Laplace transform is viewed as a certain real specializatio n. Both transforms change differ- entiation into multiplication, thereby converting linear differential equations into algebraic equations. The Fourier transform is primarily used for solv ing boundary value problems on the real line, while initial value problems, particularl y those involving discontinuous forcing terms, are effectively handled by the Laplace transf orm. 13.1. Discrete Fourier Analysis and the Fast Fourier Transf orm. In modern digital media — audio, still images or video — conti nuous 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, upo n which all of modern signal processing is based, is that the sampled complex exponentia ls 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 comp uting the discrete Fourier representation and reconstructing the signal from its Four ier coefficients. We will concentrate on the one-dimensional version here. Le tf(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 xj=a+jh, j = 0,...,n, where h=b−a n indicates the sample rate. In signal processing applicatio ns,xrepresents time instead of space, and the xjare the times at which we sample the signal f(x). Sample rates can be very high, e.g., every 10–20 milliseconds in current speech recognition systems. For simplicity, we adopt the “standard” interval of 0 ≤x≤2π, and the nequally spaced sample points† x0= 0, x1=2π n, x2=4π n, ... xj=2jπ n, ... xn−1=2(n−1)π n.(13.1) (Signals defined on other intervals can be handled by simply r escaling the interval to have length 2 π.) Sampling a (complex-valued) signal or function f(x) produces the sample †We will find it convenient to omit the final sample point xn= 2πfrom consideration. 12/11/12 679 c/circleco√yrt2012 Peter J. Olver 1 2 3 4 5 6 -1-0.50.51 1 2 3 4 5 6 -1-0.50.51 1 2 3 4 5 6 -1-0.50.51 1 2 3 4 5 6 -1-0.50.51 1 2 3 4 5 6 -1-0.50.51 1 2 3 4 5 6 -1-0.50.51 Figure 13.1. Sampling e−ixande7ixonn= 8 sample points. vector f=/parenleftbig f0,f1,...,fn−1/parenrightbigT=/parenleftbig f(x0),f(x1),...,f(xn−1)/parenrightbigT, where fj=f(xj) =f/parenleftbigg2jπ n/parenrightbigg . (13.2) Sampling cannot distinguish between functions that have th e same values at all of the sample points — from the sampler’s point of view they are iden tical. For example, the periodic complex exponential function f(x) =einx= cosnx+ i sinnx has sampled values fj=f/parenleftbigg2jπ n/parenrightbigg = exp/parenleftbigg in2jπ n/parenrightbigg =e2jπi= 1 for all j= 0,...,n−1, and hence is indistinguishable from the constant function c(x)≡1 — both lead to the samesample vector (1 ,1,...,1)T. This has the important implication that sampling at n equallyspaced samplepoints cannotdetect periodicsignalsoffrequency n. Moregenerally, the two complex exponential signals ei(k+n)xand eikx are also indistinguishable when sampled. This has the impor tant consequence that we need only use the first nperiodic complex exponential functions f0(x) = 1, f1(x) =eix, f2(x) =e2ix, ... fn−1(x) =e(n−1)ix,(13.3) in order to represent any 2 πperiodic sampled signal. In particular, exponentials e−ikxof “negative” frequency can all be converted into positive ver sions, namely ei(n−k)x, by the same sampling argument. For example, e−ix= cosx−i sinxande(n−1)ix= cos(n−1)x+ i sin(n−1)x 12/11/12 680 c/circleco√yrt2012 Peter J. Olver have identical values on the sample points (13.1). However, off of the sample points, they arequitedifferent; the former isslowlyvarying, while thel atter represents a high frequency oscillation. In Figure 13.1, we compare e−ixande7ixwhen there are n= 8 sample values, indicated by thedots on thegraphs. The top row compares the r eal parts, cos xand cos7 x, while the bottom row compares the imaginary parts, sin xand−sin7x. Note that both functions have the same pattern of sample values, even thoug h 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 imp osed by the frames of the film as an equivalently discretized low frequency motion in r everse. 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 ad equately represented. In fact, human appreciation of music also relies on inaudible high fr equency 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 t hat it was not set high enough to fully reproduce the musical quality of an analog LP record ! Thediscrete Fourier representation decomposes a sampled function f(x) into a linear combination of complex exponentials. Since we cannot disti nguish sampled exponentials of frequency higher than n, we only need consider a finite linear combination f(x)∼p(x) =c0+c1eix+c2e2ix+···+cn−1e(n−1)ix=n−1/summationdisplay k=0ckeikx(13.4) of the first nexponentials (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)canbeviewedasa(complex-valued) interpolating trigonometric polynomial of degree ≤n−1 for the sample data fj=f(xj). Remark: Iff(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 st ate 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 reta in the complex exponential form (13.4) of the discrete Fourier sum. †In computer graphics, the term “aliasing” is used in a much broader sens e that covers a variety of artifacts introduced by discretization — particularly, t he jagged appearance of lines and smooth curves on a digital monitor. 12/11/12 681 c/circleco√yrt2012 Peter J. Olver Since we are working in the finite-dimensional vector space Cnthroughout, we may reformulate the discrete Fourier series in vectorial form. Sampling the basic exponentials (13.3) produces the complex vectors ωk=/parenleftbig eikx0,eikx1,eikx2,...,eikxn−1/parenrightbigT =/parenleftig 1,e2kπi/n,e4kπi/n,...,e2(n−1)kπi/n/parenrightigT ,k= 0,...,n−1.(13.6) The interpolation conditions (13.5) can be recast in the equ ivalent vector form f=c0ω0+c1ω1+···+cn−1ωn−1. (13.7) In other words, to compute the discrete Fourier coefficients c0,...,cn−1off, all we need to do is rewrite its sample vector fas a linear combination of the sampled exponential vectorsω0,...,ωn−1. Now, aswithcontinuous Fourier series, theabsolutelycruc ial property istheorthonor- mality of the basis elements ω0,...,ωn−1. Were it not for the power of orthogonality, Fourier analysis might have remained a mere mathematical cu riosity, rather than today’s indispensable tool. Proposition 13.1. Thesampledexponentialvectors ω0,...,ωn−1formanorthonor- mal basis of Cnwith respect to the inner product /angbracketleftf;g/angbracketright=1 nn−1/summationdisplay j=0fjgj=1 nn−1/summationdisplay j=0f(xj)g(xj),f,g∈Cn.(13.8) Remark: The inner product (13.8) is a rescaled version of the standa rd Hermitian dot product (3.92) between complex vectors. We can interpret th e inner product between the sample vectors f,gas theaverageof 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 probl em, so their discrete sampled counterparts are eigenvectors for a self-adjoint matrix ei genvalue problem; details can be found in Exercise . Here, though, to keep the discussion on track, we shall outl ine a direct proof. Proof: The crux of the matter relies on properties of the remarkabl e complex numbers ζn=e2πi/n= cos2π n+ i sin2π n,where n= 1,2,3,... . (13.9) Particular cases include ζ2=−1, ζ3=−1 2+√ 3 2i, ζ4= i,andζ8=√ 2 2+√ 2 2i.(13.10) Thenthpower of ζnis ζn n=/parenleftig e2πi/n/parenrightig n=e2πi= 1, 12/11/12 682 c/circleco√yrt2012 Peter J. Olver ζ0 5= 1ζ5 ζ2 5 ζ3 5 ζ4 5 Figure 13.2. The Fifth Roots of Unity. and hence ζnis one of the complex nthroots of unity :ζn=n√ 1. There are, in fact, n different complex nthroots of 1, including 1 itself, namely the powers of ζn: ζk n=e2kπi/n= cos2kπ n+ i sin2kπ n, k = 0,...,n−1. (13.11) Since it generates all the others, ζnis known as the primitive nthroot of unity . Geometri- cally, the nthroots (13.11) lie on the vertices of a regular unit n–gon in the complex plane; see Figure 13.2. The primitive root ζnis the first vertex we encounter as we go around the n–gon in a counterclockwise direction, starting at 1. Contin uing around, the other roots appear in their natural order ζ2 n,ζ3 n,...,ζn−1 n, and finishing back at ζn n= 1. The complex conjugate of ζnis the “last” nthroot e−2πi/n=ζn=1 ζn=ζn−1 n=e2(n−1)πi/n. (13.12) The complex numbers (13.11) are a complete set of roots of the polynomial zn−1, which can therefore be factored: zn−1 = (z−1)(z−ζn)(z−ζ2 n)···(z−ζn−1 n). On the other hand, elementary algebra provides us with the re al factorization zn−1 = (z−1)(1+z+z2+···+zn−1). Comparing the two, we conclude that 1+z+z2+···+zn−1= (z−ζn)(z−ζ2 n)···(z−ζn−1 n). Substituting z=ζk ninto both sides of this identity, we deduce the useful formul a 1+ζk n+ζ2k n+···+ζ(n−1)k n=/braceleftbiggn, k= 0, 0,0< k < n.(13.13) Sinceζn+k n=ζk n, this formula can easily be extended to general integers k; the sum is equal to nifnevenly divides kand is 0 otherwise. 12/11/12 683 c/circleco√yrt2012 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 w ritten in terms of the nthroots of unity: ωk=/parenleftbig 1,ζk n,ζ2k n,ζ3k n, ... ,ζ(n−1)k n/parenrightbigT, k = 0,...,n−1. (13.14) Therefore, applying (13.12,13), we conclude that /angbracketleftωk;ωl/angbracketright=1 nn−1/summationdisplay j=0ζjk nζjl n=1 nn−1/summationdisplay j=0ζj(k−l) n=/braceleftigg 1, k=l, 0, k/negationslash=l,0≤k,l < n, which establishes orthonormality of the sampled exponenti al vectors. Q.E.D. Orthonormality of the basis vectors implies that we can imme diately compute the Fourier coefficients in the discrete Fourier sum (13.4) by tak ing inner products: ck=/angbracketleftf;ωk/angbracketright=1 nn−1/summationdisplay j=0fjeikxj=1 nn−1/summationdisplay j=0fje−ikxj=1 nn−1/summationdisplay j=0ζ−jk nfj.(13.15) In other words, the discrete Fourier coefficient ckis obtained by averaging the sampled val- ues ofthe product function f(x)e−ikx. Thepassage from a signal toitsFourier coefficients is known as the Discrete Fourier Transform or DFT for short. The reverse procedure of reconstructing a signal from its discrete Fourier coefficien ts via the sum (13.4) (or (13.7)) is known as the Inverse Discrete Fourier Transform or IDFT. The Discrete Fourier Trans- form and its inverse define mutually inverse linear transfor mations on the space Cn, whose matrix representations can be found in Exercise . Example 13.2. Ifn= 4, then ζ4= i. The corresponding sampled exponential vectors ω0= 1 1 1 1 ,ω1= 1 i −1 −i ,ω2= 1 −1 1 −1 ,ω3= 1 −i −1 i , form an orthonormal basis of C4with respect to the averaged Hermitian dot product /angbracketleftv;w/angbracketright=1 4/parenleftbig v0w0+v1w1+v2w2+v3w3/parenrightbig ,where v= v0 v1 v2 v3 ,w= w0 w1 w2 w3 . Given the sampled function values f0=f(0), f1=f/parenleftbig1 2π/parenrightbig , f2=f(π), f3=f/parenleftbig3 2π/parenrightbig , we construct the discrete Fourier representation f=c0ω0+c1ω1+c2ω2+c3ω3, (13.16) where c0=/angbracketleftf;ω0/angbracketright=1 4(f0+f1+f2+f3), c1=/angbracketleftf;ω1/angbracketright=1 4(f0−if1−f2+ if3), c2=/angbracketleftf;ω2/angbracketright=1 4(f0−f1+f2−f3), c3=/angbracketleftf;ω3/angbracketright=1 4(f0+ if1−f2−if3). 12/11/12 684 c/circleco√yrt2012 Peter J. Olver 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 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+c1eix+c2e2ix+c3e3ix that agrees with f(x) on the sample points. For instance, if f(x) = 2πx−x2, then f0= 0., f1= 7.4022, f2= 9.8696, f3= 7.4022, and hence c0= 6.1685, c1=−2.4674, c2=−1.2337, c3=−2.4674. Therefore, the interpolating trigonometric polynomial is given by the real part of p(x) = 6.1685−2.4674eix−1.2337e2ix−2.4674e3ix, (13.17) namely, Rep(x) = 6.1685−2.4674 cos x−1.2337 cos2 x−2.4674 cos3 x. (13.18) In Figure 13.3 we compare the function, with the interpolati on points indicated, and dis- crete Fourier representations (13.18) for both n= 4 and n= 16 points. The resulting graphs point out a significant difficulty with the Discrete Fou rier Transform as developed so far. While the trigonometric polynomials do indeed corre ctly match the sampled func- tion values, their pronounced oscillatory behavior makes t hem 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 freque ncies that are represented in the Fourier sum. Indeed, the graphs in Figure 13.3 might re mind you of our earlier observation that, due to aliasing, low and high frequency ex ponentials can have the same sample data, but differ wildly in between the sample points. W hile the first half of the 12/11/12 685 c/circleco√yrt2012 Peter J. Olver 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 1 2 3 4 5 6246810 Figure 13.4. The Low Frequency Discrete Fourier Representation of x2−2πx. summands in (13.4) represent relatively low frequencies, t he second half do not, and can be replaced by equivalent lower frequency, and hence less osci llatory exponentials. Namely, if 0< k≤1 2n, thene−ikxandei(n−k)xhave the same sample values, but the former is of lower frequency than the latter. Thus, for interpolatory pu rposes, we should replace the second half of the summands in the Fourier sum (13.4) by their low frequency alternatives. Ifn= 2m+1 is odd, then we take /hatwidep(x) =c−me−imx+···+c−1e−ix+c0+c1eix+···+cmeimx=m/summationdisplay k=−mckeikx(13.19) as the equivalent low frequency interpolant. If n= 2mis even — which is the most common case occurring in applications — then /hatwidep(x) =c−me−imx+···+c−1e−ix+c0+c1eix+···+cm−1ei(m−1)x=m−1/summationdisplay k=−mckeikx (13.20) will be our choice. (It is a matter of personal taste whether t o usee−imxoreimxto represent the highest frequency term.) In both cases, the Fo urier coefficients with negative indices are the same as their high frequency alternatives: c−k=cn−k=/angbracketleftf;ωn−k/angbracketright=/angbracketleftf;ω−k/angbracketright, (13.21) whereω−k=ωn−kis the sample vector for e−ikx∼ei(n−k)x. Returning to the previous example, for interpolating purpo ses, we should replace (13.17) by the equivalent low frequency interpolant /hatwidep(x) =−1.2337e−2ix−2.4674e−ix+6.1685−2.4674eix, (13.22) with real part Re/hatwidep(x) = 6.1685−4.9348 cos x−1.2337 cos2 x. Graphs of the n= 4 and 16 low frequency trigonometric interpolants can be se en in Figure 13.4. Thus, by utilizing only the lowest frequency ex ponentials, we successfully 12/11/12 686 c/circleco√yrt2012 Peter J. Olver suppress the aliasing artifacts, resulting in a quite reaso nable trigonometric interpolant to the given function. Remark: The low frequency version also serves to unravel the realit y 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−ke−ikx+ckeikx=akcoskx+bksinkx add up to a real trigonometric function. Therefore, the odd ninterpolant (13.19) is a real trigonometric polynomial, whereas in the even version (13. 20) only the highest frequency termc−me−imxproduces a complex term — which is, in fact, 0 on the sample poi nts. Compression and Noise Removal In a typical experimental signal, noise primarily affects th e high frequency modes, while the authentic features tend to appear in the low freque ncies. Think of the hiss and static you hear on an AM radio station or a low quality audio re cording. Thus, a very simple, but effective, method for denoising a corrupted sign al is to decompose it into its Fourier modes, as in (13.4), and then discard the high freque ncy constituents. A similar idea underlies the Dolby/circleco√†rtTrecording system used on most movie soundtracks: during the recording process, the high frequency modes are artifici ally boosted, so that scaling them back when showing the movie in the theater has the effect o f eliminating much of the extraneous noise. The one design issue is the specificati on of a cut-off between low and high frequency, that is, between signal and noise. This c hoice will depend upon the properties of the measured signal, and is left to the discret ion of the signal processor. A correct implementation of the denoising procedure is faci litated by using the un- aliased 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 ql(x) =l/summationdisplay k=−lckeikx, (13.23) wherel <1 2(n+1) specifies the selected cut-off frequency between signal a nd noise. The 2l+1≪nlow frequency Fourier modes retained in (13.23) will, in fav orable 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= 29= 512 sample points in the discrete Fourier repre- sentation, 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−5e−5ix+···+c5e5ix. To compare, we plot both the original signal and the denoised version on the same grap h. In this case, the maximal deviation is less than .15 over the entire interval [0 ,2π]. 12/11/12 687 c/circleco√yrt2012 Peter J. Olver 1 2 3 4 5 61234567 The Original Signal1 2 3 4 5 61234567 The Noisy Signal 1 2 3 4 5 61234567 The Denoised Signal1 2 3 4 5 61234567 Comparison of the Two Figure 13.5. Denoising a Signal. 1 2 3 4 5 61234567 The Original Signal1 2 3 4 5 61234567 Moderate Compression1 2 3 4 5 61234567 High Compression Figure 13.6. Compressing a Signal. The same idea underlies many data compression algorithms fo r audio recordings, digi- tal images and, particularly, video. The goal is efficient sto rage 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 wi ll, in favorable situations, not lead to any noticeable degradation of the signal or image. Th us, to compress a signal (and, simultaneously, remove high frequency noise), we ret ain only its low frequency dis- crete Fourier coefficients. The signal is reconstructed by su mming the associated truncated discrete Fourier series (13.23). A mathematical justificat ion 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 importa nce. In Figure 13.6, the same signal is compressed by retaining, r espectively, 2 l+1 = 21 and 2l+1 = 7 Fourier coefficients only instead of all n= 512 that would be required for complete accuracy. For the case of moderate compression, th e maximal deviation between the signal and the compressed version is less than 1 .5×10−4over the entire interval, 12/11/12 688 c/circleco√yrt2012 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 detail ed the original signal, the more Fourier modes need to be retained for accurate reproduc tion. The Fast Fourier Transform While one may admire an algorithm for its intrinsic beauty, i n the real world, the bottom line is always efficiency of implementation: the less t otal computation, the faster the processing, and hence the more extensive the range of app lications. 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 Cool ey and John Tukey, [ 44], discovered†a much more efficient approach to the Discrete Fourier Transfo rm, exploiting the rather special structure of the sampled exponential vec tors. 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 , [29,30]. Ingeneral, computingallthediscreteFouriercoefficients( 13.15)ofan ntimessampled signal requires a total of n2complex multiplications and n2−ncomplex additions. Note also that each complex addition z+w= (x+ iy)+(u+ iv) = (x+u)+ i(y+v) (13 .24) generally requires two real additions, while each complex m ultiplication zw= (x+ iy)(u+ iv) = (xu−yv)+ i(xv+yu) (13 .25) requires4realmultiplicationsand2realadditions, or, by employingthealternativeformula xv+yu= (x+y)(u+v)−xu−yv (13.26) for the imaginary part, 3 real multiplications and 5 real add itions. (The choice of formula (13.25) or (13.26) will depend upon the processor’s relativ e speeds of multiplication and addition.) Similarly, given the Fourier coefficients c0,...,cn−1, reconstruction of the sam- pled signal via(13.4) requires n2−ncomplex multiplications and n2−ncomplex additions. As a result, both computations become quite labor intensive for large n. Extending these ideas to multi-dimensional data only exacerbates the probl em. In order to explain the method without undue complication, w e 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 frequ ency 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 earl y 1800’s, but his insight was not fully appreciated until modern computers arrived on t he scene. 12/11/12 689 c/circleco√yrt2012 Peter J. Olver is even, then the primitive mthroot of unity ζm=m√ 1 equals the square of the primitive nthroot: ζm=ζ2 n. We use this fact to split the summation (13.15) for the order ndiscrete Fourier coefficients into two parts, collecting together the even and the odd powe rs ofζk n: ck=1 n/parenleftbig f0+f1ζ−k n+f2ζ−2k n+···+fn−1ζ−(n−1)k n/parenrightbig =1 n/parenleftbig f0+f2ζ−2k n+f4ζ−4k n+···+f2m−2ζ−(2m−2)k n/parenrightbig + +ζ−k n1 n/parenleftbig f1+f3ζ−2k n+f5ζ−4k n+···+f2m−1ζ−(2m−2)k n/parenrightbig =1 2/braceleftbigg1 m/parenleftbig f0+f2ζ−k m+f4ζ−2k m+···+f2m−2ζ−(m−1)k m/parenrightbig/bracerightbigg + +ζ−k n 2/braceleftbigg1 m/parenleftbig f1+f3ζ−k m+f5ζ−2k m+···+f2m−1ζ−(m−1)k m/parenrightbig/bracerightbigg .(13.27) Now, observe that the expressions in braces are the order mFourier coefficients for the sample data fe=/parenleftbig f0,f2,f4,...,f2m−2/parenrightbigT=/parenleftbig f(x0),f(x2),f(x4),...,f(x2m−2)/parenrightbigT, fo=/parenleftbig f1,f3,f5,...,f2m−1/parenrightbigT=/parenleftbig f(x1),f(x3),f(x5),...,f(x2m−1)/parenrightbigT.(13.28) Note that feis obtained by sampling f(x) on the evensample points x2j, whilefois obtained by sampling the same function f(x), but now at the oddsample points x2j+1. In other words, we are splitting the original sampled signal in to two “half-sampled” signals obtained by sampling on every other point. The even and odd Fo urier coefficients are ce k=1 m/parenleftbig f0+f2ζ−k m+f4ζ−2k m+···+f2m−2ζ−(m−1)k m/parenrightbig , co k=1 m/parenleftbig f1+f3ζ−k m+f5ζ−2k m+···+f2m−1ζ−(m−1)k m/parenrightbig ,k= 0,...,m−1. (13.29) Since they contain just mdata values, both the even and odd samples require only m distinct Fourier coefficients, and we adopt the identificatio n ce k+m=ce k, co k+m=co k, k = 0,...,m−1. (13.30) Therefore, the order n= 2mdiscrete Fourier coefficients (13.27) can be constructed fro m a pair of order mdiscrete Fourier coefficients via ck=1 2/parenleftbig ce k+ζ−k nco k/parenrightbig , k = 0,...,n−1. (13.31) Now ifm= 2lis also even, then we can play the same game on the order mFourier coefficients (13.29), reconstructing each of them from a pair of order ldiscrete Fourier coefficients — obtained by sampling the signal at every fourth point. If n= 2ris 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 sample s the function at a single point. 12/11/12 690 c/circleco√yrt2012 Peter J. Olver The result is the desired algorithm. After some rearrangeme nt of the basic steps, we arrive at the Fast Fourier Transform, which we now present in its fina l form. We begin with a sampled signal on n= 2rsample points. To efficiently program the Fast Fourier Transform, it helps to write out each index 0 ≤j <2rin its binary (as opposed to decimal) representation j=jr−1jr−2... j2j1j0,where jν= 0 or 1; (13 .32) the notation is shorthand for its rdigit binary expansion j=j0+2j1+4j2+8j3+···+2r−1jr−1. We then define the bit reversal map ρ(jr−1jr−2... j2j1j0) =j0j1j2... jr−2jr−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 especia lly that the bit reversal map ρ=ρrdepends upon the original choice of r= log2n. Secondly, for each 0 ≤k < r, define the maps αk(j) =jr−1... jk+10jk−1... j0, βk(j) =jr−1... jk+11jk−1... j0=αk(j)+2k,forj=jr−1jr−2... j1j0. (13.34) In other words, αk(j) sets the kthbinary digit of jto 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 ea sy to implement on modern binary computers. Given a sampled signal f0,...,fn−1, its discrete Fourier coefficients c0,...,cn−1are computed by the following iterative algorithm: c(0) j=fρ(j), c(k+1) j=1 2/parenleftbig c(k) αk(j)+ζ−j 2k+1c(k) βk(j)/parenrightbig ,j= 0,...,n−1, k= 0,...,r−1,(13.35) in which ζ2k+1is the primitive 2k+1root of unity. The final output of the iterative proce- dure, namely cj=c(r) j, j = 0,...,n−1, (13.36) are the discrete Fourier coefficients of our signal. The prepr ocessing step of the algorithm, where we define c(0) j, produces a more convenient rearrangement of the sample val ues. The subsequent steps successively combine the Fourier coeffi cients of the appropriate even and odd sampled subsignals together, reproducing (13.27) i n a different notation. The following example should help make the overall process clea rer. Example 13.3. Consider the case r= 3, and so our signal has n= 23= 8 sampled valuesf0,f1,...,f7. We begin the process by rearranging the sample values c(0) 0=f0, c(0) 1=f4, c(0) 2=f2, c(0) 3=f6, c(0) 4=f1, c(0) 5=f5, c(0) 6=f3, c(0) 7=f7, 12/11/12 691 c/circleco√yrt2012 Peter J. Olver 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 c(1) 0=1 2(c(0) 0+c(0) 1), c(1) 1=1 2(c(0) 0−c(0) 1), c(1) 2=1 2(c(0) 2+c(0) 3), c(1) 3=1 2(c(0) 2−c(0) 3), c(1) 4=1 2(c(0) 4+c(0) 5), c(1) 5=1 2(c(0) 4−c(0) 5), c(1) 6=1 2(c(0) 6+c(0) 7), c(1) 7=1 2(c(0) 6−c(0) 7), where we combine successive pairs of the rearranged sample v alues. The second stage of the iteration has k= 1 with ζ4= i. We find c(2) 0=1 2(c(1) 0+c(1) 2), c(2) 1=1 2(c(1) 1−ic(1) 3), c(2) 2=1 2(c(1) 0−c(1) 2), c(2) 3=1 2(c(1) 1+ ic(1) 3), c(2) 4=1 2(c(1) 4+c(1) 6), c(2) 5=1 2(c(1) 5−ic(1) 7), c(2) 6=1 2(c(1) 4−c(1) 6), c(2) 7=1 2(c(1) 5+ ic(1) 7). Note that the indices of the combined pairs of coefficients diff er by 2. In the last step, wherek= 2 and ζ8=√ 2 2(1+ i), we combine coefficients whose indices differ by 4 = 22; the final output c0=c(3) 0=1 2(c(2) 0+c(2) 4), c4=c(3) 4=1 2(c(2) 0−c(2) 4), c1=c(3) 1=1 2/parenleftbig c(2) 1+√ 2 2(1−i)c(2) 5/parenrightbig , c5=c(3) 5=1 2/parenleftbig c(2) 1−√ 2 2(1−i)c(2) 5/parenrightbig , c2=c(3) 2=1 2/parenleftbig c(2) 2−ic(2) 6/parenrightbig , c6=c(3) 6=1 2/parenleftbig c(2) 2+ ic(2) 6/parenrightbig , c3=c(3) 3=1 2/parenleftbig c(2) 3−√ 2 2(1+ i)c(2) 7/parenrightbig , c7=c(3) 7=1 2/parenleftbig c(2) 3+√ 2 2(1+ i)c(2) 7/parenrightbig , is the complete set of discrete Fourier coefficients. Let us count the number of arithmetic operations required in the Fast Fourier Trans- form algorithm. At each stage in the computation, we must per formn= 2rcomplex additions/subtractions and the same number of complex mult iplications. (Actually, the number of multiplications is slightly smaller since multip lications by ±1 and±i are ex- tremely simple. However, this does not significantly alter t he final operations count.) There are r= log2nstages, and so we require a total of rn=nlog2ncomplex addi- tions/subtractions and the same number of multiplications . Now, when nis large,nlog2n issignificantly 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 nlog2n= 10,240 — a net savings of 99%. As a result, many large scale computati ons that would be in- tractable using the direct approach are immediately brough t into the realm of feasibility. This is the reason why all modern implementations of the Disc rete Fourier Transform are based on the FFT algorithm and its variants. The reconstruction of the signal from the discrete Fourier c oefficients c0,...,cn−1is speeded up in exactly the same manner. The only differences ar e that we replace ζ−1 n=ζn byζn, and drop the factors of1 2since there is no need to divide by nin the final result (13.4). Therefore, we apply the slightly modified iterative procedure f(0) j=cρ(j), f(k+1) j=f(k) αk(j)+ζj 2k+1f(k) βk(j),j= 0,...,n−1, k= 0,...,r−1,(13.37) 12/11/12 692 c/circleco√yrt2012 Peter J. Olver and finish with f(xj) =fj=f(r) j, j = 0,...,n−1. (13.38) Example 13.4. The reconstruction formulae in the case of n= 8 = 23Fourier coef- ficientsc0,...,c7, which were computed in Example 13.3, can be implemented as f ollows. First, we rearrange the Fourier coefficients in bit reversed o rder: f(0) 0=c0, f(0) 1=c4, f(0) 2=c2, f(0) 3=c6, f(0) 4=c1, f(0) 5=c5, f(0) 6=c3, f(0) 7=c7, Then we begin combining them in successive pairs: f(1) 0=f(0) 0+f(0) 1, f(1) 1=f(0) 0−f(0) 1, f(1) 2=f(0) 2+f(0) 3, f(1) 3=f(0) 2−f(0) 3, f(1) 4=f(0) 4+f(0) 5, f(1) 5=f(0) 4−f(0) 5, f(1) 6=f(0) 6+f(0) 7, f(1) 7=f(0) 6−f(0) 7. Next, f(2) 0=f(1) 0+f(1) 2, f(2) 1=f(1) 1+ if(1) 3, f(2) 2=f(1) 0−f(1) 2, f(2) 3=f(1) 1−if(1) 3, f(2) 4=f(1) 4+f(1) 6, f(2) 5=f(1) 5+ if(1) 7, f(2) 6=f(1) 4−f(1) 6, f(2) 7=f(1) 5−if(1) 7. Finally, the sampled signal values are f(x0) =f(3) 0=f(2) 0+f(2) 4, f (x4) =f(3) 4=f(2) 0−f(2) 4, f(x1) =f(3) 1=f(2) 1+√ 2 2(1+ i)f(2) 5, f (x5) =f(3) 5=f(2) 1−√ 2 2(1+ i)f(2) 5, f(x2) =f(3) 2=f(2) 2+ if(2) 6, f (x6) =f(3) 6=f(2) 2−if(2) 6, f(x3) =f(3) 3=f(2) 3−√ 2 2(1−i)f(2) 7, f (x7) =f(3) 7=f(2) 3+√ 2 2(1−i)f(2) 7. 13.2. Wavelets. Trigonometric Fourier series, both continuous and discret e, are amazingly powerful, but they do suffer from one potentially serious defect. The ba sis functions eikx= coskx+ i sinkxare spread out over the entire interval [ −π,π], and so are not well-suited to processing localized signals — meaning data that are concen trated 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 o f the Fourier trigonometric functions, but, in addition, adapts to localized structure s in signals. This dream was the inspiration for the development of the modern theory of wave lets. 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 mathemati cian Alfr´ ed Haar in 1910, [88]. We consider the space of functions (signals) defined the in terval [0,1], equipped with the standard L2inner product /angbracketleftf;g/angbracketright=/integraldisplay1 0f(x)g(x)dx. (13.39) 12/11/12 693 c/circleco√yrt2012 Peter J. Olver -0.2 0.2 0.4 0.6 0.8 1 1.2 -1-0.50.51 ϕ1(x)-0.2 0.2 0.4 0.6 0.8 1 1.2 -1-0.50.51 ϕ2(x) -0.2 0.2 0.4 0.6 0.8 1 1.2 -1-0.50.51 ϕ3(x)-0.2 0.2 0.4 0.6 0.8 1 1.2 -1-0.50.51 ϕ4(x) Figure 13.7. The First Four Haar Wavelets. This choice is merely for convenience, being slightly bette r suited to our construction than [−π,π]or [0,2π]. Moreover, the usual scaling arguments can beused to adapt the wavelet formulas to any other interval. TheHaar wavelets are certainpiecewise constant functions. The first four are graphed in Figure 13.7 The first is the box function ϕ1(x) =ϕ(x) =/braceleftbigg1,0< x≤1, 0,otherwise ,(13.40) 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 v alue, it will be more convenient to consistently choose the left hand limi ting value. The second Haar function ϕ2(x) =w(x) =  1,0< x≤1 2, −1,1 2< x≤1, 0,otherwise ,(13.41) is known as the mother wavelet . The third and fourth Haar functions are compressed 12/11/12 694 c/circleco√yrt2012 Peter J. Olver versions of the mother wavelet: ϕ3(x) =w(2x) =  1,0< x≤1 4, −1,1 4< x≤1 2, 0,otherwise ,ϕ4(x) =w(2x−1) =  1,1 2< x≤3 4, −1,3 4< x≤1, 0,otherwise , calleddaughter wavelets . One can easily check, by direct evaluation of the integrals , that thefourHaarwaveletfunctionsareorthogonalwithrespect totheL2innerproduct(13.39). The scaling transformation x/ma√sto→2xserves to compress the wavelet function, while the translation 2 x/ma√sto→2x−1 moves the compressed version to the right by a half a unit. Furthermore, we can represent the mother wavelet by compres sing and translating the scaling function,: w(x) =ϕ(2x)−ϕ(2x−1). (13.42) It is these two operations of scaling and compression — coupl ed with the all-important orthogonality — that underlies the power of wavelets. The Haar wavelets have an evident discretization. If we deco mpose the interval (0 ,1] into the four subintervals /parenleftbig 0,1 4/bracketrightbig ,/parenleftbig1 4,1 2/bracketrightbig ,/parenleftbig1 2,3 4/bracketrightbig ,/parenleftbig3 4,1/bracketrightbig , (13.43) on which the four wavelet functions are constant, then we can represent each of them by a vector in R4whose entries are the values of each wavelet function sample d at the left endpoint of each subinterval. In this manner, we obtain the w avelet sample vectors v1= 1 1 1 1 ,v2= 1 1 −1 −1 ,v3= 1 −1 0 0 ,v4= 0 0 1 −1 ,(13.44) that form the orthogonal wavelet basis of R4we encountered in Examples 2.35 and 5.10. Orthogonality of the vectors (13.44) with respect to the sta ndard Euclidean dot product is equivalent to orthogonality of the Haar wavelet functions w ith respect to the inner product (13.39). Indeed, if f(x)∼f= (f1,f2,f3,f4) and g(x)∼g= (g1,g2,g3,g4) arepiecewise constant real functions that achieve the indicated values on the four subin- tervals (13.43), then their L2inner product /angbracketleftf;g/angbracketright=/integraldisplay1 0f(x)g(x)dx=1 4/parenleftbig f1g1+f2g2+f3g3+f4g4/parenrightbig =1 4f·g, is equal to the averaged dot product of their sample values — t he real form of the inner product (13.8) that was used in the discrete Fourier transfo rm. Since the vectors (13.44) form an orthogonal basis of R4, we can uniquely decompose any such piecewise constant function as a linear combinatio n of wavelets f(x) =c1ϕ1(x)+c2ϕ2(x)+c3ϕ3(x)+c4ϕ4(x), 12/11/12 695 c/circleco√yrt2012 Peter J. Olver or, equivalently, in terms of the sample vectors, f=c1v1+c2v2+c3v3+c4v4. The required coefficients ck=/angbracketleftf;ϕk/angbracketright /bardblϕk/bardbl2=f·vk /bardblvk/bardbl2 are fixed by our usual orthogonality formula (5.7). Explicit ly, c1=1 4(f1+f2+f3+f4), c2=1 4(f1+f2−f3−f4),c3=1 2(f1−f2), c4=1 2(f3−f4). Before proceeding to the more general case, let us introduce an important analytical definition that quantifies precisely how localized a functio n is. Definition 13.5. Thesupportof a function f(x), written supp f, is the closure of the set where f(x)/negationslash= 0. Thus, a point will belong to the support of f(x), provided fis not zero there, or at least is not zero at nearby points. More precisely: Lemma 13.6. Iff(a)/negationslash= 0, thena∈suppf. More generally, a point a∈suppfif and only if there exist a convergent sequence xn→asuch that f(xn)/negationslash= 0. Conversely, a/negationslash∈suppfif and only if f(x)≡0on an interval a−δ < x < a +δfor some δ >0. Intuitively, thesmaller the support ofa function, themore localizeditis. Forexample, 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)/negationslash= 0 at nearby points. The two daughter wavelets have smaller support: suppϕ3=/bracketleftbig 0,1 2/bracketrightbig , suppϕ4=/bracketleftbig1 2,1/bracketrightbig , and so are twice as localized. An extreme case is the delta fun ction, whose support is a single point. In contrast, the support of the Fourier trigon ometric basis functions is all of R, since they only vanish at isolated points. The effect of scalings and translations on the support of a fun ction is easily discerned. Lemma 13.7. Ifsuppf= [a,b], and g(x) =f(rx−δ),then suppg=/bracketleftbigga+δ r,b+δ r/bracketrightbigg . In other words, scaling xby a factor rcompresses the support of the function by a factor 1/r, while translating xtranslates the support of the function. The key requirement for a wavelet basis is that it contains fu nctions with arbitrarily small support. To this end, the full Haar wavelet basis is obt ained from the mother wavelet by iterating the scaling and translation processes. We begi n with the scaling function ϕ(x), (13.45) 12/11/12 696 c/circleco√yrt2012 Peter J. Olver from which we construct the mother wavelet via (13.42). For a ny “generation” j≥0, we form the wavelet offspring by first compressing the mother wav elet so that its support fits into an interval of length 2−j, wj,0(x) =w(2jx),so that supp wj,0= [0,2−j], (13.46) and then translating wj,0so as to fill up the entire interval [0 ,1] by 2jsubintervals, each of length 2−j, defining wj,k(x) =wj,0(x−k) =w(2jx−k),where k= 0,1,...2j−1.(13.47) Lemma 13.7 implies that supp wj,k= [ 2−jk,2−j(k+1) ], and so the combined supports of all the jthgeneration of wavelets is the entire interval:2j−1/uniondisplay k=0suppwj,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 introduce d as ϕ3andϕ4, namely w1,0(x) =w(2x), w1,1(x) =w(2x−1). The second generation, j= 2, appends four additional granddaughter wavelets to our basis: w2,0(x) =w(4x), w2,1(x) =w(4x−1), w2,2(x) =w(4x−2), w2,3(x) =w(4x−3). The 8 Haar wavelets ϕ,w0,0,w1,0,w1,1,w2,0,w2,1,w2,2,w2,3are constant on the 8 subin- tervals of length1 8, taking the successive sample values indicated by the colum ns of the matrix W8= 1 1 1 0 1 0 0 0 1 1 1 0 −1 0 0 0 1 1 −1 0 0 1 0 0 1 1 −1 0 0 −1 0 0 1−1 0 1 0 0 1 0 1−1 0 1 0 0 −1 0 1−1 0 −1 0 0 0 1 1−1 0 −1 0 0 0 −1 . (13.48) Orthogonality of the wavelets is manifested in the orthogon ality of the columns of W8. (Unfortunately, the usual terminological constraints of D efinition 5.18 prevent us from callingW8an orthogonal matrix because its columns are not orthonorma l!) Thenthstage consists of 2n+1different wavelet functions comprising the scaling func- tions and all the generations up to the nth:w0(x) =ϕ(x) andwj,k(x) for 0≤j≤nand 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). 12/11/12 697 c/circleco√yrt2012 Peter J. Olver Proof: First, note that each wavelet wj,k(x) is equal to +1 on an interval of length 2−j−1and to−1 on an adjacent interval of the same length. Therefore, /angbracketleftwj,k;ϕ/angbracketright=/integraldisplay1 0wj,k(x)dx= 0, (13.49) since the +1 and −1 contributions cancel each other. If two different wavelets wj,kand wl,mwith, 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: /angbracketleftwj,k;wl,m/angbracketright=/integraldisplay1 0wj,k(x)wl,m(x)dx= 0. Otherwise, except in the case when the two wavelets are ident ical, the support of wl,mis entirely contained inan interval where wj,kis constant and so wj,k(x)wl,m(x) =±wl,m(x). Therefore, by (13.49), /angbracketleftwj,k;wl,m/angbracketright=/integraldisplay1 0wj,k(x)wl,m(x)dx=±/integraldisplay1 0wl,m(x)dx= 0. Finally, we compute /bardblϕ/bardbl2=/integraldisplay1 0dx= 1, /bardblwj,k/bardbl2=/integraldisplay1 0wj,k(x)2dx= 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, th ewavelet series of a signal f(x) is given by f(x)∼c0ϕ(x) +∞/summationdisplay j=02j−1/summationdisplay k=0cj,kwj,k(x). (13.51) Orthogonality implies that the wavelet coefficients c0,cj,kcan be immediately computed using the standard inner product formula coupled with (13.5 0): c0=/angbracketleftf;ϕ/angbracketright /bardblϕ/bardbl2=/integraldisplay1 0f(x)dx, cj,k=/angbracketleftf;wj,k/angbracketright /bardblwj,k/bardbl2= 2j/integraldisplay2−jk+2−j−1 2−jkf(x)dx−2j/integraldisplay2−j(k+1) 2−jk+2−j−1f(x)dx.(13.52) The convergence properties of the wavelet series (13.51) ar e 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 th e first plot. The next plots show the partial sums over j= 0,...,rwithr= 2,3,4,5,6. We have used a discontinuous signal to demonstrate that there is no n onuniform Gibbs phenomenon in a Haar wavelet expansion. Indeed, since the wavelets are t hemselves discontinuous, they do not have any difficult uniformly converging to a discontinu ous function. On the other hand, it takes quite a few wavelets to begin to accurately rep roduce 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 (exc luding points very close to the discontinuity). 12/11/12 698 c/circleco√yrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 0.2 0.4 0.6 0.8 11234567 Figure 13.8. Haar Wavelet Expansion. Remark: To the novice, there may appear to be many more wavelets than trigono- metric functions. But this is just another illusion of the ma gic show of infinite dimensional space. The point is that they both form a countably infinite se t 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 ma ny more Haar wavelets. This handicap makes the Haar system of less use in practical situa tions, and served to motivate the search for a more sophisticated choice of wavelet basis. Just as the discrete Fourier representation arises from a sa mpled version of the full Fourier series, so there is a discrete wavelet transformati on for suitably sampled signals. To take full advantage of the wavelet basis, we sample the sig nalf(x) atn= 2requally spaced sample points xk=k/2r, fork= 0,...,n−1, on the interval [0 ,1]. As before, we can identify the sampled signal with the vector f=/parenleftbig f(x0),f(x1),...,f(xn−1)/parenrightbigT=/parenleftbig f0,f1,...,fn−1/parenrightbigT∈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= 2rwavelets ϕ(x) andwj,k(x) forj= 0,...,r−1 andk= 0,...,2j−1. Letw0∼ϕ(x) andwj,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).) Mor eover, orthogonality of the wavelets immediately implies that the wavelet vectors f orm an orthogonal basis of Rn withn= 2r. We can decompose our sample vector (13.53) as a linear combi nation of the sampled wavelets, f=/hatwidec0w0+r−1/summationdisplay j=02j−1/summationdisplay k=0/hatwidecj,kwj,k, (13.54) 12/11/12 699 c/circleco√yrt2012 Peter J. Olver where, by our usual orthogonality formulae, /hatwidec0=/angbracketleftf;w0/angbracketright /bardblw0/bardbl2=1 2r2r−1/summationdisplay i=0fi, /hatwidecj,k=/angbracketleftf;wj,k/angbracketright /bardblwj,k/bardbl2= 2j−r k+2r−j−1−1/summationdisplay i=kfi−k+2r−j−1/summationdisplay i=k+2r−j−1fi .(13.55) These are the basic formulae connecting the functions f(x), or, rather, its sample vector f, and its discrete wavelet transform consisting of the 2rcoefficients/hatwidec0,/hatwidecj,k. The recon- structed function /tildewidef(x) =/hatwidec0ϕ(x)+r−1/summationdisplay j=02j−1/summationdisplay k=0/hatwidecj,kwj,k(x) (13 .56) is constant on each subinterval of length 2−r, and has the same value /tildewidef(x) =/tildewidef(xi) =f(xi) =fi, xi= 2−ri≤x < xi+1= 2−r(i+1), as our signal at the left hand endpoint of the interval. In oth er 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 provid e 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 f or this is that the Haar wavelets are piecewise constant, and so even an affine functio ny=αx+βrequires many sample values, and hence a relativelyextensive collection of Haarwavelets, to beaccurately reproduced. In particular, compression and denoising algo rithms based on Haar wavelets are either insufficiently precise or hopelessly inefficient, a nd hence of minor practical value. For a long time it was thought that it was impossible to simult aneously 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 In- grid Daubechies produced the first examples of wavelet bases that realized all three basic criteria. Since then, wavelets have developed into a sophis ticated and burgeoning industry with major impact on modern technology. Significant applica tions include compression, storage and recognition of fingerprints in the FBI’s data bas e, and the JPEG2000 image format, which, unlike earlier Fourier-based JPEG standard s, incorporates wavelet tech- nology in its image compression and reconstruction algorit hms. In this section, we will present a brief outline of the basic ideas underlying Daubec hies’ remarkable construction. The recipe for any wavelet system involves two basic ingredi ents — a scaling function and a mother wavelet. The latter can be constructed from the s caling function by a prescription similar to that in (13.42), and therefore we fir st concentrate on the properties of the scaling function. The key requirement is that the scal ing function must solve a 12/11/12 700 c/circleco√yrt2012 Peter J. Olver -2 -1 1 2 3 4 -0.20.20.40.60.811.2 Figure 13.9. The Hat Function. dilation equation of the form ϕ(x) =p/summationdisplay k=0ckϕ(2x−k) =c0ϕ(2x)+c1ϕ(2x−1)+···+cpϕ(2x−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,...,cpare not arbitrary, since the properties of orthogonality and lo calization will impose certain rather stringent requirements. Example 13.10. The Haar or box scaling function (13.40) satisfies the dilati on equation (13.57) with c0=c1= 1, namely ϕ(x) =ϕ(2x)+ϕ(2x−1). (13.58) We recommend that you convince yourself of the validity of th is identity before continuing. Example 13.11. Another example of a scaling function is the hat function ϕ(x) =  x, 0≤x≤1, 2−x,1≤x≤2, 0, otherwise ,(13.59) graphed in Figure 13.9, whose variants play a starring role i n the finite element method, cf. (11.167). The hat function satisfies the dilation equati on ϕ(x) =1 2ϕ(2x)+ϕ(2x−1)+1 2ϕ(2x−2), (13.60) which is (13.57) with c0=1 2,c1= 1,c2=1 2. Again, the reader should be able to check this identity by hand. 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 equa tions remains much less well developed than that of differential equations or integral eq uations. Even to prove that (nonzero) solutions exist is a nontrivial analytical probl em. Since we already know two 12/11/12 701 c/circleco√yrt2012 Peter J. Olver explicit examples, let us defer the discussion of solution t echniques until we understand how the dilation equation can be used to construct a wavelet b asis. Given a solution to the dilation equation, we define the mother wavelet to be w(x) =p/summationdisplay k=0(−1)kcp−kϕ(2x−k) =cpϕ(2x)−cp−1ϕ(2x−1)+cp−2ϕ(2x−2)+··· ±c0ϕ(2x−p),(13.61) This formula directly generalizes the Haar wavelet relatio n (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 wavele t: wj,k(x) =w(2jx−k). (13.62) In the general framework, we do not necessarily restrict our attention to the interval [0 ,1] and sojandkcan, in principle, be arbitrary integers. Let us investigate what sort of conditions should be imposed on the dilation coeffi- cientsc0,...,cpin order that we obtain a viable wavelet basis by this constru ction. First, localization of the wavelets requires that the scaling func tion has bounded support, and soϕ(x)≡0 whenxlies outside some bounded interval [ a,b]. If we integrate both sides of (13.57), we find /integraldisplayb aϕ(x)dx=/integraldisplay∞ −∞ϕ(x)dx=p/summationdisplay k=0ck/integraldisplay∞ −∞ϕ(2x−k)dx. (13.63) Now using the change of variables y= 2x−k, withdx=1 2dy, we find /integraldisplay∞ −∞ϕ(2x−k)dx=1 2/integraldisplay∞ −∞ϕ(y)dy=1 2/integraldisplayb aϕ(x)dx, (13.64) where we revert to xas our (dummy) integration variable. We substitute this res ult back into (13.63). Assuming that/integraldisplayb aϕ(x)dx/negationslash= 0, we discover that the dilationcoefficients must satisfy c0+···+cp= 2. (13.65) Example 13.12. Once we impose the constraint (13.65), the very simplest ver sion of the dilation equation is ϕ(x) = 2ϕ(2x) (13 .66) wherec0= 2istheonly(nonzero)coefficient. Uptoconstantmultiple, theonly“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 . Other solutions, such as ϕ(x) = 1/x, are not localized, and thus not useful for constructing a wa velet basis. 12/11/12 702 c/circleco√yrt2012 Peter J. Olver The second condition we require is orthogonality of the wave lets. For simplicity, we only consider the standard L2inner product† /angbracketleftf;g/angbracketright=/integraldisplay∞ −∞f(x)g(x)dx. It turns out that the orthogonality of the complete wavelet s ystem is guaranteed once we know that the scaling function ϕ(x) is orthogonal to all its integer translates: /angbracketleftϕ(x);ϕ(x−m)/angbracketright=/integraldisplay∞ −∞ϕ(x)ϕ(x−m)dx= 0 for all m/negationslash= 0.(13.67) We first note the formula /angbracketleftϕ(2x−k);ϕ(2x−l)/angbracketright=/integraldisplay∞ −∞ϕ(2x−k)ϕ(2x−l)dx (13.68) =1 2/integraldisplay∞ −∞ϕ(x)ϕ(x+k−l)dx=1 2/angbracketleftϕ(x);ϕ(x+k−l)/angbracketright follows from the same change of variables y= 2x−kused in (13.64). Therefore, since ϕ satisfies the dilation equation (13.57), /angbracketleftϕ(x);ϕ(x−m)/angbracketright=/angbracketleftiggp/summationdisplay j=0cjϕ(2x−j);p/summationdisplay k=0ckϕ(2x−2m−k)/angbracketrightigg (13.69) =p/summationdisplay j,k=0cjck/angbracketleftϕ(2x−j);ϕ(2x−2m−k)/angbracketright=1 2p/summationdisplay j,k=0cjck/angbracketleftϕ(x);ϕ(x+j−2m−k)/angbracketright. If we require orthogonality (13.67) of all the integer trans lates ofϕ, then the left hand side of this identity will be 0 unless m= 0, while only the summands with j= 2m+kwill be nonzero on the right. Therefore, orthogonality requires th at /summationdisplay 0≤k≤p−2mc2m+kck=/braceleftbigg2, m= 0, 0, m/negationslash= 0.(13.70) The algebraic equations (13.65,70) for the dilation coeffici ents 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 0+c2 1= 2, and soc0=c1= 1 is the only solution, resulting in the Haar dilation equat ion (13.58). If we have three coefficients c0,c1,c2, then (13.65), (13.70) require c0+c1+c2= 2, c2 0+c2 1+c2 2= 2, c0c2= 0. †In all instances, the functions have bounded support, and so the inn er product integral can be reduced to an integral over a finite interval where both fandgare nonzero. 12/11/12 703 c/circleco√yrt2012 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 reformulati on of the Haar case; see Exercise . In particular, the hat function (13.59) does notgive rise to orthogonal wavelets. The remarkable fact, discovered by Daubechies, is that ther eisa nontrivial solution for four (and, indeed, any even number) of nonzero coefficient sc0,c1,c2,c3. The basic equations (13.65), (13.70) require c0+c1+c2+c3= 2, c2 0+c2 1+c2 2+c2 3= 2, c0c2+c1c3= 0.(13.71) The particular values c0=1+√ 3 4, 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ϕ(2x)+3+√ 3 4ϕ(2x−1)+3−√ 3 4ϕ(2x−2)+1−√ 3 4ϕ(2x−3).(13.73) Any nonzero solution of bounded support to this remarkable f unctional equation will give rise to a scaling function ϕ(x), a mother wavelet w(x) =1−√ 3 4ϕ(2x)−3−√ 3 4ϕ(2x−1)+3+√ 3 4ϕ(2x−2)−1+√ 3 4ϕ(2x−3),(13.74) and then, by compression and translation (13.62), the compl ete system of orthogonal wavelets wj,k(x). Before explaining how to solve the Daubechies dilation equa tion, let us complete the proof of orthogonality. It is easy to see that, by translatio n invariance, since ϕ(x) and ϕ(x−m) are orthogonal for any m/negationslash= 0, so are ϕ(x−k) andϕ(x−l) for any k/negationslash=l. Next we prove orthogonality of ϕ(x−m) andw(x): /angbracketleftw(x);ϕ(x−m)/angbracketright=/angbracketleftiggp/summationdisplay j=0(−1)j+1cjϕ(2x−1+j);p/summationdisplay k=0ckϕ(2x−2m−k)/angbracketrightigg =p/summationdisplay j,k=0(−1)j+1cjck/angbracketleftϕ(2x−1+j);ϕ(2x−2m−k)/angbracketright =1 2p/summationdisplay j,k=0(−1)j+1cjck/angbracketleftϕ(x);ϕ(x−1+j−2m−k)/angbracketright, using (13.68). By orthogonality (13.67) of the translates o fϕ, the only summands that are nonzero are when j= 2m+k+1; the resulting coefficient of /bardblϕ(x)/bardbl2is /summationdisplay k(−1)kc1−2m−kck= 0, where the sum is over all 0 ≤k≤psuch that 0 ≤1−2m−k≤p. Each term in the sum appears twice, with opposite signs, and hence the result is a lways zero — no matter what the coefficients c0,...,cpare! The proof of orthogonality of the translates w(x−m) of the mother wavelet, along with all her wavelet descendants w(2jx−k), relies on a similar argument, and the details are left as an exercise for the read er. 12/11/12 704 c/circleco√yrt2012 Peter J. Olver Figure 13.10. Approximating the Daubechies Wavelet. Solving the Dilation Equation Let us next discuss how to solve the dilation equation (13.57 ). The solution we are after doesnothaveanelementary formula, andwerequireasl ightlysophisticatedapproach to recover it. The key observation is that (13.57) has the for m of a fixed point equation ϕ=F[ϕ], not in ordinary Euclidean space, but in an infinite-dimensio nal function space. With luck, the fixed point (or, more correctly, fixed function) will be st able, 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 ϕn+1(x) =p/summationdisplay k=0ckϕn(2x−k), n = 0,1,2,... . (13.75) Before attempting to prove convergence of this iterative pr ocedure to the Daubechies scal- ing function, let us experimentally investigate what happe ns. A reasonable choice for the initial guess might be the Haar sc aling or box function ϕ0(x) =/braceleftbigg1,0< t≤1. 0,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 ca n now try to prove convergence of the iterative scheme. This turns out to be true; a fully rig orous proof relies on the Fourier transform, [ 52], but is a little too advanced for this text and will be omitte d. 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 . 12/11/12 705 c/circleco√yrt2012 Peter J. Olver Once we have established convergence, we are now able to veri fy that the scaling function and consequential system of wavelets form an ortho gonal system of functions. Proposition 13.14. All integer translates ϕ(x−k), fork∈Zof the Daubechies scaling function, and all wavelets wj,k(x) =w(2jx−k),j≥0, are mutually orthogonal functions withrespect tothe L2inner product. Moreover, /bardblϕ/bardbl2= 1, while/bardblwj,k/bardbl2= 2−j. Proof: Asnoted earlier, theorthogonalityof the entire wavelet s ystem willfollow once we know the orthogonality (13.67) of the scaling function an d 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 t his property. Details are relegated to Exercise . Q.E.D. In practical computations, the limiting procedure for cons tructing the scaling function is not so convenient, and an alternative means of computing i ts 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∈Zexceptϕ0(1) = 1. The iterative functional equation(13.75)willthen producethe valuesoftheiterate sϕ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 itera tive system v(n+1)=Av(n)(13.76) with coefficient matrix A=/parenleftigg 3+√ 3 41+√ 3 4 1−√ 3 43−√ 3 4/parenrightigg and where v(n)=/parenleftbigg ϕn(1) ϕn(2)/parenrightbigg . Referring back to Chapter 10, the solution to such an iterati ve system is specified by the eigenvalues and eigenvectors of the coefficient matrix, w hich are λ1= 1,v1=/parenleftigg 1+√ 3 4 1−√ 3 4/parenrightigg , λ2=1 2,v2=/parenleftbigg −1 1/parenrightbigg . We write the initial condition as a linear combination of the eigenvectors v(0)=/parenleftbigg ϕ0(1) ϕ0(2)/parenrightbigg =/parenleftbigg 1 0/parenrightbigg = 2v1−1−√ 3 2v2. The solution is v(n)=Anv(0)= 2Anv1−1−√ 3 2Anv2= 2v1−1 2n1−√ 3 2v2. The limiting vector/parenleftigg ϕ(1) ϕ(2)/parenrightigg = lim n→∞v(n)= 2v1=/parenleftigg 1+√ 3 2 1−√ 3 2/parenrightigg 12/11/12 706 c/circleco√yrt2012 Peter J. Olver -0.5 0.5 11.5 22.5 33.5 -1.5-1-0.50.511.52 -0.5 0.5 11.5 22.5 33.5 -1.5-1-0.50.511.52 Figure 13.11. The Daubechies Scaling Function and Mother Wavelet. gives the desired values of the scaling function: ϕ(1) =1−√ 3 2= 1.366025... , ϕ (2) =−1+√ 3 2=−.366025... , ϕ(m) = 0,for all m/negationslash= 1,2.(13.77) With this in hand, the Daubechies dilation equation (13.73) then prescribes the func- tion values ϕ/parenleftbig1 2m/parenrightbig at all half integers, because when x=1 2mthen 2x−k=m−kis an integer. Once we know its values at the half integers, we can r e-use equation(13.73) to give its values at quarter integers1 4m. Continuing onwards, we determine the values of ϕ(x) at alldyadic points , meaning rational numbers of the form x=m/2jform,j∈Z. Continuity will then prescribe its value at any other x∈Rsincexcan be written as the limit of dyadic numbers xn— namely those obtained by truncating its binary (base 2) exp ansion at the nthdigit beyond the decimal (or, rather “binary”) point. But, i n 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 memor y. 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 Dau bechies scaling function in Figure 13.11. It is continuous, but non-differen tiable function — and its graph has a very jagged, fractal-like appearance when viewed at cl ose range. The Daubechies scaling function is, in fact, a close relative of the famous e xample of a continuous, nowhere differentiable function originally due to Weierstrass, [ 128,158], whose construction also relies on a similar scaling argument. With the values of the Daubechies scaling function on a suffici ently dense set of dyadic pointsinhand, theconsequential valuesofthemother wavel et aregivenby formula (13.74). Note that supp ϕ= suppw= [0,3]. The daughter wavelets are then found by the usual compression and translation procedure (13.62). The Daubechies wavelet expansion of a function whose suppor t is contained in†[0,1] †For functions with larger support, one should include additional terms in the expansion corresponding to further translates of the wavelets so as to cover t he entire support of the function. Alternatively, one can translate and rescale xto fit the function’s support inside [0 ,1]. 12/11/12 707 c/circleco√yrt2012 Peter J. Olver 0.2 0.4 0.6 0.8 12468 0.2 0.4 0.6 0.8 12468 0.2 0.4 0.6 0.8 12468 0.2 0.4 0.6 0.8 12468 0.2 0.4 0.6 0.8 12468 0.2 0.4 0.6 0.8 12468 Figure 13.12. Daubechies Wavelet Expansion. is then given by f(x)∼c0ϕ(x)+∞/summationdisplay j=02j−1/summationdisplay k=−2cj,kwj,k(x). (13.78) The inner summation begins at k=−2 so asto include allthe wavelet offspring wj,kwhose supporthasanontrivialintersectionwiththeinterval[0 ,1]. Thewaveletcoefficients c0,cj,k are computed by the usual orthogonality formula c0=/angbracketleftf;ϕ/angbracketright=/integraldisplay3 0f(x)ϕ(x)dx, cj,k=/angbracketleftf;wj,k/angbracketright= 2j/integraldisplay2−j(k+3) 2−jkf(x)wj,k(x)dx=/integraldisplay3 0f/parenleftbig 2−j(x+k)/parenrightbig w(x)dx,(13.79) where we agree that f(x) = 0 whenever x <0 orx >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 complet eness 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 ar e left as exercises for the reader to implement. Example 13.15. In Figure 13.12, we plot the Daubechies wavelet expansions o f the same signal for Example 13.9. The first plot is the origina l signal, and the following show the partial sums of (13.78) over j= 0,...,rwithr= 2,3,4,5,6. Unlike the Haar expansion, the Daubechies wavelets do exhibit the nonunifo rm Gibbs phenomenon at the interior discontinuity as well as the endpoints, since the f unction is set to 0 outside the interval [0 ,1]. Indeed, the Daubechies wavelets are continuous, and so c annot converge uniformly to a discontinuous function. 13.3. The Fourier Transform. Fourier series and their ilk are designed to solve boundary v alue problems on bounded intervals. The extension of Fourier methods to unbounded in tervals — the entire real line 12/11/12 708 c/circleco√yrt2012 Peter J. Olver — leads naturally to the Fourier transform, which is a powerf ul 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 p artial 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 ca se of Fourier series. Althoughtherigorousdetailsareratherexacting,theunde rlyingideaisnotsodifficult. 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 ofbasic trigonometricfunctions. One evident app roach is to construct its Fourier series on progressively larger a nd larger intervals, and then take thelimitastheintervals’ lengthbecomes infinite. Thelimi tingprocess convertstheFourier sums into integrals, and the resulting representation of a f unction is renamed the Fourier transform. Sincewearedealing withaninfiniteinterval, th ereareno longerany periodicity requirements on the function f(x). Moreover, the frequencies represented in the Fourier transform are no longer constrained by the length of the inte rval, and so we are effectively decomposing a quite general, non-periodic function into a c ontinuous superposition of trigonometric functions of all possible frequencies. Let us present the details of this construction in a more conc rete form. The computa- tions will be significantly simpler if we work with the comple x version of the Fourier series from the outset. Our starting point is the rescaled Fourier s eries (12.85) on a symmetric interval [ −ℓ,ℓ] of length 2 ℓ, which we rewrite in the adapted form f(x)∼∞/summationdisplay ν=−∞/radicalbiggπ 2/hatwidefℓ(kν) ℓeikνx. (13.80) The sum is over the discrete collection of frequencies kν=πν ℓ, ν = 0,±1,±2,... , (13.81) corresponding to those trigonometric functions that have p eriod 2ℓ. For reasons that will soon become apparent, the Fourier coefficients of fare now denoted as cν=/angbracketleftf;eikνx/angbracketright=1 2ℓ/integraldisplayℓ −ℓf(x)e−ikνxdx=/radicalbiggπ 2/hatwidefℓ(kν) ℓ, so that /hatwidefℓ(kν) =1√ 2π/integraldisplayℓ −ℓf(x)e−ikνxdx. (13.82) This reformulation of the basic Fourier series formula allo ws us to smoothly pass to the limit when the interval’s length ℓ→ ∞. On an interval of length 2 ℓ, the frequencies (13.81) required to represent a function in Fourier series form are equally distributed, with interf requency spacing ∆k=kν+1−kν=π ℓ. (13.83) 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 12/11/12 709 c/circleco√yrt2012 Peter J. Olver anticipate that allpossible frequencies will be represented. Indeed, letting kν=kbe arbitrary in (13.82), and sending ℓ→ ∞, results in the infinite integral /hatwidef(k) =1√ 2π/integraldisplay∞ −∞f(x)e−ikxdx (13.84) known as the Fourier transform of the function f(x). Iff(x) is a reasonably nice function, e.g., piecewise continuous and decaying to 0 reasonably qui ckly as|x| → ∞, its Fourier transform/hatwidef(k) is defined for all possible frequencies −∞< k <∞. This formula will sometimes conveniently be abbreviated as /hatwidef(k) =F[f(x)], (13.85) whereFis theFourier transform operator . To reconstruct the function from its Fourier transform, we e mploy a similar limiting procedure on the Fourier series (13.80), which we first rewri te in a more suggestive form: f(x)∼1√ 2π∞/summationdisplay ν=−∞/hatwidefℓ(kν)eikνx∆k. (13.86) For each fixed value of x, the right hand side has the form of a Riemann sum, [ 9,159], over the entire frequency space −∞< k <∞, for the function gℓ(k) =/hatwidefℓ(k)eikx. Under reasonable hypotheses, as ℓ→ ∞, the functions /hatwidefℓ(k)→/hatwidef(k) converge to the Fourier transform; moreover, the interfrequency spacing ∆ k→0 and so one expects the Riemann sums to converge to the integral f(x)∼1√ 2π/integraldisplay∞ −∞/hatwidef(k)eikxdk. (13.87) The result is the inverse Fourier transform , and serves to recover the original signal from its Fourier transform. In abbreviated form f(x) =F−1[/hatwidef(k)] (13 .88) is the inverse of the Fourier transform operator (13.85). In this manner, the Fourier series (13.86) becomes a Fourier integral that reconstructs the fu nctionf(x) as a (continuous) superposition of complex exponentials eikxofallpossible frequencies. The contribution of each such exponential is the value of the Fourier transfor m/hatwidef(k), evaluated at the exponential’s frequency. It is worth pointing out that both the Fourier transform (13. 85) and its inverse (13.88) define linear maps on function space. This means that the Four ier transform of the sum of two functions is the sum of their individual transforms, w hile multiplying a function by a constant multiplies its Fourier transform by the same fact or: F[f(x)+g(x)] =F[f(x)]+F[g(x)] =/hatwidef(k)+/hatwideg(k), F[cf(x)] =cF[f(x)] =c/hatwidef(k).(13.89) A similar statement hold for the inverse Fourier transform F−1. 12/11/12 710 c/circleco√yrt2012 Peter J. Olver -2 -1 1 2 -0.20.20.40.60.811.2 -3 -2 -1 1 2 3 -0.50.511.52 -2 -1 1 2 -0.20.20.40.60.811.2 -2 -1 1 2 -0.20.20.40.60.811.2 Figure 13.13. Fourier Transform of a Rectangular Pulse. Recapitulating, by lettingthe 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 a complete function /hatwidef(k) defined on all of frequency space k∈R. The reconstruction of f(x) from its Fourier transform/hatwidef(k) via (13.87) can be rigorously justified under suitable hypo theses. For example, if f(x) is piecewise C1on all of Rand decays reasonably rapidly, f(x)→0 as |x| → ∞, in order that its Fourier integral (13.84) converges, then it can be proved, [ 63], that the inverse Fourier integral (13.87) will converge to f(x) at all points of continuity, and to the midpoint1 2(f(x−)+f(x+)) at jump discontinuities — just like a Fourier series. In particular, its Fourier transform /hatwidef(k)→0 must also decay as |k| →0, implying that (as with Fourier series) the very high frequency modes make n egligible contributions to the reconstruction of the signal. A more precise, general re sult 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) =/braceleftbigg1,−a < x < a, 0,|x|> a.(13.90) of width 2 ais easily computed: /hatwidef(k) =1√ 2π/integraldisplaya −ae−ikxdx=eika−e−ika √ 2πik=/radicalbigg 2 πsinak k. (13.91) On the other hand, the reconstruction of the pulse via the inv erse transform (13.87) tells us that 1 π/integraldisplay∞ −∞eikxsinak kdk=f(x) =  1,−a < x < a, 1 2, x=±a, 0,|x|> a.(13.92) 12/11/12 711 c/circleco√yrt2012 Peter J. Olver Note the convergence to the middle of the jump discontinuiti es atx=±a. Splitting this complex integral into its real and imaginary parts, we deduc e a pair of striking trigono- metric integral identities 1 π/integraldisplay∞ −∞coskxsinak kdk=  1,−a < x < a, 1 2, x=±a, 0,|x|> a,1 π/integraldisplay∞ −∞sinkxsinak kdk= 0. (13.93) Just as many Fourier series yield nontrivial summation form ulae, the reconstruction of a function from its Fourier transform often leads to nontr ivial integral identities. One cannotcompute the integral (13.92) by the Fundamental Theorem of C alculus, since there is no elementary function†whose derivative equals the integrand. Moreover, it is not e ven 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,171], fails to apply. Thus, the convergence of the integral is marginal at best: the trigonometric oscil lations somehow overcome the slow rate of decay of 1 /kand thereby induce the (conditional) convergence of the int egral. In Figure 13.13 we display the box function with a= 1, its Fourier transform, along with a reconstruction obtained by numerically integrating (13. 93). Since we are dealing with an infinite integral, we must break off the numerical integrat or 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 what we observed in Fourier series. Example 13.17. Consider an exponentially decaying right-handed pulse‡ fr(x) =/braceleftbigge−ax, x > 0, 0, x < 0,(13.94) wherea >0. We compute its Fourier transform directly from the definit ion: /hatwidef(k) =1√ 2π/integraldisplay∞ 0e−axe−ikxdx=−1√ 2πe−(a+ik)x a+ ik/vextendsingle/vextendsingle/vextendsingle/vextendsingle∞ x=0=1√ 2π(a+ ik). As in the preceding example, the inverse Fourier transform p roduces a nontrivial complex †One can use Euler’s formula (12.51) to reduce the integral to a version of theexponential in- tegral/integraldisplay (eαk/k)dk, but it can be proved, [ 31], that there is no formula for in terms of elementary functions. ‡Note that we can’t Fourier transform the entire exponential function e−axbecause it does not go to zero at both ±∞, which is required for the integral (13.84) to converge. 12/11/12 712 c/circleco√yrt2012 Peter J. Olver -3 -2 -1 1 2 3 -0.20.20.40.60.81 Right Pulse fr(x)-3 -2 -1 1 2 3 -0.20.20.40.60.81 Left Pulse fl(x) -3 -2 -1 1 2 3 -0.20.20.40.60.81 Even Pulse fe(x)-3 -2 -1 1 2 3 -1-0.50.51 Odd Pulse fo(x) Figure 13.14. Exponential Pulses. integral identity: 1 2π/integraldisplay∞ −∞eikx a+ ikdk=  e−ax, x > 0, 1 2, x = 0, 0, x < 0.(13.95) Similarly, a pulse that decays to the left, fl(x) =/braceleftbiggeax, x < 0, 0, x > 0,(13.96) wherea >0 is still positive, has Fourier transform /hatwidef(k) =1√ 2π(a−ik). (13.97) This also follows from the general fact that the Fourier tran sform of f(−x) is/hatwidef(−k); see Exercise . The even exponentially decaying pulse fe(x) =e−a|x|(13.98) is merely the sum of left and right pulses: fe=fr+fl. Thus, by linearity, /hatwidefe(k) =1√ 2π(a+ ik)+1√ 2π(a−ik)=/radicalbigg 2 πa k2+a2, (13.99) The result is real and even because fe(x) is a real even function; see Exercise . The inverse Fourier transform (13.87) produces another nontri vial integral identity: e−a|x|=1 π/integraldisplay∞ −∞aeikx k2+a2dk=a π/integraldisplay∞ −∞coskx k2+a2dk. (13.100) 12/11/12 713 c/circleco√yrt2012 Peter J. Olver (The imaginary part of the integral vanishes because the int egrand is odd.) On the other hand, the odd exponentially decaying pulse, fo(x) = (sign x)e−a|x|=/braceleftbigge−ax, x > 0, −eax, x < 0,(13.101) is the difference of the right and left pulses, fo=fr−fl, and has purely imaginary and odd Fourier transform /hatwidefo(k) =1√ 2π(a+ ik)−1√ 2π(a−ik)=−i/radicalbigg 2 πk k2+a2. (13.102) The inverse transform is (signx)e−a|x|=−i π/integraldisplay∞ −∞keikx k2+a2dk=1 π/integraldisplay∞ −∞ksinkx k2+a2dk. (13.103) As a final example, consider the rational function f(x) =1 x2+a2,where a >0. (13.104) Its Fourier transform requires integrating /hatwidef(k) =1√ 2π/integraldisplay∞ −∞e−ikx x2+a2dx. (13.105) The indefinite integral (anti-derivative) does not appear i n basic integration tables, and, in fact, cannot be done in terms of elementary functions. How ever, we have just managed to evaluate this particular integral! Look at (13.100). If w e change xtokandkto−x, then we exactly recover the integral (13.105) up to a factor o fa/radicalbig 2/π. We conclude that the Fourier transform of (13.104) is /hatwidef(k) =/radicalbiggπ 2e−a|k| a. (13.106) This last example is indicative of an important general fact . The reader has no doubt already noted the remarkable similarity between the Fourie r transform (13.84) and its inverse (13.87). Indeed, the only difference is that the form er 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/hatwidef(k), then the Fourier transform of /hatwidef(x)isf(−k). The Symmetry Principle allows us to reduce the tabulation of Fourier transforms by half. For instance, referring back to Example 13.16, we dedu ce that the Fourier transform of the function f(x) =/radicalbigg 2 πsinax x 12/11/12 714 c/circleco√yrt2012 Peter J. Olver is /hatwidef(k) =σ(−k+a)−σ(−k−a) =σ(k+a)−σ(k−a) =  1,−a < k < a, 1 2, k=±a, 0,|k|> a.(13.107) Note that, by linearity, we can divide both f(x) and/hatwidef(k) by/radicalbig 2/πto deduce the Fourier transform ofsinax x. Warning : Some authors leave out the√ 2πfactor in the definition (13.84) of the Fourier transform /hatwidef(k). This alternative convention does have a slight advantage of elimi- nating many√ 2πfactors in the Fourier transforms of elementary functions. However, this necessitates an extra such factor in the reconstruction for mula (13.87), which is achieved by replacing√ 2πby 2π. A significant disadvantage is that the resulting formulae f or the Fourier transform and its inverse are not as close, and so the Symmetry Principle of Theo- rem 13.18 requires some modification. (On the other hand, con volution — to be discussed below — is a little easier without the extra factor.) Yet anot her, more recent convention can be found in Exercise . When consulting any particular reference book, the reader alwaysneeds to check which convention is being used. All of the functions in Example 13.17 required a >0 for the Fourier integrals to con- verge. The functions that emerge inthe limitas agoes to 0 are offundamental importance. Let us start with the odd exponential pulse (13.101). When a→0, the function fo(x) converges to the sign function f(x) = signx=σ(x)−σ(−x) =/braceleftbigg+1, x > 0, −1, x < 0.(13.108) Taking the limit of the Fourier transform (13.102) leads to /hatwidef(k) =−i/radicalbigg 2 π1 k. (13.109) The nonintegrable singularity of /hatwidef(k) atk= 0 is indicative of the fact that the sign func- tion does notdecay as |x| → ∞. In this case, neither the Fourier transform integral nor it s inverse are well-defined as standard (Riemann, or even Lebes gue, [159]) integrals. Never- theless, 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.110) The limit of the Fourier transform (13.99) is lim a→0/radicalbigg 2 π2a k2+a2=/braceleftbigg0, k/negationslash= 0, ∞, k= 0.(13.111) 12/11/12 715 c/circleco√yrt2012 Peter J. Olver This limiting behavior should remind the reader of our const ruction (11.31) of the delta function as the limit of the functions δ(x) = lim n→∞n π(1+n2x2)= lim a→0a π(a2+x2), Comparing with (13.111), we conclude that the Fourier trans form of the constant function (13.110) is a multiple of the delta function in the frequency variable: /hatwidef(k) =√ 2π δ(k). (13.112) The direct transform integral δ(k) =1 2π/integraldisplay∞ −∞e−ikxdx is, strictly speaking, not defined because the infinite integ rals of the oscillatory sine and cosine functions don’t converge! However, this identity ca n be validly interpreted within the framework of weak convergence and generalized function s. On the other hand, the inverse transform formula (13.87) yields /integraldisplay∞ −∞δ(k)eikxdk=eik0= 1, which is in accord with the basic definition (11.37) of the del ta 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 /hatwideδ(k) =1√ 2π/integraldisplay∞ −∞δ(x)e−ikxdx=e−ik0 √ 2π≡1√ 2π. (13.113) a result that also follows from the symmerty principle of The orem 13.18. To determine the Fourier transform of a delta spike δy(x) =δ(x−y) concentrated at position x=y, we compute /hatwideδy(k) =1√ 2π/integraldisplay∞ −∞δ(x−y)e−ikxdx=e−iky √ 2π. (13.114) Theresult isapure exponential infrequency space. Applyin g theinverse Fouriertransform (13.87) leads, formally, to the remarkable identity δy(x) =δ(x−y) =1 2π/integraldisplay∞ −∞e−ik(x−y)dk=1 2π/angbracketlefteiky;eikx/angbracketright, (13.115) where/angbracketleft·;·/angbracketrightdenotes the usual L2Hermitian inner product on R. Since the delta function vanishes for x/negationslash=y, this identity is telling us that complex exponentials of di ffering frequen- cies are mutually orthogonal. However, this statement must be taken with a grain of salt, since the integral does not converge in any standard sense (e ither Riemann or Lebesgue). 12/11/12 716 c/circleco√yrt2012 Peter J. Olver But it is possible to make sense of this identity within the la nguage of generalized func- tions. Indeed, multiplying both sides by f(x), and then integrating with respect to x, we find f(y) =1 2π/integraldisplay∞ −∞/integraldisplay∞ −∞f(x)e−ik(x−y)dxdk. (13.116) Thisisa perfectly valid formula, being a restatement (or, rather, combination) of the basic formulae (13.84,87) connecting the direct and inverse Four ier transforms of the function f(x). Conversely, the Symmetry Principle tells us that the Fourie r transform of a pure exponential eilxwill be a shifted delta spike√ 2π δ(k−l), concentrated in frequency space. Both results are particular cases of the general Shif t Theorem, whose proof is left as an exercise for the reader. Theorem 13.19. Iff(x)has Fourier transform /hatwidef(k), then the Fourier transform of the shifted function f(x−y)ise−iky/hatwidef(k). Similarly, the transform of the product function eilxf(x)forlreal is the shifted transform /hatwidef(k−l). Since the Fourier transform uniquely associates a function /hatwidef(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 (o r, even better, computer algebra systems such as Mathematica orMaple) that recognize a wide variety of transforms of basic functions of importance in applications. The accom panying table lists some of the most important examples of functions and their Fourier t ransforms. Note that, by applying the Symmetry Principle of Theorem 13.18, each tabu lar entry can be used to deduce two different Fourier transforms. A more extensive co llection of Fourier transforms can be found in [ 142]. Derivatives and Integrals One of the most remarkable and important properties of the Fo urier transform is that it converts calculus into algebra! More specifically, the tw o basic operations in calculus — differentiation and integration of functions — are realized as algebraic operations on their Fourier transforms. (The downside is that algebraic operat ions become more complicated in the transform domain.) Let us begin with derivatives. If we differentiate†the basic inverse Fourier transform formula f(x)∼1√ 2π/integraldisplay∞ −∞/hatwidef(k)eikxdk. with respect to x, we obtain f′(x)∼1√ 2π/integraldisplay∞ −∞ik/hatwidef(k)eikxdk. (13.117) †We are assuming the integrand is sufficiently nice so that we can bring the derivative under the integral sign; see [ 63,199] for a fully rigorous derivation. 12/11/12 717 c/circleco√yrt2012 Peter J. Olver Short Table of Fourier Transforms f(x) /hatwidef(k) 1√ 2π δ(k) δ(x)1√ 2π σ(x)/radicalbiggπ 2δ(k)−i√ 2π k signx −i/radicalbigg 2 π1 k σ(x+a)−σ(x−a)/radicalbigg 2 πsinak k e−axσ(x)1√ 2π(a+ ik) eax(1−σ(x))1√ 2π(a−ik) e−a|x|/radicalbigg 2 πa k2+a2 e−ax2 e−k2/4a √ 2a tan−1xπ3/2 √ 2δ(k)−i/radicalbiggπ 2e−|k| k f(cx+d)eikd/c |c|/hatwidef/parenleftbiggk c/parenrightbigg f(x) /hatwidef(−k) /hatwidef(x) f(−k) f′(x) ik/hatwidef(k) xf(x) i/hatwidef′(k) f∗g(x)√ 2π/hatwidef(k)/hatwideg(k) Note: The parameters a,c,dare real, with a >0 andc/negationslash= 0. 12/11/12 718 c/circleco√yrt2012 Peter J. Olver Theresultingintegralisitselfintheform ofaninverseFou riertransform, namelyof i k/hatwidef(k) which immediately implies the following key result. Proposition 13.20. The Fourier transform of the derivative f′(x)of a function is obtained by multiplication of its Fourier transform by ik: F[f′(x)] = ik/hatwidef(k). (13.118) Similarly, the Fourier transform of xf(x)is obtained by differentiating the Fourier trans- form off(x): F[xf(x)] = id/hatwidef dk. (13.119) The second statement follows from the first by use of the Symme try 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|: f′ e(x) =−a(signx)e−a|x|=−afo(x). Proposition 13.20 says that their Fourier transforms are re lated by ik/hatwidefe(k) = i/radicalbigg 2 πka k2+a2=−a/hatwidefo(k), as previously noted in (13.99,102). On the other hand, the od d exponential pulse has a jump discontinuity of magnitude 2 at x= 0, and so its derivative contains a delta function: f′ o(x) =−ae−a|x|+2δ(x) =−afe(x)+2δ(x). This is reflected in the relation between their Fourier trans forms. If we multiply (13.102) by ikwe obtain ik/hatwidefo(k) =/radicalbigg 2 πk2 k2+a2=/radicalbigg 2 π−/radicalbigg 2 πa2 k2+a2= 2/hatwideδ(k)−a/hatwidefe(k). The Fourier transform, just like Fourier series, is complet ely compatible with the calculus of generalized functions. Higher order derivatives are handled by iterating the first o rder formula (13.118). Corollary 13.22. The Fourier transform of f(n)(x)is(ik)n/hatwidef(k). This result has an important consequence: the smoothness of f(x) is manifested in the rate of decay of its Fourier transform /hatwidef(k). We already noted that the Fourier transform of a (nice) function must decay to zero at large fre quencies:/hatwidef(k)→0 as |k| → ∞. If the nthderivative f(n)(x) is also a reasonable function, then its Fourier transform/hatwidestf(n)(k) = (ik)n/hatwidef(k) must go to zero as |k| → ∞. This requires that /hatwidef(k) go to zero more rapidly than |k|−n. Thus, the smoother f(x), the more rapid the decay of its 12/11/12 719 c/circleco√yrt2012 Peter J. Olver Fourier transform. As a general rule of thumb, local feature s off(x), such as smoothness, are manifested by global features of /hatwidef(k), such as decay for large |k|. The Symmetry Principle implies that reverse is also true: global feature s off(x) correspond to local features of/hatwidef(k). This local-global duality is one of the major themes of Fou rier theory. Integration is the inverse operation to differentiation, an d so should correspond to division by 2 πikin frequency space. As with Fourier series, this is not compl etely correct; thereisanextraconstantinvolved,andthiscontributesan extradeltafunctioninfrequency space. Proposition 13.23. Iff(x)has Fourier transform /hatwidef(k), then the Fourier transform of its integral g(x) =/integraldisplayx −∞f(y)dyis /hatwideg(k) =−i k/hatwidef(k)+π/hatwidef(0)δ(k). (13.120) Proof: First notice that lim x→−∞g(x) = 0, lim x→+∞g(x) =/integraldisplay∞ −∞f(x)dx=√ 2π/hatwidef(0). Therefore, by subtracting a suitable multiple of the step fu nction from the integral, the resulting function h(x) =g(x)−√ 2π/hatwidef(0)σ(x) decays to 0 at both ±∞. Consulting our table of Fourier transforms, we find /hatwideh(k) =/hatwideg(k)−π/hatwidef(0)δ(k)+i k/hatwidef(0). (13.121) On the other hand, h′(x) =f(x)−√ 2π/hatwidef(0)δ(x). Sinceh(x)→0 as|x| → ∞, we can apply our differentiation rule (13.118), and conclud e that ik/hatwideh(k) =/hatwidef(k)−/hatwidef(0). (13.122) Combining (13.121) and (13.122) establishes the desired fo rmula (13.120). Q.E.D. Example 13.24. The Fourier transform of the inverse tangent function f(x) = tan−1x=/integraldisplayx 0dy 1+y2=/integraldisplayx −∞dy 1+y2−π 2 can be computed by combining Proposition 13.23 with (13.106 ): /hatwidef(k) =−i/radicalbiggπ 2e−|k| k+π3/2 √ 2δ(k). 12/11/12 720 c/circleco√yrt2012 Peter J. Olver Applications to Differential Equations The fact that the Fourier transform coverts differentiation in the physical domain into multiplication in the frequency domain is one of its most com pelling features. A partic- ularly important consequence is that it effectively transfo rms differential equations into algebraic equations, and thereby opens the door to their sol ution by elementary algebra! One begins by applying the Fourier transform to both sides of the differential equation under consideration. Solving the resulting algebraic equa tion will produce a formula for the Fourier transform of the desired solution, which can the n be immediately reconstructed via the inverse Fourier transform. The Fourier transform is particularly well adapted to bound ary 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 functi ons). In quantum mechanics, [124,130], these solutions are known as the bound states of the system, and correspond to subatomic particles that are trapped or localized in a regio n of space by some sort of force field. For example, the electrons in an atom are bound states l ocalized by the electrostatic attraction of the nucleus. As a specific example, consider the boundary value problem −d2u dx2+ω2u=h(x), −∞< x <∞, (13.123) whereω >0 is a positive constant. In lieu of boundary conditions, we r equire 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 k2/hatwideu(k)+ω2/hatwideu(k) =/hatwideh(k) relatingtheFouriertransformsof uandh. Unlikethedifferentialequation, thetransformed equation can be immediately solved for /hatwideu(k) =/hatwideh(k) k2+ω2. (13.124) Therefore, we can reconstruct the solution by applying the i nverse Fourier transform for- mula (13.87): u(x) =1√ 2π/integraldisplay∞ −∞/hatwideh(k)eikx k2+ω2dk. (13.125) For example, if the forcing function is an even exponential p ulse, h(x) =e−|x|with/hatwideh(k) =/radicalbigg 2 π1 k2+1, then (13.125) writes the solution as a Fourier integral: u(x) =1 π/integraldisplay∞ −∞eikx (k2+ω2)(k2+1)dk=1 π/integraldisplay∞ −∞coskx (k2+ω2)(k2+1)dk. 12/11/12 721 c/circleco√yrt2012 Peter J. Olver noting that the imaginary part of the complex integral vanis hes because the integrand is an odd function. (Indeed, if the forcing function is real, th e solution must also be real.) The Fourier integral can be explicitly evaluated by using pa rtial fractions to rewrite /hatwideu(k) =/radicalbigg 2 π1 (k2+ω2)(k2+1)=/radicalbigg 2 π1 ω2−1/parenleftbigg1 k2+1−1 k2+ω2/parenrightbigg , ω2/negationslash= 1, Thus, according to our Fourier Transform Table, the solutio n to this boundary value problem is u(x) =e−|x|−1 ωe−ω|x| ω2−1when ω2/negationslash= 1. (13.126) The reader may wish to verify that this function is indeed a so lution, meaning that it is twice continuously differentiable (which is not so immediat ely apparent from the formula), decays to 0 as |x| → ∞, and satisfies the differential equation everywhere. The “re sonant” caseω2= 1 is left as an exercise. Remark: The method of partial fractions that you learned in first yea r calculus is often an effective tool for evaluating (inverse) Fourier tra nsforms 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.124), its Fourier transform with respect to xis /hatwideG(k,y) =1√ 2πe−iky k2+ω2. which is the product of an exponential factor e−iky, representing the Fourier transform ofδy(x), times a multiple of the Fourier transform of the even expon ential pulse e−ω|x|. We apply Theorem 13.19, and conclude that the Green’s functi on for this boundary value problem is an exponential pulse centered at y, namely G(x,y) =1 2ωe−ω|x−y|. Observe that, as with other self-adjoint boundary value pro blems, the Green’s function is symmetric under interchange of xandy. As a function of x, it satisfies the homogeneous differential equation −u′′+ω2u= 0, except at the point x=ywhen its derivative has a jump discontinuity of unit magnitude. It also decays as |x| → ∞, as required by the boundary conditions. The Green’s function superpositi on principle tells us that the solution to the inhomogeneous boundary value problem (13.1 23) under a general forcing can be represented in the integral form u(x) =/integraldisplay∞ −∞G(x,y)h(y)dy=1 2ω/integraldisplay∞ −∞e−ω|x−y|h(y)dy. (13.127) The reader may enjoy recovering the particular exponential solution (13.126) from this integral formula. 12/11/12 722 c/circleco√yrt2012 Peter J. Olver Convolution The final Green’s function formula (13.127) 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. Theconvolution of scalar functions f(x) andg(x) is the scalar function h=f∗gdefined by the formula h(x) =f(x)∗g(x) =/integraldisplay∞ −∞f(x−y)g(y)dy. (13.128) We record the basic properties of the convolution product, l eaving their verification as exercises for the reader. All of these assume that the impl ied convolution integrals converge. (a)Symmetry : f∗g=g∗f, (b)Bilinearity :/braceleftbiggf∗(ag+bh) =a(f∗g)+b(f∗h), (af+bg)∗h=a(f∗h)+b(g∗h),a,b∈C, (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 nota unit for the convolution product; indeed, f∗1 =/integraldisplay∞ −∞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) =/integraldisplay∞ −∞f(x−y)δ(y)dy=f(x), which follows from the basic property (11.37) of the delta fu nction. In particular, our solution (13.127) has the form of a convol ution product between an even exponential pulse g(x) =1 2ωe−ω|x|and the forcing function: u(x) =g(x)∗h(x). On the other hand, its Fourier transform (13.124) is, up to a f actor, the ordinary multi- plicative product /hatwideu(k) =√ 2π/hatwideg(k)/hatwideh(k) of the Fourier transforms of gandh. In fact, this is a general property of the Fourier trans- form: convolution in the physical domain corresponds to mul tiplication in the frequency domain, and conversely. 12/11/12 723 c/circleco√yrt2012 Peter J. Olver Theorem 13.26. The Fourier transform of the convolution h(x) =f(x)∗g(x)of two functions is a multiple of the product of their Fourier tr ansforms : /hatwideh(k) =√ 2π/hatwidef(k)/hatwideg(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 : /hatwideh(k) =1√ 2π/hatwidef(k)∗/hatwideg(k). Proof: Combining the definition of the Fourier transform with the c onvolution for- mula (13.128), we find /hatwideh(k) =1√ 2π/integraldisplay∞ −∞h(x)e−ikxdx=1√ 2π/integraldisplay∞ −∞/integraldisplay∞ −∞f(x−y)g(y)e−ikxdxdy, where we are assuming that the integrands are sufficiently nic e to allow us to interchange the order of integration, [ 9]. Applying the change of variables z=x−yin the inner integral produces /hatwideh(k) =1√ 2π/integraldisplay∞ −∞/integraldisplay∞ −∞f(z)g(y)e−ik(y+z)dydz =√ 2π/parenleftbigg1√ 2π/integraldisplay∞ −∞f(z)e−ikzdz/parenrightbigg/parenleftbigg1√ 2π/integraldisplay∞ −∞g(y)e−ikydy/parenrightbigg =√ 2π/hatwidef(k)/hatwideg(k). The second statement can be proved similarly, or by simply no ting that it follows directly from the Symmetry Principle of Theorem 13.18. Q.E.D. Example 13.27. We already know, (13.107), that the Fourier transform of f(x) =sinx x is the box function /hatwidef(k) =/radicalbiggπ 2/bracketleftbig σ(k+1)−σ(k−1)/bracketrightbig =  /radicalbiggπ 2,−1< k <1, 0,|k|>1, We also know that the Fourier transform of g(x) =1 xis/hatwideg(k) =−i/radicalbiggπ 2signk. Therefore, the Fourier transform of their product h(x) =f(x)g(x) =sinx x2 12/11/12 724 c/circleco√yrt2012 Peter J. Olver -3 -2 -1 1 2 3 -1-0.50.51 Figure 13.15. The Fourier transform ofsinx x2. can be obtained by convolution: /hatwideh(k) =1√ 2π/hatwidef∗/hatwideg(k) =1√ 2π/integraldisplay∞ −∞/hatwidef(l)/hatwideg(k−l)dl =−i/radicalbiggπ 8/integraldisplay1 −1sign(k−l)dl=  i/radicalbiggπ 2k <−1, −i/radicalbiggπ 2k,−1< k <1, −i/radicalbiggπ 2k >1. A graph of/hatwideh(k) appears in Figure 13.15. The Fourier Transform on Hilbert Space While we do not have the space to embark on a fully rigorous tre atment of the theory underlying the Fourier transform, it is worth outlining a fe w of the most basic features. We have already noted that the Fourier transform, when define d, is a linear map, taking functions f(x)onphysical spacetofunctions /hatwidef(k)onfrequency space. Acriticalquestionis 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.84) is required to converge, and thisplacesrestrictions onthefunction and i tsasymptoticsat largedistances. It turns out the proper setting for the rigorous theory is the Hilbert space of complex- valued square-integrable functions — the same infinite-dim ensional 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 ent ire real line. Thus, the Hilbert space L2= L2(R) is the infinite-dimensional vector space †We leave aside the more advanced issues involving generalized funct ions in this subsection. 12/11/12 725 c/circleco√yrt2012 Peter J. Olver consisting of all complex-valued functions f(x) which are defined for all x∈Rand have finite L2norm: /bardblf/bardbl2=/integraldisplay∞ −∞|f(x)|2dx <∞. (13.129) For example, any piecewise continuous function that satisfi es the decay criterion |f(x)| ≤M |x|1/2+δ,for all sufficiently large |x| ≫0, (13.130) 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 an d identification of its elements is quite subtle. The inner product on the Hilbert space L2is prescribed in the usual manner, /angbracketleftf;g/angbracketright=/integraldisplay∞ −∞f(x)g(x)dx, (13.131) so that/bardblf/bardbl2=/angbracketleftf;f/angbracketright. The Cauchy–Schwarz inequality |/angbracketleftf;g/angbracketright| ≤ /bardblf/bardbl/bardblg/bardbl (13.132) ensures that the inner product integral is finite whenever f,g∈L2. Let us state the fundamental theorem governing the effect of t he Fourier transform on functions in Hilbert space. It can be regarded as a direct ana log of the pointwise conver- gence Theorem 12.7 for Fourier series. A fully rigorous proo f can be found in [ 179,199]. Theorem 13.28. Iff(x)∈L2issquare-integrable, then itsFouriertransform /hatwidef(k)∈ L2is a well-defined, square-integrable function of the freque ncy variable k. Iff(x)is continuously differentiable at a point x, then its inverse Fourier transform (13.87)equals its value f(x). More generally, if the right and left hand limits f(x−),f′(x−), andf(x+), f′(x+)exist, then the inverse Fourier transform integral converg es to the average value 1 2/bracketleftbig f(x−)+f(x+)/bracketrightbig . Thus, the Fourier transform /hatwidef=F[f] defines a linear transformation from L2func- tions ofxto L2functions 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/hatwidef(k) =F[f(x)]and/hatwideg(k) =F[g(x)], then/angbracketleftf;g/angbracketright=/angbracketleft/hatwidef;/hatwideg/angbracketright, i.e., /integraldisplay∞ −∞f(x)g(x)dx=/integraldisplay∞ −∞/hatwidef(k)/hatwideg(k)dk. (13.133) Proof: Let us sketch a formal proof that serves to motivate why this result is valid. We use the definition (13.84) of the Fourier transform to eval uate /integraldisplay∞ −∞/hatwidef(k)/hatwideg(k)dk=/integraldisplay∞ −∞/parenleftbigg1√ 2π/integraldisplay∞ −∞f(x)e−ikxdx/parenrightbigg/parenleftbigg1√ 2π/integraldisplay∞ −∞g(y)e+ikydy/parenrightbigg dk =/integraldisplay∞ −∞/integraldisplay∞ −∞f(x)g(y)/parenleftbigg1 2π/integraldisplay∞ −∞e−ik(x−y)dk/parenrightbigg dxdy. 12/11/12 726 c/circleco√yrt2012 Peter J. Olver Now according to (13.115), the inner kintegral can be replaced by a delta function δ(x−y), and hence /integraldisplay∞ −∞/hatwidef(k)/hatwideg(k)dk=/integraldisplay∞ −∞/integraldisplay∞ −∞f(x)g(y)δ(x−y)dxdy=/integraldisplay∞ −∞f(x)g(x)dx. This completes our “proof”; see [ 179,199] for a rigorous version. Q.E.D. In particular, orthogonal functions, /angbracketleftf;g/angbracketright= 0, will have orthogonal Fourier trans- forms,/angbracketleft/hatwidef;/hatwideg/angbracketright= 0. Choosing f=gin (13.133) results in the Plancherel formula /bardblf/bardbl2= /bardbl/hatwidef/bardbl2, or, explicitly,/integraldisplay∞ −∞|f(x|2dx=/integraldisplay∞ −∞|/hatwidef(k)|2dk. (13.134) We conclude that the Fourier transform map F:L2→L2defines a unitaryor norm- preserving linear transformation on Hilbert space, mappin g L2functions of the physical variable xto L2functions of the frequency variable k. Quantum Mechanics and the Uncertainty Principle In quantum mechanics, the wave functions or probability den sities of a quantum sys- tem are characterized as the elements of unit norm, /bardblϕ/bardbl= 1, belonging to the under- lying state space, which, in a one-dimensional model of a sin gle particle, is the Hilbert space L2= L2(R). All measurable physical quantities are represented by li near operators A:L2→L2acting on (a suitable dense subspace of) the state space. The se are obtained by the rather mysterious process of “quantizing” their clas sical counterparts, which are ordinary functions, [ 59,124,130]. In particular, the particle’s position, usually denoted byQ, is represented by the operation of multiplication by x, soQ[ϕ] =xϕ(x), whereas itsmomentum , denoted by Pfor historical reasons, is represented by the differentiati on operator P= id/dx, so that P[ϕ] = iϕ′(x). In its popularized form, Heisenberg’s Uncertainty Princip le is a by now familiar philo- sophical concept. It was first formulated in the 1920’s by the German physicist Werner Heisenberg, one of the founders of modern quantum mechanics , and states that, in a physi- cal system, certain quantities cannot be simultaneously me asured with complete accuracy. For instance, the more precisely one measures the position o f 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 u ncertainty couples energy and time. Experimental verification of the uncertainty princip le can be found even in fairly simple situations. Consider a light beam passing through a s mall hole. The position of the photons is constrained by the hole; the effect of their moment a is in the pattern of light diffused on a screen placed beyond the hole. The smaller the ho le, the more constrained the position, and the wider the image on the screen, meaning t he less certainty there is in the observed momentum. This is not the place to discuss the philosophical and experi mental consequences of Heisenberg’s principle. What we will show is that the Uncert ainty Principle is, in fact, a rigorous theorem concerning the Fourier transform! In quan tum theory, each of the paired quantities, e.g., positionandmomentum, areinterrelated bytheFouriertransform. Indeed, 12/11/12 727 c/circleco√yrt2012 Peter J. Olver Proposition 13.20 says that the Fourier transform of the diff erentiation operator represent- ing monmentum is a multiplication operator representing po sition on the frequency space and vice versa. This Fourier transform-based duality betwe en position and momentum, or multiplication and differentiation, lies at the heart of the Uncertainty Principle. In general, suppose Ais a linear operator on Hilbert space representing a physica l quantity. If the quantum system is in a state represented by a particular wave function ϕ, then the localization of the quantity Ais measured by the norm /bardblA[ϕ]/bardbl. The smaller this norm, the more accurate the measurement. Thus, /bardblQ[ϕ]/bardbl=/bardblxϕ(x)/bardblmeasures the localization of the position of the particle represented by ϕ; the smaller /bardblQ[ϕ]/bardbl, 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 Plan cherel’s formula (13.134), /bardblP[ϕ]/bardbl=/bardblϕ′(x)/bardbl=/bardblik/hatwideϕ(k)/bardbl=/bardblk/hatwideϕ(k)/bardbl, measures the localization in the momentum of the particle ne ar 0, 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 Uncer tainty Principle states that these two quantities cannot simultaneously be arbitrarily small. Theorem 13.30. Ifϕ(x)is a wave function, so /bardblϕ/bardbl= 1, then /bardblQ[ϕ]/bardbl /bardblP[ϕ]/bardbl ≥1 2. (13.135) Proof: The proof rests on the Cauchy–Schwarz inequality /vextendsingle/vextendsingle/vextendsingle/angbracketleftxϕ(x);ϕ′(x)/angbracketright/vextendsingle/vextendsingle/vextendsingle≤ /bardblxϕ(x)/bardbl /bardblϕ′(x)/bardbl=/bardblQ[ϕ]/bardbl /bardblP[ϕ]/bardbl. (13.136) On the other hand, writing out the inner product term /angbracketleftxϕ(x);ϕ′(x)/angbracketright=/integraldisplay∞ −∞xϕ(x)ϕ′(x)dx. Let us integrate by parts, using the fact that ϕ(x)ϕ′(x) =d dx/bracketleftbig1 2ϕ(x)2/bracketrightbig . Sinceϕ(x)→0 as|x| → ∞, the boundary terms vanish, and hence /angbracketleftxϕ(x);ϕ′(x)/angbracketright=/integraldisplay∞ −∞xϕ(x)ϕ′(x)dx=−/integraldisplay∞ −∞1 2ϕ(x)2dx=−1 2, since/bardblϕ/bardbl= 1. Substituting back into (13.136) completes the proof. Q.E.D. The inequality (13.135) quantifies the statement that the mo re 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 shou ld consult an introductory text on mathematical quantum mechanics, e.g., [ 124,130]. †To measure it at another position a, one replaces xbyx−a. 12/11/12 728 c/circleco√yrt2012 Peter J. Olver 13.4. The Laplace Transform. In engineering applications, the Fourier transform is ofte n overshadowed by a close relative. The Laplace transform plays an essential role in c ontrol theory, linear systems analysis, electronics, and many other fields of practical en gineering and science. How- ever, the Laplace transform is most properly interpreted as a particular real form of the more fundamental Fourier transform. When the Fourier trans form is evaluated along the imaginary axis, the complex exponential factor becomes rea l, and the result is the Laplace transform, which maps real-valued functions to real-value d functions. Since it is so closely allied to the Fourier transform, the Laplace transform enjo ys many of its featured proper- ties, including linearity. Moreover, derivatives are tran sformed into algebraic operations, which underlies its applications to solving differential eq uations. The Laplace transform is one-sided; it only looks forward in time and prefers funct ions that decay — transients. The Fourier transform looks in both directions and prefers o scillatory functions. For this reason, while the Fourier transform is used to solve boundar y value problems on the real line, the Laplace transform is much better adapted to initia l value problems. Since we will be applying the Laplace transform to initial va lue problems, we switch our variable from xtotto 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 /hatwidef(k) =1√ 2π/integraldisplay∞ 0f(t)e−iktdt, since, by our assumption, its negative tvalues make no contribution to the integral. The Laplace transform of such a function is obtained by replacing i kby a real†variable s, leading to F(s) =L[f(t)] =/integraldisplay∞ 0f(t)e−stdt, (13.137) where, in accordance with the standard convention, the fact or of√ 2πhas been omitted. ByallowingcomplexvaluesoftheFourierfrequency variabl ek, wemayidentifytheLaplace transform with√ 2πtimes the evaluation of the Fourier transform for values of k=−is on the imaginary axis: F(s) =√ 2π/hatwidef(−is). (13.138) Since the exponential factor in the integral has become real , the Laplace transform Ltakes real functions to real functions. Moreover, since the integ ral kernel e−stis exponentially decaying for s >0, we are no longer required to restrict our attention to func tions 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) =/integraldisplay∞ 0e(α−s)tdt=1 s−α. (13.139) †One can also define the Laplacetransform at complexvalues of s, but this willnot berequired in the applications discussed here. 12/11/12 729 c/circleco√yrt2012 Peter J. Olver Note that the integrand is exponentially decaying, and henc e the integral converges, if and only if Re( α−s)<0. Therefore, the Laplace transform (13.139) is, strictly s peaking, only defined at sufficiently large s >Reα. In particular, for an oscillatory exponential, L[eiωt] =1 s−iω=s+ iω s2+ω2provided s >0. Taking real and imaginary parts of this identity, we discove r the formulae for the Laplace transforms of the simplest trigonometric functions: L[cosωt] =s s2+ω2, L[sinωt] =ω s2+ω2. (13.140) Two additional important transforms are L[1] =/integraldisplay∞ 0e−stdt=1 s,L[t] =/integraldisplay∞ 0te−stdt=1 s/integraldisplay∞ 0e−stdt=1 s2.(13.141) The second computation relies on an integration by parts, ma king sure that the boundary terms at s= 0,∞vanish. Remark: In every case, we really mean the Laplace transform of the fu nction 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) =/braceleftbigg1, t > 0, 0, t < 0,(13.142) and so the first formula in (13.141) should more properly be wr itten L[σ(t)] =1 s. (13.143) However, in the traditional approach to the Laplace transfo rm, one only considers the functionsonthepositive taxis, andsothestepfunctionandtheconstantfunctionare, from this viewpoint, indistinguishable. However, once one move s beyond a purely mechanistic approach, any deeper understanding of the properties of the Laplace transform requires keeping this distinction firmly in mind. 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 oforderaif |f(t)|< M eat,for all t > t0, (13.144) for some M >0 andt0>0. Notethattheexponential growthconditiononlydepends upo nthefunction’s behavior for large values of t. Ifa <0, thenfis, in fact, exponentially decaying as x→ ∞. Since eat< ebtfora < band allt >0, iff(t) has exponential growth of order a, it automat- ically has exponential growth of any higher order b > a. All polynomial, trigonometric, 12/11/12 730 c/circleco√yrt2012 Peter J. Olver and exponential functions (with linear argument) have expo nential growth. The simplest example of a function that does not satisfy any exponential g rowth bound is f(t) =et2, since it grows faster than any simple exponential eat. The following result guarantees the existence of the Laplac e 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 appli cations. Theorem 13.33. Iff(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.144) implies that w e can bound the in- tegrand in (13.137) by |f(t)e−st|< 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 pra ctice, we may not be able to explicitly evaluate the Laplace transform integral. Nev ertheless, the Laplace transforms of most common functions are not hard to find, and extensive li sts have been tabulated, [143]. An abbreviated table of Laplace transforms can be found on the following page. Nowadays, the most convenient sources of transform formula s are computer algebra pack- ages, including Mathematica andMaple. According to Theorem 13.28, when it exists, the Fourier tran sform uniquely specifies the function, except possibly at jump discontinuities wher e the limiting value must be half way in between. An analogous result can be established for th e Laplace transform. Lemma 13.34. Iff(t)andg(t)are piecewise continuous functions that are of ex- ponential growth, and L[f(t)] =L[g(t)]for allssufficiently 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.138), with the Fourier transfo rm 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/integraldisplayi∞ −i∞F(s)estds, t > 0. (13.145) In practice, one hardly ever uses this complicated formula t o compute the inverse Laplace transform. Rather, one simply relies on tables of known Lapl ace transforms, coupled with a few basic rules that will be covered in the following subsec tion. †See Section 16.5 for details on complex integration. The stated formula do esn’t apply to all functions of exponential growth. A more universally valid inverse Lap lace transform formula is obtained by shifting the complex contour to run from b−i∞tob+ i∞for some b > a, the order of exponential growth of f. 12/11/12 731 c/circleco√yrt2012 Peter J. Olver Table of Laplace Transforms f(t) F(s) 11 s t1 s2 tn n! sn+1 δ(t−c) e−sc eαt 1 s−α cosωts s2+ω2 sinωtω s2+ω2 ectf(t) F(s−c) σ(t−c)f(t−c) e−scF(s) tf(t) −F′(s) f′(t) sF(s)−f(0) f(n)(t) snF(s)−sn−1f(0)− −sn−2f′(0)−···−f(n−1)(0) f(t)∗g(t) F(s)G(s) In this table, nis a non-negative integer, ωis any real number, while c≥0 is any non-negative real number. 12/11/12 732 c/circleco√yrt2012 Peter J. Olver The Laplace Transform Calculus The first observation is that the Laplace transform is a linea r operator, and so L[f+g] =L[f]+L[g], L[cf] =cL[f], (13.146) for any constant c. Moreover, just like its Fourier progenitor, the Laplace tr ansform converts calculus into algebra. In particular, differentia tion turns into multiplication by the transform variable, but with one additional term that de pends upon the value of the function at t= 0. Theorem 13.35. Letf(t)be continuously differentiable for t >0and have expo- nential growth of order a. IfL[f(t)] =F(s)then, for s > a, L[f′(t)] =sF(s)−f(0). (13.147) Proof: The proof relies on an integration by parts: L[f′(t)] =/integraldisplay∞ 0f′(t)e−stdt=f(t)e−st/vextendsingle/vextendsingle/vextendsingle/vextendsingle∞ t=0+s/integraldisplay∞ 0f(t)e−stdt = lim t→∞f(t)e−st−f(0)+sF(s). The exponential growth inequality (13.144) implies that fir st term vanishes for s > a, and the remaining terms agree with (13.147). Q.E.D. Example 13.36. According toExample13.31, theLaplace transform of thefun ction sinωtis L[sinωt] =ω s2+ω2. Its derivative is d dtsinωt=ωcosωt, and therefore L[ωcosωt] =sL[sinωt] =ωs s2+ω2, since sin ωtvanishes at t= 0. The result agrees with (13.140). On the other hand, d dtcosωt=−ωsinωt, and so L[−ωsinωt] =sL[cosωt]−1 =s2 s2+ω2−1 =−ω2 s2+ω2, again in agreement with the known formula. 12/11/12 733 c/circleco√yrt2012 Peter J. Olver Remark: The final term −f(0) in (13.147) is a manifestation of the discontinuity in f(t) att= 0. Keep in mind that the Laplace transform only applies to fu nctions with f(t) = 0 for all t <0, and so if f(0)/negationslash= 0, then f(t) has a jump discontinuity of magnitude f(0)att= 0. Therefore, bythecalculusofgeneralizedfunctions, it sderivative 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 Lapla ce transform calculus, one suppresses the delta function when computing the derivativ ef′(t). However, its effect must reappear on the other side of the differentiation formul a (13.147), and the upshot is the extra term −f(0). Laplace transforms of higher order derivatives are found by iterating the first order formula (13.147). For example, if f∈C2, then L[f′′(t)] =sL[f′(t)]−f′(0) =s2F(s)−sf(0)−f′(0). (13.148) In general, for an ntimes continuously differentiable function, L[f(n)(t)] =snF(s)−sn−1f(0)−sn−2f′(0)− ··· − f(n−1)(0). (13.149) On the other hand, multiplying the function by tcorresponds to differentiation of its Laplace transform (up to a change of sign): L[tf(t)] =−F′(s). (13.150) The proof of this formula is left as an exercise for the reader . Conversely, integration corresponds to dividing the Lapla ce transform by s, so L/bracketleftbigg/integraldisplayt 0f(τ)dτ/bracketrightbigg =F(s) s. (13.151) 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[t2] =1 sL[2t] =2 s3,and, more generally, L[tn] =n! sn+1.(13.152) There is also a shift formula, analogous to Theorem 13.19 for Fourier transforms, but with one important caveat. Since all functions must vanish f ort <0, we are only allowed to shift them to the right, a shift to the left would produce no nzero function values for somet <0. In general, the Laplace transform of the function f(t−c) shifted to the right by an amount c >0, is L[f(t−c)] =/integraldisplay∞ 0f(t−c)e−stdt=/integraldisplay∞ −cf(t)e−s(t+c)dt (13.153) =/integraldisplay0 −cf(t)e−s(t+c)dt+/integraldisplay∞ 0f(t)e−s(t+c)dt=e−sc/integraldisplay∞ 0f(t)e−stdt=e−scF(s). In this computation, we first used a change of variables in the integral, replacing t−cby t; then, the fact that f(t)≡0 fort <0 was used to eliminate the integral from −cto 0. When using the shift formula in practice, it is important to k eep 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. 12/11/12 734 c/circleco√yrt2012 Peter J. Olver Example 13.37. Consider the square wave pulse f(t) =/braceleftbigg1, b < t < c, 0,otherwise ,for some 0< b < c. To compute its Laplace transform, we write it as the differen ce f(t) =σ(t−b)−σ(t−c) of shifted versions of the step function (13.142). Combinin g the shift formula (13.153) and the formula (13.143) for the Laplace transform of the step fu nction, we find L[f(t)] =L[σ(t−b)]−L[σ(t−c)] =e−sb−e−sc s. (13.154) We already noted that the Fourier transform of the convoluti on product of two func- tions is realized as the ordinary product of their individua l 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.128) reduces to a finite integral h(t) =f(t)∗g(t) =/integraldisplayt 0f(t−τ)g(τ)dτ. (13.155) 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. IfL[f(t)] =F(s)andL[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 transfo rm proceeds along the same lines as its Fourier transform version Theorem 13.26, a nd is left as Exercise for the reader. Applications to Initial Value Problems The keyapplicationoftheLaplace transform istofacilitat ethe solutionofinitialvalue problems for linear, constant coefficient ordinary different ial equations. As a prototypical example, consider the second order initial value problem ad2u dt2+bdu dt+cu=f(t), u (0) =α,du dt(0) =β, (13.156) inwhich a,b,careconstant. Wewillsolvetheinitialvalueproblem by appl yingtheLaplace transform to both sides of the differential equation. In view of the differentiation formulae (13.147,148), a/parenleftbig s2L[u(t)]−su(0)−/squaresmallsolidu(0)/parenrightbig +b/parenleftbig sL[u(t)]−u(0)/parenrightbig +cL[u(t)] =L[f(t)]. SettingL[u(t)] =U(s) andL[f(t)] =F(s), and making use of the initial conditions, the preceding equation takes the form (as2+bs+c)U(s) =F(s)+(as+b)α+aβ. (13.157) 12/11/12 735 c/circleco√yrt2012 Peter J. Olver Thus, by applying the Laplace transform, we have effectively reduced the entire initial value problem to a single elementary algebraic equation! So lving for U(s) =F(s)+(as+b)α+aβ as2+bs+c, (13.158) 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 fou nd by suitably massaging the answer (13.158) to be a combination of known tr ansforms. Example 13.39. Consider the initial value problem /squaresmallsolid/squaresmallsolidu+u= 10e−3t, u (0) = 1,/squaresmallsolidu(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 s2+1+10 (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 partia l fractions to write it as a sum U(s) =s+2 s2+1+1 s+3+3−s s2+1=1 s+3+5 s2+1 of terms appearing in the table. We conclude that the solutio n to our initial value problem is u(t) =e−3t+5sint. Of course, the last example is a problem that you can easily so lve directly. The stan- dard method learned in your first course on differential equat ions is just as effective in finding the final solution, and does not require all the extra L aplace transform machin- ery! The Laplace transform method is, however, particularl y 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 uni t force over time interval 1 2π < t <2π, after which it left to vibrate on its own. The initial value p roblem is d2u dt2+4u=/braceleftigg 1,1 2π < t <2π, 0,otherwise ,u(0) =/squaresmallsolidu(0) = 0. Taking the Laplace transform of the differential equation, a nd using (13.154), 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.153) u(t) =h/parenleftbig t−1 2π/parenrightbig −h(t−2π), 12/11/12 736 c/circleco√yrt2012 Peter J. Olver whereh(t) is the function with Laplace transform L[h(t)] =H(s) =1 s(s2+4)=1 4/parenleftbigg1 s−s s2+4/parenrightbigg , which has been conveniently rewritten using partial fracti ons. Referring to our table of Laplace transforms, h(t) =/braceleftbigg1 4−1 4cos2t, t > 0, 0, t < 0. Therefore, our desired solution is u(t) =  0, 0≤t≤1 2π, 1 4+1 4cos2t,1 2π≤t≤2π, 1 2cos2t, 2π≤t. Note that the solution u(t) is only C1at 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 t he forcing function, leading to a solution on that interval depending upon two integration c onstants. The integration constants are then adjusted so that the solution satisfies th e initial conditions and is continuouslydifferentiableateachpointofdiscontinuity oftheforcingfunction. Thedetails are straightforward, but messy. The Laplace transform meth od successfully bypasses these intervening manipulations. Example 13.41. Consider the initial value problem d2u dt2+ω2u=f(t), u (0) =/squaresmallsolidu(0) = 0, involvingageneralforcing function f(t). ApplyingtheLaplacetransform tothedifferential equation and using the initial conditions, (s2+ω2)U(s) =F(s),and hence U(s) =F(s) s2+ω2. The right hand side is the product of the Laplace transform of the forcing function f(t) and that of the trigonometric function k(t) =sinωt ω. Therefore, Theorem 13.38 implies that the solution can be written as their convolution u(t) =f∗k(t) =/integraldisplayt 0k(t−τ)f(τ)dτ=/integraldisplayt 0sinω(t−τ) ωf(τ)dτ. (13.159) where k(t) =/braceleftiggsinωt ω, t > 0, 0, t < 0. 12/11/12 737 c/circleco√yrt2012 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 i mpulse force that is applied instantaneously atthetime t=τ. Noteparticularlythat(unlikeboundary valueproblems) 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.159) can be viewed as a linear superposition of fundamental solution responses induced by expressing the f orcing function f(t) =/integraldisplay∞ 0f(τ)δ(t−τ)dτ as a superposition of concentrated impulses over time. This concludes our brief introduction to the Laplace transf orm and a few of its many applications to physical problems. More details can be foun d in almost all applied texts on mechanics, electronic circuits, signal processing, con trol theory, and many other areas. 12/11/12 738 c/circleco√yrt2012 Peter J. Olver