Monday, February 9, 2009

The discrete wavelet transform

Introduction

For processing classical images the sampling is made in accordance with Shannon's [32] well-known theorem. The discrete wavelet transform (DWT) can be derived from this theorem if we process a signal which has a cut-off frequency. For such images the frequency band is always limited by the size of the camera aperture.

A digital analysis is provided by the discretisation of formula 14.1, with some simple considerations on the modification of the wavelet pattern by dilation. Usually the wavelet function $\psi^*(x)$ has no cut-off frequency and it is necessary to suppress the values outside the frequency band in order to avoid aliasing effects. We can work in Fourier space, computing the transform scale by scale. The number of elements for a scale can be reduced, if the frequency bandwidth is also reduced. This is possible only for a wavelet which also has a cut-off frequency. The decomposition proposed by Littlewood and Paley [22] provides a very nice illustration of the reduction of elements scale by scale. This decomposition is based on an iterative dichotomy of the frequency band. The associated wavelet is well localized in Fourier space where it allows a reasonable analysis to be made although not in the original space. The search for a discrete transform which is well localized in both spaces leads to multiresolution analysis.

Examples of Wavelets

Morlet's Wavelet

The wavelet defined by Morlet is [16]:
$\displaystyle \hat{g}(\omega) = e^{ -2 \pi^2(\nu - \nu_0)^2}$ (14.7)

it is a complex wavelet which can be decomposed in two parts, one for the real part, and the other for the imaginary part.
\begin{eqnarray*}g_r(x) & = & \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}} \cos(2\p... ... = & \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}} \sin(2\pi\nu_0 x) \end{eqnarray*}

where $\nu_0$ is a constant. The admissibility condition is verified only if 0.8$" align="middle" border="0" height="39" width="82">. Figure 14.1 shows these two functions.
Figure 14.1: Morlet's wavelet: real part at left and imaginary part at right.
\begin{figure} \centerline{ \hbox{\psfig{figure=fig_morlet.ps,bbllx=1cm,bblly=13.5cm,bburx=20.5cm,bbury=27cm,height=5cm,width=15cm,clip=} }} \end{figure}

Mexican Hat

The Mexican hat defined by Murenzi [30] is:
$\displaystyle g(x) = (1 - x^2) e^{-\frac{1}{2}x^2}$ (14.8)

it is the second derivative of a Gaussian (see figure 14.2).
Figure 14.2: Mexican Hat
\begin{figure} \centerline{ \hbox{\psfig{figure=fig_mexicain.ps,bbllx=1cm,bblly=13.5cm,bburx=10.5cm,bbury=27cm,height=5cm,width=10cm,clip=} }} \end{figure}

The continuous wavelet transform

The Morlet-Grossmann definition of the continuous wavelet transform [17] for a 1D signal $f(x)\in L^2(R)$ is:

$\displaystyle W(a,b)=\frac{1}{\sqrt a}\int_{-\infty}^{+\infty}f(x) \psi^*(\frac{x-b}{a}) dx$ (14.1)

where z* denotes the complex conjugate of z, $\psi^*(x)$ is the analyzing wavelet, a (>0) is the scale parameter and b is the position parameter. The transform is characterized by the following three properties:
1.
it is a linear transformation,
2.
it is covariant under translations:
$\displaystyle f(x) \longrightarrow f(x-u) \qquad W(a,b)\longrightarrow W(a,b-u)$ (14.2)

3.
it is covariant under dilations:
$\displaystyle f(x) \longrightarrow f(sx) \qquad W(a,b)\longrightarrow s^{-\frac{1} {2}}W(sa,sb)$ (14.3)

The last property makes the wavelet transform very suitable for analyzing hierarchical structures. It is like a mathematical microscope with properties that do not depend on the magnification.

In Fourier space, we have:

$\displaystyle \hat W(a,\nu)=\sqrt a \hat f(\nu)\hat{\psi}^*(a\nu)$ (14.4)

When the scale a varies, the filter $\hat{\psi}^*(a\nu)$ is only reduced or dilated while keeping the same pattern.

Now consider a function W(a,b) which is the wavelet transform of a given function f(x). It has been shown [#grossmann<#14252,#holschn<#14253] that f(x) can be restored using the formula:

$\displaystyle f(x)=\frac{1}{C_{\chi}} \int_0^{+\infty}\int_{-\infty}^{+\infty} \frac{1}{\sqrt a}W(a,b)\chi(\frac{x-b}{a})\frac{da.db}{a^2}$ (14.5)

where:
$\displaystyle C_{\chi}=\int_0^{+\infty} \frac{\hat{\psi}^*(\nu)\hat{\chi}(\nu)}{\nu}d\nu =\int_{-\infty}^0 \frac{\hat{\psi}^*(\nu)\hat{\chi}(\nu)}{\nu}d\nu$ (14.6)

Generally $\chi(x)=\psi(x)$, but other choices can enhance certain features for some applications.

The reconstruction is only available if $C_{\chi}$ is defined (admissibility condition). In the case of $\chi(x)=\psi(x)$, this condition implies $\hat \psi(0)=0$, i.e. the mean of the wavelet function is 0.

Endnote

Most of basic wavelet theory has been done. The mathematics have been worked out in excruciating detail and wavelet theory is now in the refinement stage. The refinement stage involves generalizations and extensions of wavelets, such as extending wavelet packet techniques.

The future of wavelets lies in the as-yet uncharted territory of applications. Wavelet techniques have not been thoroughly worked out in applications such as practical data analysis, where for example discretely sampled time-series data might need to be analyzed. Such applications offer exciting avenues for exploration.

Musical Tones

Victor Wickerhauser has suggested that wavelet packets could be useful in sound synthesis (13). His idea is that a single wavelet packet generator could replace a large number of oscillators. Through experimentation, a musician could determine combinations of wave packets that produce especially interesting sounds.

Wickerhauser feels that sound synthesis is a natural use of wavelets. Say one wishes to approximate the sound of a musical instrument. A sample of the notes produced by the instrument could be decomposed into its wavelet packet coefficients. Reproducing the note would then require reloading those coefficients into a wavelet packet generator and playing back the result. Transient characteristics such as attack and decay- roughly, the intensity variations of how the sound starts and ends- could be controlled separately (for example, with envelope generators), or by using longer wave packets and encoding those properties as well into each note. Any of these processes could be controlled in real time, for example, by a keyboard.

Notice that the musical instrument could just as well be a human voice, and the notes words or phonemes.

A wavelet-packet-based music synthesizer could store many complex sounds efficiently because

  • wavelet packet coefficients, like wavelet coefficients, are mostly very small for digital samples of smooth signals; and
  • discarding coefficients below a predetermined cutoff introduces only small errors when we are compressing the data for smooth signals.

Similarly, a wave packet-based speech synthesizer could be used to reconstruct highly compressed speech signals. Figure 8 illustrates a wavelet musical tone or toneburst.

Fig8

Fig. 8. Wavelets for music: a graphical representation of a Wickerhauser toneburst. This screenshot of a toneburst was taken while it was playing in the Macintosh commercial sound program Kaboom! Factory. (Toneburst courtesy Victor Wickerhauser)

Press Speaker Button to hear this wavelet toneburst (110 kB).

Detecting Self-similarity in a Time Series

Wavelet analysis is proving to be a very powerful tool for characterizing behavior, especially self-similar behavior, over a wide range of time scales.

In 1993, Scargle and colleagues at NASA-Ames Research Center and elsewhere investigated the quasiperiodic oscillations (QPOs) and very low-frequency noise (VLFN) from an astronomical X-ray accretion source, Sco X-1 as possibly being caused by the same physical phenomenon (12). Sco X-1 is part of a close binary star system in which one member is a late main sequence star and the other member (Sco X-1) is a compact star generating bright X rays. The causes for QPOs in X-ray sources have been actively investigated in the past, but other aperiodic phenomena such as VLFNs have not been similarly linked in the models. Their Sco X-1 data set was an interesting 5-20 keV EXOSAT satellite time-series consisting of a wide-range of time scales, from 2 ms to almost 10 hours.

Galactic X-ray sources are often caused by the accretion of gas from one star to another in a binary star system. The accreted object is usually a compact star such as a white dwarf, neutron star, or black hole. Gas from the less massive star flows to the other star via an accretion disk (that is, a disk of matter around the compact star flowing inward) around the compact star. The variable luminosities are caused by irregularities in the gas flow. The details of the gas flow are not well-known.

The researchers noticed that the luminosity of Sco X-1 varied in a self-similar manner, that is, the statistical character of the luminosities examined at different time resolutions remained the same. Since one of the great strengths of wavelets is that they can process information effectively at different scales, Scargle used a wavelet tool called a scalegram to investigate the time-series.

Scargle defines a scalegram of a time series as the average of the squares of the wavelet coefficients at a given scale. Plotted as a function of scale, it depicts much of the same information as does the Fourier power spectrum plotted as a function of frequency. Implementing the scalegram involves summing the product of the data with a wavelet function, while implementing the Fourier power spectrum involves summing the data with a sine or cosine function. The formulation of the scalegram makes it a more convenient tool than the Fourier transform because certain relationships between the different time scales become easier to see and correct, such as seeing and correcting for photon noise.

The scalegram for the time-series clearly showed the QPOs and the VLFNs, and the investigators were able to calculate a power-law to the frequencies. Subsequent simulations suggested that the cause of Sco-X1's luminosity fluctuations may be due to a chaotic accretion flow.

Denoising Noisy Data

In diverse fields from planetary science to molecular spectroscopy, scientists are faced with the problem of recovering a true signal from incomplete, indirect or noisy data. Can wavelets help solve this problem? The answer is certainly "yes," through a technique called wavelet shrinkage and thresholding methods, that David Donoho has worked on for several years (10).

The technique works in the following way. When you decompose a data set using wavelets, you use filters that act as averaging filters and others that produce details (11). Some of the resulting wavelet coefficients correspond to details in the data set. If the details are small, they might be omitted without substantially affecting the main features of the data set. The idea of thresholding, then, is to set to zero all coefficients that are less than a particular threshold. These coefficients are used in an inverse wavelet transformation to reconstruct the data set. Figure 6 is a pair of "before" and "after" illustrations of a nuclear magnetic resonance (NMR) signal. The signal is transformed, thresholded and inverse-transformed. The technique is a significant step forward in handling noisy data because the denoising is carried out without smoothing out the sharp structures. The result is cleaned-up signal that still shows important details.

Fig6

Fig. 6. "Before" and "after" illustrations of a nuclear magnetic resonance signal. The original signal is at the top, the denoised signal at the bottom. (Images courtesy David Donoho, Stanford University, NMR data courtesy Adrian Maudsley, VA Medical Center, San Francisco).

Figure 7 displays an image created by Donoho of Ingrid Daubechies (an active researcher in wavelet analysis and the inventor of smooth orthonormal wavelets of compact support), and then several close-up images of her eye: an original, an image with noise added, and finally denoised image. To denoise the image, Donoho:

  1. transformed the image to the wavelet domain using Coiflets with three vanishing moments,
  2. applied a threshold at two standard deviations, and
  3. inverse-transformed the image to the signal domain.

Fig7

Fig. 7. Denoising an image of Ingrid Daubechies' left eye. The top left image is the original. At top right is a close-up image of her left eye. At bottom left is a close-up image with noise added. At bottom right is a close-up image, denoised. The photograph of Daubechies was taken at the 1993 AMS winter meetings with a Canon XapShot video still-frame camera. (Courtesy David Donoho)