An Overview Of Wavelet Transform Concepts And Applications

An overview of wavelet transform concepts and applicationsChristopher Liner, University of HoustonFebruary 26, 2010AbstractThe continuous wavelet transform utilizing a complex Morlet analyzing wavelet has aclose connection to the Fourier transform and is a powerful analysis tool for decomposingbroadband wavefield data. A wide range of seismic wavelet applications have beenreported over the last three decades, and the free Seismic Unix processing system nowcontains a code (succwt) based on the work reported here.IntroductionThe continuous wavelet transform (CWT) is one method of investigating the time-frequencydetails of data whose spectral content varies with time (non-stationary time series). Motivation for the CWT can be found in Goupillaud et al. [12], along with a discussion of itsrelationship to the Fourier and Gabor transforms.As a brief overview, we note that French geophysicist J. Morlet worked with nonstationary time series in the late 1970’s to find an alternative to the short-time Fouriertransform (STFT). The STFT was known to have poor localization in both time and frequency, although it was a first step beyond the standard Fourier transform in the analysis ofsuch data. Morlet’s original wavelet transform idea was developed in collaboration with theoretical physicist A. Grossmann, whose contributions included an exact inversion formula.A series of fundamental papers flowed from this collaboration [16, 12, 13], and connectionswere soon recognized between Morlet’s wavelet transform and earlier methods, includingharmonic analysis, scale-space representations, and conjugated quadrature filters. For further details, the interested reader is referred to Daubechies’ [7] account of the early historyof the wavelet transform.Here we describe implementation details of the CWT as applied to wavefield measurements, using reflection seismic data as a specific example. A classic paper in the general field of time-frequency decomposition of reflection seismic data is Chakraborty andOkaya [4]. They discuss the short-time Fourier transform, continuous wavelet transform,discrete wavelet transform, and matched pursuit decomposition, but application to realseismic data is limited to the matched pursuit method.The CWT implementation described here uses a complex Morlet analyzing wavelet. Asa time-domain gaussian-tapered complex exponential, this particular wavelet has a naturaland compelling connection to the venerable Fourier transform. An inventory of several1

other analyzing wavelets can be found in the excellent summary paper of Torrence andCompo [35].Fourier TransformWe will use the convention that a time function, g(t), and the Fourier Transform (FT) ofthat function, g(ω), are in the time or frequency domain as indicated by the argument listrather than some variation on the function symbol.The forward FT is defined as usualZ g(t) eiωt dtg(ω) ,(1) where scaling constants have been omitted, i 1, and ω is angular frequency related tolinear frequency (f in Hz) by ω 2πf . The function g(ω) is complex, and can be expressedin polar form g(ω) Aeiθ to reveal the amplitude spectrum, A(ω), and phase spectrum,θ(ω). The inverse FT is given byZ g(t) g(ω) e iωt dω.(2) We define the Dirac delta function heuristically as the unit spike function(δ(t t0 ) 1,0,t t0t 6 t0.(3)A key feature of δ(t) is the sifting property it exhibits under the action of integrationZ δ(t t0 ) g(t) dt g(t0 ),(4)It is customary to remark that the FT decomposes a transient time signal (data) intoindependent harmonic components and, therefore, the function g(ω) has exact frequencylocalization and no time localization. In other words, the FT can precisely detect whichfrequencies reside in the data, but yields no information about the time position of signalfeatures. We can investigate this claim by considering the FT impulse response.Let g(t) δ(t t0 ), take the FT, and apply the sifting property to yieldZg(ω) δ(t t0 ) eiωt dt eiωt0.(5)Recognizing the result is already in polar form, we see the amplitude spectrum, A 1,contains no information about the time location of the spike. The phase, however,θ ω t0 2πf t02,(6)

is linear, and the spike location, t0 , is encoded in the phase slope. For this elementary casethe spike location can be recovered exactly byt0 1 dθ2π df.(7)More challenging is the case of two spikesg(t) δ(t t1 ) δ(t t2 )(8)that transforms tog(ω) eiωt1 eiωt2(9)and which can be written in terms of t t2 t1 as g(ω) eiωt1 1 eiω t .(10)The accessible spike time information in this function lies in the amplitude rather than thephase spectrum. Specifically, consider the zeros of the amplitude spectrum g(ω) eiωt11 eiω t 0.(11)By definition, the first right-side term is never zero, meaning the equality can only hold ifeiω t 1(12)iω t ln( 1)(13)i2πfn t i(2n 1)π2n 1fn ; n 0, 1, 2, .2 t(14)(15)where fn is the nth notch frequency. This relationship says that if we can observe thefrequency associated with, say, the first spectral notch, f0 , then we can calculate the timeseparation of the spikes using t t2 t1 12f0.(16)This accomplishment is a muted victory since we do not find the absolute spike times, onlythe delay between them. Furthermore, even this weak result breaks down if the spikes havedifferent amplitudes (no hard zeros develop in the spectrum), or we move on to the threeor N-spike case. We conclude that although there is some time localization information inthe Fourier transform, it quickly becomes an unreasonable exercise to extract it for evensimple impulsive time functions.Something better is needed, and that something is a wavelet transform.3

Wavelet TransformThe continuous wavelet transform can be defined in a variety of ways that differ with respectto normalization constants and conjugation. We define the transform asg(a, b) apZ g(t) ψ t ba dt,(17)where ψ(t) is the analyzing wavelet, b is a time-like translation variable, a is a dimensionlessfrequency scale variable, and p is a real normalization parameter. As with the FT, we usethe argument list to indicate whether g() is in the physical domain, g(t), or the transformdomain, g(a, b). This convention is routinely used in geophysics where multidimensionaltransforms are often encountered [6, 23].For reconstruction of the original time series, the inverse CWT is in principle given bythe double integralZ Z g(t) g(a, b) ψ 0 t ba da db,(18)where the inverse analyzing wavelet ψ 0 need not be the same as the forward transformwavelet. As discussed by Torrence and Compo [35], if the inverse wavelet is chosen to be adelta function the b integral can be done analytically via the sifting property of the deltafunction. The inverse is then simplified to the real part of the summation over scales Z g(t) Reg(a, t) da.(19) Complex Morlet WaveletAs discussed earlier, the Fourier transform kernel is given byKF T ei2πf t.(20) where i 1 and f is frequency in Hertz. The kernel in a continuous wavelet transformis simply the analyzing wavelet,KCW T ψ(t) .(21)For wave-like data a good choice for the analyzing wavele is the complex Morlet wavelet2ψ(t) e (t/c) ei2πf0 t,(22)where t is time, f0 is a frequency parameter, and c is a damping parameter with units oftime. This equation describes a time-domain function that is the product of a gaussian anda complex exponential, and whose center frequency is f0 .The frequency domain representation of the complex Morlet wavelet is2ψ(f ) e (cf ) δ(f f0 )4,(23)

where normalization factors have been omitted, and denotes convolution. This is, again,a gaussian function, now centered in the frequency domain on f0 .There are several choices that can be made with respect to wavelet normalization. Ourdefinition, equation 22 along with p 0 in equation 17, normalizes the time-domain peakamplitude.Comparing the leading terms in equations 22 and 23, we see a duality that expressesthe characteristic CWT resolution trade-off. The damping parameter, c, controls the rateat which the time-domain wavelet, and frequency domain spectrum, is driven toward zero.As time localization increases (large c), frequency localization decreases, and vice versa.We have chosen to write the complex Morlet wavelet in this particular form and notationto emphasize the central role of the damping factor. The value of this parameter has afirst-order effect on any CWT result. In fact, in the limit as c the Morlet waveletdefined here becomes equal to the Fourier kernel.Figure 1 illustrates the time-frequency resolution properties of the complex Morletwavelet. In this example, a 60 Hz Morlet wavelet (A) is shown with the real part as asolid line and imaginary part dashed. The damping parameter is c 1/120, resulting ina gaussian-tapered amplitude spectrum (B). Lowering the center frequency to 30 Hz andusing c 1/60 gives a wider time-domain wavelet (C) and a narrower frequency spectrum(D). It should be noted, that these results are not dependent on the absolute size of thefrequencies being considered, similar plots could be constructed for 60 and 30 MHz.Note the choice of damping parameter in Figure 1 as the inverse of twice the centerfrequency of the wavelet. Using this value will preserve a wavelet time span equal two timesits period.When used to compute the CWT,the argument of the analyzing wavelet (equation 22) t bis translated and scaled, ψ a . The effect of scaling is two-fold. In the real exponential,the scale value multiplies the damping parameter, effectively broadening the time-domainextent of the wavelet. In the complex exponential, the scale divides the frequency, loweringit precisely enough to maintain the number of time-domain oscillations.Scales and FrequenciesThe continuous wavelet transform represented by equation 17 is related to convolution asevidenced by the appearance of (t b) in the integral. Thus, the CWT can be written asg(a, t) g(t) ψ(t/a),(24)showing the CWT can be implemented as a series of complex-valued time domain convolutions. Each convolution involves the input data and a version of the analyzing waveletmodified by the scaling variable. Strictly speaking, the CWT is an inner product, or thezero lag of a complex-valued correlation, that can be expressed as a convolution under certain symmetry conditions. The CWT as a cross-correlation is understandable since we arematching similarities between the signal and the analyzing wavelet.From a computational point of view, complex time-domain convolutions are highly inefficient. The mathematical form of equation 24 suggests implementation in the Fourier5

Figure 1: Complex Morlet wavelet in the time and frequency domain. (A) 60 Hz Morlet waveletshowing real (solid) and imaginary (dash) parts. The damping parameter is c 1/120. (B) Fourieramplitude spectrum of the 60 Hz Morlet wavelet showing characteristic gaussian shape. (C) A 30 HzMorlet wavelet (c 1/60) has longer time duration, but the same number of oscillations. (D) Theamplitude spectrum of the 30 Hz wavelet is narrower because the time-domain function is wider.domain where the convolution will become multiplication. Recalling the Fourier transformscaling property for a general time function g(t)F T {g(t/a)} a g(af ) ,(25)we can write equation 24 in the Fourier domain as the compact and efficient resultg(a, t) F T 1 {a g(f ) ψ(af )}.(26)A fundamental contribution of Morlet’s early work [12] was recognition that the naturalsampling of this scale variable is dyadic (i.e., logarithmic, base 2). If the wavelet is dilatedby a factor of 2, this means the frequency content has been shifted by one octave. Theanalogy with music and singing clear, and is perpetuated by defining the scale range interms of octaves and voices. The number of octaves determines the span of frequenciesbeing analyzed, while the number of voices per octave determines the number of samples(scales) across this span.Specifically, let the number of octaves in a CWT be No and the number of voices per6