VIEWS: 9 PAGES: 30 POSTED ON: 8/1/2012
March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets WAVELET ANALYSIS AND SOME OF ITS APPLICATIONS IN PHYSICS J.-P. ANTOINE e e Institut de Physique Th´orique, Universit´ Catholique de Louvain B-1348 Louvain-la-Neuve, Belgium E-mail: Antoine@fyma.ucl.ac.be We review the general properties of the wavelet transform, both in its continuous and its discrete versions, in one or two dimensions, and we describe some of its applications in signal and image processing. We also consider its extension to higher dimensions. 1. What is wavelet analysis? Wavelet analysis is a particular time- or space-scale representation of signals which has found a wide range of applications in physics, signal processing and applied mathematics in the last few years. In order to get a feeling for it and to understand its success, let us consider ﬁrst the case of one- dimensional signals. It is a fact that most real life signals are nonstationary. They often contain transient components, sometimes very signiﬁcant physically, and mostly cover a wide range of frequencies. In addition, there is frequently a direct correlation between the characteristic frequency of a given segment of the signal and the time duration of that segment. Low frequency pieces tend to last a long interval, whereas high frequencies occur in general for a short moment only. Human speech signals are typical in this respect. Vowels have a relatively low mean frequency and last quite long, whereas consonants contain a wide spectrum, up to very high frequencies, especially in the attack, but they are very short. Clearly standard Fourier analysis is inadequate for treating such signals, since it loses all information about the time localization of a given frequency component. In addition, it is very uneconomical. If a segment of the signal is almost ﬂat, thus carries little information, one still has to sum an (alternating) inﬁnite series for reproducing it. Worse yet, Fourier analysis is 1 March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 2 2 4 Figure 1. A traditional time-frequency representation of a signal (from Mozart’s Don Giovanni, Act 1). highly unstable with respect to perturbation, because of its global character. For instance, if one adds an extra term, with a very small amplitude, to a linear superposition of sine waves, the signal will barely be modiﬁed, but the Fourier spectrum will be completely perturbed. Such pathologies do not occur if the signal is represented in terms of localized components. Therefore, signal analysts turn to time-frequency (TF) representations. The idea is that one needs two parameters. One, called a, characterizes the frequency, the other one, b, indicates the position in the signal. This concept of a TF representation is in fact quite old and familiar. The most obvious example is simply a musical score (Figure 1)! If one requires, in addition, the transform to be linear, a general TF transform will take the form: ∞ s(x) → S(b, a) = ψb,a (x) s(x) dx, (1.1) −∞ where s is the signal and ψb,a the analyzing function (we denote the time variable by x, in view of the extension to higher dimensions). The as- sumption of linearity is actually nontrivial, for there exists a whole class of quadratic or, more properly, sesquilinear time-frequency representations. The prototype is the so-called Wigner-Ville transform, introduced origi- nally by E.P. Wigner1 in quantum mechanics (in 1932!) and extended by J. Ville2 to signal analysis: +∞ x x Ws (b, ξ) = dx e−iξx s(b − ) s(b + ), ξ = 1/a. (1.2) −∞ 2 2 Notice that, contrary to the linear version (1.1), this transform is intrinsic, it contains no external probe ψ, which inavoidably inﬂuences the result. Within the class of linear TF transforms, two stand out as particularly simple and eﬃcient, the windowed or short time Fourier transform (WFT) and the wavelet transform (WT). For both of them, the analyzing function ψb,a is obtained by a group action on a basic (or mother) function ψ, only the group diﬀers. The essential diﬀerence between the two is in the way March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 3 1/a ≈ frequency ψb,a (x) a<1 high a=1 medium a>1 low x Figure 2. The function ψb,a (x) for diﬀerent values of the scale parameter a : in the case of the Windowed Fourier Transform (left); in the case of the wavelet transform (right). the frequency parameter a is introduced. For the WT, one takes: 1 x−b ψb,a (x) = √ ψ . (1.3) a a The action of a on the function ψ is a dilation (a > 1) or a contraction (a < 1): The shape of the function is unchanged, it is simply spread out or squeezed. As for b, it is simply a translation. By contrast, the WFT takes for ψb,a the function ψb,a (x) = eix/a ψ(x − b). This means that the a-dependence is a modulation (1/a ∼ frequency); the window has constant width, but the lower a, the larger the number of oscillations in the window ψ. All this is illustrated on Figure 2. March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 4 Actually one should distinguish between two radically diﬀerent versions of the wavelet transform, the continuous WT (CWT) and the discrete WT (DWT), although they are based on the same basic formula, namely: ∞ x−b S(b, a) = a−1/2 ψ s(x) dx, a > 0, b ∈ R. (1.4) −∞ a In this relation, s is a ﬁnite energy signal and S(b, a) is its WT with respect to the analyzing wavelet ψ. Now, in the CWT, all values of a and b are o considered. The CWT therefore plays the same rˆle as the Fourier trans- form and is mostly used for analysis and feature detection in signals. Of course, the CWT must be discretized for its numerical implementation, but the user may choose the sampling grid according to his needs. The price to pay, however, is that no orthonormal bases are obtained in this way, only frames: the CWT is a redundant representation. This may be seen as a defect, but redundancy actually has the advantage of ensuring the stability of the representation. The DWT, on the other hand, is based on a preselected grid (often dyadic) and is explicitly designed for generating (bi)orthonormal bases, starting from multiresolution analysis. It may be viewed as the analogue of the Discrete Fourier Transform (see, for instance, Ref. 3) and is more appropriate for data compression and signal reconstruction. As a matter of fact, the (discretized) CWT is incompatible with the DWT, they are based on totally diﬀerent philosophies and also have diﬀerent purposes. In these lectures, we will review the CWT, both from the theoretical side and in its practical applications, in one and two dimensions, plus some extensions. We will also give some indications on the Discrete WT. Finally a word about references. As a ﬁrst contact, introductory articles such as Refs. 3 or 4 may be a good suggestion, followed by the popular book of B. Burke Hubbard5 or the elementary book of Y. Meyer.6 As for textbooks, we note in particular the celebrated volume of I. Daubechies7 and the treatise of S. Mallat.8 Our main source of information, however, especially concerning the applications of the CWT in physics, will be van den Berg’s survey9 and our own upcoming monograph.10 In addition, we ought to mention the popular newsletter The Wavelet Digest, originally founded 12 years ago by Wim Sweldens, and nowadays edited by Michael Unser at the Swiss Federal Institute of Technology in Lausanne, Switzerland (EPFL). The circulation has by now passed the 20,000 mark and is still growing! For further information, we refer to the website <http://www.wavelet.org/>. March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 5 2. The one-dimensional CWT 2.1. Basic formulas As announced above, the basic formulas read, in time domain and in fre- quency domain, respectively, ∞ S(b, a) = a−1/2 dx ψ (a−1 (x − b)) s(x) (2.1) −∞ ∞ = a1/2 dω ψ(aω) s(ω) eiωb , (2.2) −∞ where a > 0 is a scale parameter, b ∈ R is a translation parameter and the hat denotes a Fourier transform. Thus the pair (b, a) runs over the time-scale half-plane R2 . Sometimes one uses also the full range a = 0, but + this is physically less natural. In these relations, s is a square integrable function and the analyzing wavelet ψ assumed to be well localized both in the space (or time) domain and in the frequency domain. Here the minimal requirement is that ψ ∈ L1 ∩L2 , but in practice stronger conditions are usually imposed (like ψ ∈ S, Schwartz’s space of fast decreasing functions). In consequence, the CWT provides good bandpass ﬁltering both in x and in ω. Moreover, ψ must satisfy the following admissibility condition, which guarantees the invertibility of the WT (see (2.13) below): ∞ |ψ(ω)|2 cψ ≡ 2π dω < ∞. (2.3) −∞ |ω| In most cases, this condition may be reduced to the requirement that ψ has zero mean (hence it must be oscillating): ∞ ψ(0) = 0 ⇐⇒ dx ψ(x) = 0. (2.4) −∞ In addition, ψ is often required to have a certain number of vanishing moments: ∞ dx xn ψ(x) = 0, n = 0, 1, . . . N. (2.5) −∞ This property improves the eﬃciency of ψ at detecting singularities in the signal, since it is blind to polynomials up to order N . Finally, it is often useful, but not essential, that ψ be progressive, which means that ψ is real and ψ(ω) = 0 for ω < 0 (in the language of signal processing, ψ is an analytic signal). March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 6 −10 0 10 −10 0 10 Figure 3. (left) The Mexican hat or Marr wavelet; (right) The real part of the Morlet wavelet, for ωo = 5.6. In order to ﬁx ideas, we exhibit here two commonly used wavelets (see Figure 3). (1) The Mexican hat or Marr wavelet This wavelet is simply the second derivative of a Gaussian: 2 2 ψH (x) = (1 − x2 ) e−x /2 , ψH (ω) = ω 2 e−ω /2 . (2.6) It is real and admissible, it has 2 vanishing moments n = 0, 1, but it is not progressive. (2) The Morlet wavelet This wavelet is essentially a plane wave within a Gaussian window: 2 2 ψM (x) = eiko x e−x /2σo + c(x) (2.7) −[(ω−ωo )σo ]2 /2 ψM (ω) = σo e + c(ω). (2.8) It is complex, but not progressive, strictly speaking (numerically it is). Here the correction term c must be added in order to satisfy the admissibility condition (2.4), but in practice one will arrange that this term be numerically negligible ( 10−4 ) and thus can be omitted (it suﬃces to choose the basic frequency |ωo | large enough, namely, one has to take |σo ωo | > 5.5). These two wavelets have very diﬀerent properties and, naturally, they will be used in quite diﬀerent situations. Typically, the Mexican hat is sensitive to singularities in the signal, and it yields a genuine time-scale analysis. On the other hand, since the Morlet wavelet is complex, it will catch the phase of the signal, hence will be sensitive to frequencies, and will lead to a time-frequency analysis, somewhat closer to a Gabor analysis. In both cases, additional ﬂexibility is obtained by adding a width parameter to the Gaussian. March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 7 ω ∼ 1/a 6 a<1: ωo /a – Ω/a aL a=1: ωo – Ω L a>1: ωo /a – Ω/a aL - b x Figure 4. d Support properties of ψb,a and ψb,a . 2.2. Localization properties and interpretation We must now make more precise the localization conditions on the wavelet ψ. Assume that the numerical support of ψ(x) is an interval of length L around 0, and that of ψ(ω) is an interval of length Ω, centered around the mean frequency ωo (by numerical support, we mean the smallest set outside of which the function is numerically negligible). Then the numerical support of ψb,a (x) is an interval of length aL around b, while that of ψ(ω) is an interval of length Ω/a, centered around ωo /a. Therefore (see Figure 4): • if a 1, ψb,a is a wide window (long duration) and ψb,a is peaked around the small frequency ωo /a: for large scales, the CWT is sensitive to low frequencies, and thus yields a rough analysis. • if a 1, ψb,a is a narrow window (short duration), and ψb,a is wide and centered around the high frequency ωo /a: for very small scales, the CWT is sensitive to high frequencies (small details). These support properties are illustrated in Figure 5 for the case of the Morlet wavelet ψM . Combining now these localization properties with the fact that the CWT is a convolution with a zero mean ﬁlter, we conclude that: March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 8 −10 0 10 −10 0 10 −10 0 10 0 10 20 0 10 20 0 10 20 Figure 5. Support properties of the Morlet wavelet ψM : for a = 0.5, 1, 2 (left to right), d ψb,a has width 3, 6, 12, respectively (top), while ψb,a has width 3, 1.5, 0.75, and peaks at 12, 6, 3 (bottom). • the CWT provides a local ﬁltering in time (b) and scale (a): S(b, a) ≈ 0 ⇐⇒ ψb,a (x) ≈ s(x). • the CWT may be interpreted as a mathematical microscope, with optics ψ, position b, and magniﬁcation 1/a. • the CWT works at constant relative bandwidth : ∆ω/ω = const. As a consequence, one may consider the CWT as a singularity detector and analyzer. 2.3. Mathematical properties Given an admissible wavelet ψ, the CWT Wψ : s(x) → S(b, a) is a linear map, with the following properties: (1) Covariance under translation and dilation: Wψ : s(x − xo ) → S(b − xo , a) (2.9) 1 x b a Wψ : √ s →S , . (2.10) ao ao ao ao March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 9 (2) Energy conservation : ∞ da db dx |s(x)|2 = c−1 ψ |S(b, a)|2 . (2.11) −∞ R2 + a2 Thus, |S(b, a)|2 may be interpreted as the energy density in the time-scale half-plane. The relation (2.11) means that Wψ is an isometry from the space of signals L2 (R) onto a closed subspace Hψ of L2 (R2 , da db/a2 ), namely, the + space of wavelet transforms. An equivalent statement is that the wavelet ψ generates a resolution of the identity: da db c−1 ψ |ψba ψba | = I. (2.12) R2 + a2 (3) Reconstruction formula : As a consequence, Wψ is invertible on its range Hψ by the adjoint map, i.e., we obtain an exact reconstruction formula: da db s(x) = c−1 ψ ψb,a (x) S(b, a). (2.13) R2 + a2 In other words, we have a representation of the signal s(x) as a linear superposition of wavelets ψb,a with coeﬃcients S(b, a). (4) The projection Pψ : L2 (R2 , da db/a2 ) → Hψ is an integral operator, + with kernel K(b , a ; b, a) = c−1 ψb a |ψba . ψ (2.14) The kernel K is the autocorrelation function of ψ and it is a reproducing kernel. Indeed, the function f ∈ L2 (R2 , da db/a2 ) is the WT of a certain + signal iﬀ it satisﬁes the reproduction property da db f (b , a) = c−1 ψ ψb ,a |ψb,a f (b, a). (2.15) R2 + a2 This means that the CWT is a highly redundant representation! The last statement implies that the full information must be contained is a small subset of half-plane. This property may be exploited in two ways: either one takes a discrete subset, which leads to the theory of frames, or one considers only the lines of local maxima, called ridges. We will explore these two approaches in succession. March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 10 2.4. Discretization of the CWT : Frames We will say that a discrete lattice Γ = {aj , bk , j, k ∈ Z} yields a good discretization if s(x) = ψjk |s ψjk (x), ∀ s ∈ L2 (R), (2.16) j,k ∈ Z where ψjk ≡ ψbk ,aj and ψjk is explicitly constructible from ψjk (the relation is still exact). The relevant concept here is that of a frame, that is, a family {ψjk } for which there exists two constants A > 0, B < ∞ such that 2 A s | ψjk |s |2 B s 2 , ∀ s ∈ L2 (R). (2.17) j,k ∈ Z The upper bound means that the map s → { ψjk |s } is continuous from L2 to l2 , whereas the lower bound guarantees the numerical stability of the inverse map.7 The constants A, B are called frame bounds. If A = B = 1, the frame is said to be tight. Of course, if A = B = 1 and ψjk = 1, ∀ j, k, we simply get an orthonormal basis. A detailed analysis then shows that the expansion (2.16) converges essentially as a power series in the quantity B/A − 1. Now the question is, given the wavelet ψ, can one ﬁnd a lattice Γ such that {ψjk } is a good frame ? By good frame, we mean a frame such that (2.16) converges fast, and thus can be truncated after very few terms. This requires B/A − 1 1. In order to achieve this, one usually chooses a lattice adapted to the geometry of the time-scale half-plane, for instance, the dyadic lattice ψjk (x) = 2j/2 ψ(2j x − k), j, k ∈ Z. A direct estimation shows that both the Mexican hat and the Morlet wavelets give good, but nontight frames.7 2.5. Reducing the computational cost : Ridges Real life signals are often entangled and noisy, which makes the WT diﬃcult to interpret. But the energy density |S(b, a)|2 is usually well concentrated, around lines of local maxima, called ridges. A fundamental result (based on stationary phase arguments, the integral in (1.4) being rapidly oscillating) is that the restriction of the WT S(b, a) to the set of its ridges contains essentially the whole information. Thus it suﬃcient to consider this re- striction, called the skeleton of the WT, which reduces substantially the March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 11 1 1 1.5 1.5 Log Scale Log Scale 2 2 2.5 2.5 3 3 100 200 300 400 500 600 700 800 900 1000 100 200 300 400 500 600 700 800 900 1000 Position Position Figure 6. CWT of an academic discontinous signal (left) and the corresponding skele- ton, showing the vertical ridges pointing towards each singularity (right); the signal is shown at the bottom of each panel. computational cost, and also simpliﬁes considerably the interpretation of the WT. Actually one may distinguish several types of ridges. Two particular cases are the vertical ridges, characteristic of singularities in the signal, and horizontal ridges, that signal the appearance of given prominent frequencies (see Refs. 11 and 12 for a full discussion). Let us discuss the ﬁrst type. Assume s(x) has a singularity of order α −1 at xo (here we take α ∈ Z, but the analysis can be extended to any α ∈ R): s(x) ∼ θ(x − xo ) (x − xo )α for x ∼ xo , so that dα+1 s(x) ∼ δ(x − xo ) (most singularities may be modeled in this x way). Let the wavelet be the nth derivative of a smooth function (e.g. a Gaussian), ψ = dxn φ, with n α + 1. Then the CWT of s(x) reads xo − b S(b, a) ∼ aα dxn−α−1 φ . a Assume now that the modulus of dxn−α−1 φ has N maxima {φm , m = 1, . . . , N } at positions {xm , m = 1, . . . , N }. Then, for each a, the mod- ulus |S(b, a)| has N maxima localized at positions {bm = axm + xo , m = 1, . . . , N }, which converge toward xo as a → 0. Furthermore, the maxima of |S(b, a)| lie on N lines, called vertical ridges {b = axm + xo , m = 1, . . . , N }, which converge toward the singularity xo of the signal, and the modulus of |S(b, a)| along the mth ridge behaves as aα : |S(b = axm + xo , a)| = Γ(α + 1) aα φm . (2.18) March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 12 250 (a) 200 Force [Volts] 150 100 50 500 600 700 800 900 1000 1100 1200 1300 1400 1500 Time [10−4 s] (b) (c) 1.86 1.86 3.59 3.59 6.90 6.90 Scale Scale 13.29 13.29 25.58 25.58 49.25 49.25 600 700 800 900 1000 1100 1200 1300 1400 1500 600 700 800 900 1000 1100 1200 1300 1400 1500 −4 −4 Time [10 s] Time [10 s] Figure 7. Analysis of the rebound signal, with a Mexican hat wavelet: (a) The signal and the points detected by the respective ridges; (b) The modulus of the CWT; and (c) The corresponding skeleton (from Ref. 12). Hence the strength α of the singularity may be read oﬀ a log-log plot: ln |S(axm + xo , a)| ∼ α ln a + ln φm . (2.19) o This is, of course, the analysis of the local H¨lder regularity of the signal, and it is based on the covariance property (2.10) of the CWT under dilation. We illustrate this situation with two examples. The ﬁrst one (Figure 6) is an academic signal, with singularities of various strengths, analyzed with a Mexican hat wavelet. The skeleton of the CWT clearly shows the vertical ridges pointing towards each singularity. The second example comes from the analysis of the behavior of a ma- terial under impact, as made in Ref.12. The physical context is that of a so-called ‘instrumented falling weight impact’ testing. During such a test, a striker falls from a certain height on a clamped disk, so that either the striker rebounds or the disk breaks. In both cases, one records the time and the force on the striker. This type of event occurs on a very short time scale and is thus essentially transient, so that a time-frequency method March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 13 is required for the analysis. We focus here on the case of a rebound and a wavelet analysis of the force signal. The Mexican hat detects precisely three discontinuity points, namely, ﬁrst contact, maximal penetration and last contact. The results are shown in Figure 7. The signal is given in panel (a), whereas the CWT and its skeleton are presented in (b) and (c). The latter, in particular, shows the three vertical ridges that point towards the three instants mentioned. A further analysis exploits the behavior of the modulus of the CWT along each ridge. This yields precious insight into the physics of the phenomenon, particularly in the case of the rupture of the sample, not shown here. 2.6. The discrete WT (DWT) One of the successes of the WT was the discovery that it is possible to construct functions ψ for which {ψjk , j, k ∈ Z} is indeed an orthonormal basis of L2 (R), while keeping the good properties of wavelets, including space and frequency localization. In addition, this approach yields fast algorithms, and this is the key to the usefulness of wavelets in many appli- cations. The construction is based on two facts. First, almost all examples of orthonormal bases of wavelets can be derived from a multiresolution analy- sis, and then the whole construction may be transcripted into the language of Quadrature Mirror Filters (QMF), familiar in the signal processing lit- erature. A multiresolution analysis of L2 (R) is an increasing sequence of closed subspaces . . . ⊂ V−2 ⊂ V−1 ⊂ V0 ⊂ V1 ⊂ V2 ⊂ . . . with j ∈ Z Vj = {0} and j∈Z Vj dense in L2 (R), and such that (1) f (x) ∈ Vj ⇔ f (2x) ∈ Vj+1 ; (2) There exists a function φ ∈ V0 , called a scaling function, such that the family {φ(x − k), k ∈ Z} is an orthonormal basis of V0 . Combining conditions (1) and (2), one sees that {φjk (x) ≡ 2j/2 φ(2j x − k), k ∈ Z} is an orthonormal basis of Vj . The space Vj can be interpreted as an approximation space at resolution 2j . Deﬁning now Vj ⊕ Wj = Vj+1 , (2.20) March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 14 Signal 2 1 0 1 0 −1 1 0 −1 1 0 −1 2 0 −2 2 0 −2 2 0 −2 10 0 −10 100 200 300 400 500 600 700 800 900 1000 Figure 8. Five level decomposition of the academic signal on an orthonormal basis of Daubechies d6 wavelets. The low resolution approximation (∈ V−5 ) is shown on the bottom panel and the ﬁve levels of details with increasing resolution (∈ Wj , j = −5, . . . , −1) in the next ﬁve panels. we see that Wj contains the additional details needed to improve the reso- lution from 2j to 2j+1 . Thus one gets the decomposition L2 (R) = Wj . (2.21) j∈ Z Equivalently, one may choose a lowest resolution level jo and obtain ∞ L2 (R) = Vjo ⊕ Wj , (2.22) j=jo The crucial theorem then asserts the existence of a function ψ, some- times called the mother wavelet, explicitly computable from φ, such that {ψjk (x) ≡ 2j/2 ψ(2j x − k), j, k ∈ Z} constitutes an orthonormal basis of L2 (R): these are the orthonormal wavelets. Well-known examples include the Haar wavelets, the B-splines, and the various Daubechies wavelets. In practice, one starts from a sampled signal, taken in some VJ , and then the decomposition (2.22) is replaced by the ﬁnite representation J−1 VJ = Vjo ⊕ Wj . (2.23) j=jo March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 15 Figure 8 shows an example of a decomposition of order 5, namely the academic signal of Figure 6 decomposed over an orthonormal basis of Daubechies d6 wavelets.7 Thus we take J = 0 and jo = −5 in (2.23): V0 = V−5 ⊕ W−5 ⊕ W−4 ⊕ W−3 ⊕ W−2 ⊕ W−1 . (2.24) The construction of orthogonal wavelets then proceeds with (fast) al- gorithms for generating in a hierarchical way the various coeﬃcients in a decomposition like (2.24), exploiting standard procedures from signal pro- cessing. The resulting tool is quite eﬃcient, but in fact too rigid. Indeed, once the scaling function is given, everything is ﬁxed, in particular, one has no freedom left in the design of the wavelet ψ. In addition, the DWT is no longer covariant under (discrete) translations, the so-called shift invariance is lost, which is a seroius ﬂaw for pattern recognition. Therefore, vari- ous generalizations have been proposed in the literature. To name a few: biorthogonal wavelet bases, wavelet packets and the Best Basis Algorithm, the lifting scheme and second generation wavelets, the redundant WT (here one uses a rectangular lattice, instead of the dyadic one, so as to restore a (discrete) translation covariance). We shall refrain from describing these, for lack of space. Further information may be found in Ref.6. 2.7. Applications of the 1-D CWT The CWT has found a wide variety of applications in various branches of physics and/or signal processing. We will list here a representative selection of one-dimensional applications, in order to convey to the reader a feeling about the scope and richness of the ﬁeld. Most of the early applications, and the original references, may be found in the proceedings volumes Refs. 13–15. In all cases, the CWT is primarily used for analyzing transient phenomena, detecting abrupt changes in a signal or comparing it with a given pattern. Here is the list: • Sound and acoustics: musical synthesis, speech analysis (formant detection), disentangling of an underwater acoustic wavetrain. • Geophysics: analysis of microseisms in oil prospection, gravime- try (ﬂuctuations of the local gravitational ﬁeld), seismology, geo- magnetism (ﬂuctuations of the Earth magnetic ﬁeld), astronomy (ﬂuctuations of the length of the day, variations of solar activity, measured by the sunspots, etc). • Fractals, quasicrystals, turbulence (1-D and 2-D): diﬀusion limited aggregates, arborescent growth phenomena, structure of quasiperi- March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 16 odic patterns, identiﬁcation of coherent structures in developed turbulence. • Atomic physics: analysis of harmonic generation in laser-atom in- teraction. • Spectroscopy: in NMR spectroscopy, subtraction of spectral lines, noise ﬁltering. • Medical and biological applications: analyzing or monitoring of EEG, VEP, ECG; long-range correlations in DNA sequences. • o Analysis of local singularities: determination of local H¨lder expo- nents of functions. • Shape characterization: in robotic vision, contour of an object treated as a complex curve in the plane. • Industrial applications: monitoring of nuclear, electrical or me- chanical installations; quality assessment of telephone lines; analy- sis of behavior of materials under impact. 3. The two-dimensional CWT 3.1. Basic formulas Now we turn to the two-dimensional CWT, which has become a major tool in image processing. In this context, an image is a two-dimensional signal of ﬁnite energy, represented by a complex-valued function s ∈ L2 (R2 , d2 x): 2 s = d2 x |s(x)|2 < ∞. (3.1) R2 (sometimes it is useful to take s integrable as well). In practice, a black and white image will be represented by a bounded non-negative function: 0 s(x) M, ∀x ∈ R2 (M > 0), (3.2) the discrete values of s(x) corresponding to the level of gray of each pixel. Given an image s, all the geometric operations we want to apply to it are obtained by combining three elementary transformations of the plane, namely, rigid translations in the plane of the image, dilations or scaling (global zooming in and out) and rotations. Explicitly, the transformations act on x ∈ R2 in the familiar way: (i) translation by b ∈ R2 : x → x = x + b, (ii) dilation by a factor a > 0 : x → x = ax, (iii) rotation by an angle θ : x → x = rθ (x), March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 17 where cos θ − sin θ rθ ≡ ,0 θ < 2π, sin θ cos θ is the familiar 2 × 2 rotation matrix. The combined action of these three types of transformations is real- ized by the following unitary map in the space L2 (R2 , d2 x) of ﬁnite energy signals: U (b, a, θ)s (x) ≡ sb,a,θ (x) = a−1 s(a−1 r−θ (x − b)). (3.3) In terms of this action, the basic formulas for the 2-D CWT read S(b, a, θ) = a−1 d2 x ψ(a−1 r−θ (x − b)) s(x) (3.4) R2 =a d2 k eib·k ψ(ar−θ (k)) s(k) (3.5) R2 As in 1-D, we have to impose an admissibility condition on the wavelet ψ, namely, |ψ(k)|2 cψ ≡ (2π)2 d2 k < ∞, (3.6) R2 |k|2 which again may be replaced in practice by the following necessary condi- tion: ψ(0) = 0 ⇐⇒ d2 x ψ(x) = 0. (3.7) R2 The important observation to make here is that all the formulas are almost identical in 1-D and in 2-D ! As a consequence, the interpretation of the CWT as a singularity analyzer still holds, and the mathematical properties of the 2-D transform strictly parallel those of its 1-D counterpart. 3.2. Group-theoretical justiﬁcation However, before proceeding, we should pause and ask the reason of this situation. As can be expected, the answer lies in group theory. (1) In one dimension Observe that dilations and translations are aﬃne transformations of the line: x = (b, a)y ≡ ay + b, a > 0, b ∈ R, March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 18 and they obey the following composition rule: (b, a)(b , a ) = (b + ab , aa ). Thus, the set of all these aﬃne transformations constitutes a group, the (connected) aﬃne group of the line, {(b, a)} ≡ G+ . aﬀ The natural action of (b, a) on a signal, ψ → U (b, a)ψ, reads x−b (U (b, a)ψ)(x) = a−1/2 ψ (3.8) a and U is a unitary irreducible representation of G+ in L2 (R). In addition, aﬀ U is square integrable: db da ψ admissible ⇐⇒ | U (b, a)ψ|ψ |2 < ∞ (3.9) G + a2 aﬀ This explains the mathematical meaning of the admissibility condition of wavelets. (2) In two dimensions In two dimensions, dilations, translations, and rotations together con- stitute the similitude group of the plane: SIM(2) = R2 (R+ × SO(2)), ∗ with action x = (b, a, θ)y ≡ arθ y + b. This transformation is realized by the following unitary map on ﬁnite energy signals: U (b, a, θ)s (x) = a−1 s(a−1 r−θ (x − b)) (3.10) and U is a unitary irreducible representation of SIM(2) in L2 (R2 ). In addition, U is square integrable: da 2 ψ admissible ⇐⇒ d2 b dθ U (b, a, θ)ψ|ψ <∞ (3.11) SIM(2) a3 3.3. Interpretation of the 2-D CWT The 2-D formulas being the exact analogues of the ones in 1-D, in particular the admissibility condition (3.6) or (3.7), clearly the interpretation of the CWT will be exactly the same as in 1-D. Thus, combining the localization properties of ψ with the fact that the CWT is a convolution with a zero mean function, we conclude again that the 2-D CWT • yields a local ﬁltering, this time in all three variables b, a, θ, in particular, a directional ﬁltering; March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 19 • may be seen as a mathematical directional microscope with optics ψ, global magniﬁcation 1/a, and orientation tuning parameter θ; • works with constant relative bandwidth: ∆k/k = const, k = |k|, thus it is particularly eﬃcient for large spatial frequencies. Therefore we may say that the 2-D CWT is a detector and analyzer of singularities (edges, contours, corners, . . . ). 3.4. Main properties of the 2-D CWT For the same reason, the mathematical properties of the 2-D CWT are essentially the same as in the 1-D case, so we list them without further comment: (1) Energy conservation da c−1 ψ d2 b dθ|S(b, a, θ)|2 = d2 x |s(x)|2 . (3.12) a3 (2) Reconstruction formula da s(x) = c−1 ψ d2 b dθ ψb,a,θ (x) S(b, a, θ), (3.13) a3 i.e., we have a decomposition of the signal in terms of the analyzing wavelets ψb,a,θ , with coeﬃcients S(b, a, θ). (3) Reproduction property (reproducing kernel) da S(b , a , θ ) = c−1 ψ d2 b dθ ψb ,a ,θ |ψb,a,θ S(b, a, θ). (3.14) a3 (4) the CWT is covariant under translations, dilations and rotations, that is, under SIM(2). 3.5. Choice of the analyzing wavelet Even more than in 1-D, it is important to choose a wavelet that is well adapted to the problem at hand. One can distinguish two classes, isotropic wavelets and direction sensitive wavelets. (i) Isotropic wavelets If one decides to perform a pointwise analysis, or if directions are ir- relevant, it is suﬃcient to use a rotation invariant wavelet. Two standard examples are: March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 20 20 2 40 1.5 1 60 0.5 80 0 −0.5 5 −5 100 0 0 120 5 −5 20 40 60 80 100 120 20 0.8 40 0.6 60 0.4 0.2 80 0 5 −5 100 0 0 120 5 −5 20 40 60 80 100 120 Figure 9. An isotropic wavelet: The 2-D Mexican hat: (top) In position space; (bottom) In spatial frequency space. (1) The 2-D Mexican hat wavelet This wavelet (originally introduced by Marr in his pioneering work on vision16 ) is simply the Laplacian of a Gaussian (Figure 9): 2 2 ψH (x) = (2 − |x|2 ) e−|x| /2 , ψH (k) = |k|2 e−|k| /2 (3.15) (2) The Diﬀerence-of-Gaussians or DOG wavelet 2 /2α2 2 ψD (x) = 1 2α2 exp−|x| − exp−|x| /2 (0 < α < 1) (3.16) This is a very good substitute to the Mexican hat; for α−1 = 1.6, there are almost undistinguishable. (ii) Directional wavelets On the other hand, if the goal is to detect directional features in an im- age, or to perform directional ﬁltering, clearly one should resort to a direc- tion sensitive wavelet. The most eﬃcient solution is a directional wavelet, March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 21 20 0.08 15 0.06 10 0.04 0.02 5 0 −0.02 0 −0.04 −5 −0.06 −0.08 −10 0 10 −15 20 30 30 35 40 40 15 20 25 −20 0 5 10 −20 −15 −10 −5 0 5 10 15 20 Figure 10. The 2-D Morlet wavelet, for = 2, ko = (0, 6): (left) In position space (real part); (right) In spatial frequency space. that is, a wavelet ψ(x) such that the numerical support of its Fourier trans- form ψ(k) is contained in a convex cone with apex at the origin. Two useful examples are as follows: (1) The 2-D Morlet wavelet 2 ψM (x) = eiko ·x e|x| /2 + corr. (3.17) 2 ψM (k) = e|k−ko | /2 + corr. (3.18) As in 1-D, the correction term is required in order to satisfy the admissibility condition, and here too, it is negligible for |ko | 5.6. The 2-D Morlet wavelet is shown in Figure 10 for = 2, ko = (0, 6). (2) Conical wavelets, with support in the convex cone C(−α, α) ≡ {k ∈ R2 | − α arg k α, α < π/2} 1 2 (k · e−α )m (k · eα )m e− 2 kx , k ∈ C(−α, α) ˜ ˜ ψC (k) = (3.19) 0, otherwise ˜ where eφ is the unit vector in the direction φ and α = −α + π/2, so that e−α · eα = eα · e−α = 0. The frequency space version of this so-called ˜ ˜ Gaussian conical wavelet is shown in Figure 11, in the case m = 4, α = 10◦ . March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 22 Figure 11. The Gaussian conical wavelet, in frequency space, for m = 4, α = 10◦ . 3.6. Ridges in the 2-D CWT As in 1-D, the 2-D CWT is highly redundant, so that the full information can be retrieved from small subsets of the parameter space. Here too, one can either choose a discrete subset and get a frame, or focus on the regions where most of the energy is concentrated. Thus one considers the energy density of the CWT (written here for the isotropic case): E[s](b, a) ≡ |S(b, a)|2 , (3.20) and restricts the attention to the lines of local maxima of E[s](b, a), that is, the ridges of the CWT. It should be emphasized, however, that several deﬁnitions have been proposed in the literature. The important point is that, as in 1-D, the restriction of the CWT to the set of its ridges (the skeleton of the CWT) characterizes the signal completely. Here we will deﬁne a (vertical) ridge R as a 3-D curve (r(a), a) such that, for each scale a ∈ R+ , E[s](r(a), a) is locally maximum in space and r is a continuous function of scale. Figure 12 gives a concrete example. The signal (left panel) represents a set of singularities on a smooth background, simulating a set of bright points on the surface of the Sun and modelled by a random distribution of Gaussians of small (but random) width. The corresponding vertical ridges of the CWT of that signal are shown on the right panel, clearly each ridge points towards a singularity. Given such a vertical ridge, one may distinguish three characteristic features: March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 23 Figure 12. An example of a 2-D ridge: (left) The signal: a ﬁeld of singularities, simu- lating a set of bright points on the surface of the Sun; (right) The corresponding vertical ridges of the CWT of that signal. . The amplitude of the ridge AR = lim E[s](r(a), a). a→0 . The slope of E[s] on the ridge d ln E[s](r(a), a) SR = lim . a→0 d ln a . The energy of the ridge amax da ER = E[s](r(a), a), 0 a3 where [0, amax ] is the interval of deﬁnition of r(a). All three features may be used for extracting information on the underlying physical process. The ﬁrst two, in particular, have been exploited for discriminating between sev- eral classes of small features (bright points, cosmic impacts) in images of the Sun.17 3.7. Applications of the 2-D CWT As in 1-D, even more so, the 2-D CWT has been applied to a considerable number of problems. It is useful to distinguish between two types of appli- cations, namely those pertaining to image processing and those belonging to genuine physical problems. Detailed information and references may be found in Ref. 10. In addition, two novel applications, one for each type, are presented in separate paper in this volume.18 March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 24 1 1 0.5 0.5 0 0 Z Z −0.5 −0.5 −1 −1 −1 −1 −1 −1 −0.5 −0.5 −0.5 −0.5 0 0 0 0 0.5 0.5 0.5 0.5 1 1 1 1 X X Y Y Figure 13. Two spherical wavelets, in diﬀerent positions : (top) The DOG wavelet; (bottom) The real part of the spherical Morlet wavelet at scale a = 0.03. (a) Applications in image processing • Contour detection, character recognition: detection of edges, con- tours, corners . . . • Object detection and recognition in noisy images: automatic target recognition (ATR), application to infrared radar imagery, using both position and scale-angle features. • Image retrieval: recognition of a particular image in a large data bank, characterization of images by particular features. • Medical imaging: Magnetic resonance imaging (MRI), contrast en- hancement, segmentation. • Detection of symmetries in 2-D patterns: detection of discrete inﬂa- March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 25 tion (rotation + dilation) symmetries, quasicrystals (mathematical and genuine), quasiperiodic point sets. • Image denoising: removal of noise in images using directional wavelets. • Watermarking of images: adding a robust, but invisible, signature in images (e.g. with directional wavelets). (b) Physical applications • Astronomy and astrophysics: : structure of the Universe, cosmic microwave background (CMB) radiation, feature detection in im- ages of the Sun, detection of gamma-ray sources in the Universe. • Geophysics: fault detection in geology, detection of arrival time of individual seismic waves in seismology, hurricane detection in climatology. • Fluid dynamics: detection of coherent structures in turbulent ﬂu- ids, measurement of a velocity ﬁeld, disentangling of an underwater acoustic wave train. • Fractals and the thermodynamical formalism: analysis of 2-D frac- tals by the WTMM method, determination of fractal dimension, unraveling of universal laws, shape recognition and classiﬁcation of patterns. • Texture analysis: classiﬁcation of textures, “Shape from texture” problem. 3.8. Generalizations The transition from one to two dimensions gives us a hint as how to gener- alize further the CWT: given a class of (ﬁnite energy) signals, one speciﬁes the type of transformations one wants to apply to these signals. Then, if the latter form a group and if the action of that group on the signals is a uni- tary, square integrable representation, the formalism applies exactly in the same way (square integrability may be relaxed somewhat) and an adapted wavelet transform may be constructed (the formalism in question stems from the general construction of coherent states in quantum physics19 ). In this section, we will list very brieﬂy a number of such extensions. (1) Three-dimensional wavelets This is the most immediate extension, since all formulas of the CWT extend with obvious adaptations to three or more dimensions. The new March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 26 aspect here is that axisymmetric wavelets play a prominent role, in which the parameter space is now X = R3 (R+ × S 2 ) ∼ R3 × R3 , thus we have ∗ ∗ again a doubling of dimensions (and X may be interpreted as a phase space in the sense of Hamiltonian mechanics). The same situation holds true in any dimension d 3. (2) Wavelets on the 2-sphere The CWT on the 2-sphere S 2 may be obtained, using group-theoretical techniques, from the usual 2-D CWT on the tangent plane at the north Pole. The precise link is the inverse stereographic projection θ S2 (θ, ϕ) ⇐⇒ (r, ϕ) ≡ (2 tan , ϕ) ∈ R2 . 2 In particular, all the 2-D wavelets may be lifted to the 2-sphere. We present in Figure 13 two examples of spherical wavelets obtained in this way, the spherical DOG wavelet and the spherical Morlet wavelet, respectively. Clearly the spherical CWT will ﬁnd applications in geography (maps) and in astrophysics, whenever the of the Earth, resp. the universe, has to be taken into account. An example of the latter type is the analysis of the Cosmic Microwave Background (CMB) radiation (see Refs. 10, 17 and 20 for further details and references). (3) Spatio-temporal wavelets and motion estimation Wavelets that depend on space and time can be designed starting from separate dilations and translations in both variables. Then, by a simple change of variables, the two dilations are replaced with two new opera- tions: global dilation (space-time zooming) and speed tuning; the former will catch the global size of objects, whereas the latter will detect their velocity. In more than one dimension, one simply adds the usual rotation parameter. In this way one obtains a tool for analyzing moving objects (video sequences). A systematic use of the partial energy densities then yields a tracking algorithm that outperforms many standard ones; in par- ticular, it does not lose track of the target when the latter changes its shape, like a maneuvering airplane.10 4. The 2-D discrete wavelet transform The basic recipe of the 2-D discrete wavelet transform is simply to take the tensor product of two 1-D transforms: 2-D = 1-D ⊗ 1-D. Thus, by deﬁnition, one works in Cartesian coordinates. March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 27 c3 v d3 v d2 h d d3 d3 v d1 h d d2 d2 h d d1 d1 Figure 14. Schematic three level decomposition of an image into a low resolution ap- proximation, plus increasingly ﬁner details, of the 3 types (h, v, d). The ﬁrst step is to design a multiresolution analysis in 2-D along these lines, that is, from a 1-D multiresolution analysis {Vj , j ∈ Z} of L2 (R), one (2) builds a 2-D multiresolution analysis { Vj = Vj ⊗ Vj , j ∈ Z} of L2 (R2 ). The net result is that one needs now . one scaling function Φ(x, y) = φ(x) φ(y); . three wavelets Ψh (x, y) = φ(x) ψ(y) : detects horizontal features v Ψ (x, y) = ψ(x) φ(y) : detects vertical features d Ψ (x, y) = ψ(x) ψ(y) : detects oblique features. Using this approach, one is led to the standard presentation of the decomposition of images schematized in Figure 14. A concrete example, using the familiar image lena, is shown in Figure 15. Obviously, despite its widespread use, this technique has a number of drawbacks. On one hand, it oﬀers a poor sensitivity to directions, since one is bound to the Cartesian coordinates (x, y). This is natural in certain cases, like TV, but often it is a nuisance: even displaying isotropic features is nontrivial! Worse, the discrete scheme lacks “shift invariance” (properly, translation covariance), even a discrete one, when one uses the dyadic grid. This is often disastrous, in particular in the context of form recognition. March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 28 Figure 15. Typical three level decomposition of an image into a low resolution approx- imation, plus increasingly ﬁner details, of the 3 types (h, v, d). Thus, quite naturally, a number of generalizations have been proposed in the literature, exactly as in 1-D: . Redundant transforms, using a Cartesian lattice . Biorthogonal wavelet bases . Wavelet packets . The lifting scheme : second generation wavelets. Actually there is a worse aspect. Indeed some features in an image are not really two-dimensional. a curve is in fact a one-dimensional structure, albeit not straight. The standard 2-D DWT described above completely misses that aspect, thus it is a terrible waste of computing time. This is simply visible by estimating the number of squares it takes to cover a portion of March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 29 a curve at resolution 2−j , for increasing j. This phenomenon, dubbed the “curse of dimension”, may be avoided by deﬁning various new transforms, better adapted to geometry : ridgelets, curvelets, complex wavelets, non- linear approximations. There has been a rich supply in proposals in the recent years. A look at the most recent proceedings volume is eloquent.21 References 1. E. P. Wigner, On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40, 749–759 (1932). e a 2. J. Ville, Th´orie et applications de la notion de signal analytique, Cˆbles et Transm. 2`me A, 61–74 (1948). e 3. O. Rioul and M. Vetterli, Wavelets and signal processing, IEEE SP Magazine, October 1991, 14–38. 4. C. Heil and D. Walnut, Continuous and discrete wavelet transforms, SIAM Review 31, 628–666 (1989). 5. B. B. Hubbard, The World According to Wavelets, 2nd ed. (A.K. Peters, Wellesley, MA, 1998). 6. Y. Meyer, Wavelets: Algorithms and Applications (SIAM, Philadelphia, PA, 1993). 7. I. Daubechies, Ten Lectures on Wavelets (SIAM, Philadelphia, PA, 1992). 8. S. G. Mallat, A Wavelet Tour of Signal Processing, 2nd ed. (Academic Press, San Diego, 1999). 9. J. C. van den Berg (ed.), Wavelets in Physics (Cambridge Univ. Press, Cam- bridge (UK), 1999). 10. J-P. Antoine, R. Murenzi, P. Vandergheynst, and S. T. Ali, Two-dimensional Wavelets and Their Relatives (Cambridge University Press, Cambridge (UK), 2004, in press). 11. J-P. Antoine, D. Barache, R. M. Cesar, Jr, and L. da F. Costa, Shape char- acterization with the wavelet transform, Signal Process., 62, 265–290 (1997). 12. P. Vandergheynst, E. Van Vyve, A. Goldberg, J-P. Antoine, and I. Doghri, Modeling and simulation of an impact test using wavelets, analytical solutions and ﬁnite elements, International J. of Solids and Structures 38, 5481–5508 (2001). 13. J-M. Combes, A. Grossmann and Ph. Tchamitchian (eds.), Wavelets, Time- Frequency Methods and Phase Space (Proc. Marseille 1987) (Springer, Berlin, 1989; 2d Ed. 1990). 14. Y. Meyer (ed.), Wavelets and Applications (Proc. Marseille 1989) (Springer, Berlin, and Masson, Paris, 1991). 15. Y. Meyer and S. Roques (eds.), Progress in Wavelet Analysis and Applica- e tions (Proc. Toulouse 1992) (Ed. Fronti`res, Gif-sur-Yvette, 1993). 16. D. Marr, Vision (Freeman, San Francisco, 1982). March 18, 2004 8:53 WSPC/Trim Size: 9in x 6in for Proceedings jpa˙wavelets 30 17. J-P. Antoine, L. Demanet, J-F. Hochedez, L. Jacques, R. Terrier, and E. Verwichte, Application of the 2-D wavelet transform to astrophysical images, Physicalia Magazine, 24, 93–116 (2002). 18. J-P. Antoine and L. Jacques, The 2-D wavelet transform in image processing: Two novel applications, in this volume. 19. S. T. Ali, J-P. Antoine and J-P. Gazeau, Coherent States, Wavelets and Their Generalizations (Springer, New York, 2000). o u 20. L. Cay`n, J. L. Sanz, E. Martinez-Gonzalez, A. J. Banday, F. Arg¨eso, J. E. Gallegos, K. M. Gorski, and G. Hinshaw, Spherical Mexican Hat wavelet: An application to detect non-Gaussianity in the COBE-DMR maps, Mon. Not. R. Astron. Soc., 326, 1243–1249 (2001). 21. M. A. Unser, A. Aldroubi, and A. F. Laine (eds.), Wavelet Applications in Signal and Image Processing X , Proc. SPIE, vol. 5207 (2003).