Wavelet Explorer Demonstration

Given Tuesday, October 28, 1997 
P.J. Hinton
Wolfram Research, Inc.

Initialization

Needs["Wavelets`Wavelets`"]

$TextStyle={FontFamily->"Times-Roman",FontSize->16};

Off[General::spell,General::spell1]

Introduction

General Definition

A wavelet ψ(x) is a small wave which satisfies an admissibility condition.

This requires that the square modulus of the wavelet's continuous Fourier transform must decay faster than (1/ω) as  x ->± ∞.  You can think of this as a measure of of a signal's band-limitedness.

Wavelet Transforms

Wavelet functions are useful as basis functions for transforms.  Wavelet bases are built by translating and dilating a parent function called the mother wavelet.  Two flavors of the wavelet transform exist:

Continuous Wavelet Transform

The continuous wavelet transform is defined over two real parameters a and b.

In a sense, a plays the role of the freqency parameter (ω) of the continuous Fourier transform.

The parameter b is an indicator of the region over which  (Subscript[ψ, a,b])  =   ψ(((x-b)/a))  is spatially localized.

Discrete Wavelet Transform

The discrete wavelet transform is obtained by restricting values of a to powers of two.  This transform is used widely in numerical calculations.

The underlying principle of the transform is to begin with a projection of a function onto some basis of spatially localized functions based on the dilation and translation of a scaling function.

Repeated linear transformations are applied to the original projection, resulting in a series coarsened renditions of the signal. The orthogonal complement to each projection (i.e. information lost) is an expansion of wavelet basis functions.

where (Subscript[ψ, j k ])(x) ≡ (2^(j/2))ψ((2^j)x - k).

Families of Wavelet Bases

The definition of a wavelet sufficiently broad to allow several different unrelated approaches for creating wavelet basis functions.  Some of the major ones include the following:

({{Wavelet Basis, Advantage}, {Haar, good tutorial function}, {Daubechies, compact support}, {Biorthogonal, symmetry}, {Shannon, bandlimited basis functions}, {Meyer, good smoothness properties}, {Battle Lemarié, closed form expressions}, {Local Trigonometric, highly adaptable}})

Areas of Application

Signal processing

Image compression

Acoustics

Fluid dynamics

What do you get with Wavelet Explorer ?

Documentation that provides an entry-level treatment of wavelet analysis

Several examples that are integrated with the tutorial material

A Mathematica package that contains several tools to get you going right away

Filters for defining a wide variety of common wavelet functions

Functions for to ease visualization of wavelet functions

Functions for computing derivatives of basis functions, where they exist

Wavelet transform and wavelet packet transform functions

Transforms for local trigonometric bases

Visualization functions for interpreting transform results

A Simple Comparison with Fourier Transforms

Here is a spatially localized, oscillatory function:

data = Table[{(x/256),(E^(-100((((x/256)-(3/8)))^2)))Sin[2 π 10 (x/256)]+ (E^(-200((((x/256)-(5/8)))^2)))Sin[2 π 20 (x/256)]},{x,0,255}];

ListPlot[data,PlotJoined ->True,PlotRange->All,ImageSize-> {400,250}];

Fourier Analysis

Here is an example of how one would perform a Fourier-style analysis on the data:

fourtrans = Fourier[(Transpose@ data)[[2]]];

modulus = Abs[fourtrans];

ListPlot[modulus,PlotRange->{{0,30},All},PlotJoined->True,ImageSize->{400,250}];

The modulus of the Fourier transform tells us that the function resembles a 10 Hz sinsoid somewhere, but we can't be sure exactly where.

Wavelet Analysis

Here is a wavelet analysis of the same data.

Most wavelets are referenced by a sequence of coefficients known as a "filter".

filter = DaubechiesFilter[7]

Let's take a look at the mother wavelet.

ListPlot[Wavelet[filter,8],PlotJoined->True,PlotRange -> All, ImageSize-> {400,250}];

Computing a wavelet transform requires only the sampled data and the filter description of the wavelet basis.

transform = WaveletTransform[(Transpose@ data)[[2]],filter];

Wavelet Explorer comes with tools that allow you to visualize the transform.  Here are two common methods of visualization.

PlotCoefficients[transform,Background ->RGBColor[0, 0.3, 0.6],ImageSize->{400,250},Axes->{True,True},Ticks->{ Table[{x,x/128},{x,0,255,256/8}],None},AxesLabel -> "Scale ((2^(-j)))"];

Here, we see the magnitudes of the wavelet coefficients represented as vertical bars.   The horizontal lines represent groups of wavelet basis functions at increasing levels of resolution.

This graph gives us not only an idea of the frequency content of the waveform, but it also tells us when each frequency is dominant.

PhaseSpacePlot[transform,Frame->True,FrameTicks -> None,FrameLabel->{"time","frequency"} ,LogarithmicScale-> True ,ImageSize->{400,250}];

This shows how wavelet bases carve up the time-frequency plane into tiles.  Magnitudes of wavelet coefficients are shown with grayscales, with darker values denoting larger magnitudes.

This graphic underscores the tradeoff between knowledge of time and frequency information about a function.  With wavelets we gain information about time at the expense of frequency knowledge.

Denoising a Signal

Wavelets can be used to filter out noise from signal.

data = Get["Wavelets/Data/shocknoi.dat"];

ListPlot[data, PlotJoined -> True ,AxesLabel -> {"x", "u(x)"},ImageSize-> {400,250}];

This generates a definition for a Battle-Lemarié spline wavelet.

s4 = SplineFilter[4, 8];

Here we compute the wavelet transform of the noisy signal.

wtdata = WaveletTransform[data, s4];

This performs a compression step.

Compress[Rest[wtdata], 1.5, Shrinking -> True];

The compression algorithm works according to the following rule:

Which zeroes out all wavelet coefficients less than a pre-specified threshold and expresses all other coefficients as a deviation from the threshold.

denoised = InverseWaveletTransform[
Join[{wtdata[[1]]}, %], s4];

ListPlot[denoised, PlotJoined -> True,
AxesLabel -> {"x", "u(x)"},ImageSize->{400,250}];

Here, we see the effects of increasing the threshold's magnitude.

Compress[Rest[wtdata], 3, Shrinking -> True];

ListPlot[InverseWaveletTransform[
Join[{wtdata[[1]]}, %], s4], PlotJoined -> True,
AxesLabel -> {"x", "u(x)"},ImageSize -> {400,250}];

Image Compression

Here we show how local trigonometric bases can be used to compress an image (fingerprint data).

This reads in the data. It has dimension {240,192}

(data = Get["Wavelets/Data/finger.dat"];
Dimensions[data])

We display the data showing a finger print.

ListDensityPlot[data, Mesh -> False];

We perform local cosine transform using a rectangular window, that is, using Cos-IV basis.

Here we do cosine transform on the data with smooth windows.

ct = CosTransform[data, {48, 60}, Taper -> False]; // Timing // First

ct2 = CosTransform[data, {48, 60}, Taper ->
{{Polynomial[2], 24}, {Polynomial[2], 30}}]; // Timing // First

This gives the density plot of the coefficients.

PlotCoefficients2D[ct2];

We calculate the cumulative energy function for both transforms.

(norm = Flatten[data].Flatten[data];
 ctE = CumulativeEnergy[ct, 5000]/norm;
 ct2E = CumulativeEnergy[ct2, 5000]/norm;)

We display the error versus compression ratio curve for both transforms. The transform with smooth windows (solid line) is obviously superior to that with rectangular windows.

Show[MapThread[
ListPlot[Transpose[{Table[240 192/i, {i, 1000, 5000}],
Sqrt[1-Take[#, {1000, 5000}]]}], PlotStyle -> #2,
PlotJoined -> True, DisplayFunction -> Identity,
PlotLabel -> "relative error vs. compression ratio"]&,
{{ctE, ct2E}, {Dashing[{0.02}], {}}}],
DisplayFunction -> $DisplayFunction,ImageSize->{400,250}];

Keeping 10% of the total coefficients will give an error of 9%.

Sqrt[1-ct2E[[4608]]]

We do the inverse transform on the compressed data.

(compImage = InverseCosTransform[Compress[ct2, 4608], Taper -> {{Polynomial[2], 24}, {Polynomial[2], 30}}]) // Timing // First

This is the reconstructed image. Again it resembles the original very well.

ListDensityPlot[compImage, Mesh -> False];