Wavelets and Applications - Maple Programming Help

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : Applications and Example Worksheets : Discrete Mathematics : examples/Wavelets

Wavelets and Applications

Introduction

Wavelets are powerful tools that can be used in signal processing and data compression. Wavelet transforms are an excellent alternative to Fourier transforms in many situations. In Fourier analysis, a signal is decomposed into periodic components; in wavelet analysis, a signal is decomposed into components localized in both time and frequency domains. Thus, wavelet transforms are ideal when signals are not periodic.

The Theory

In Fourier analysis, sinn π x,cosn π x is used as an orthonormal basis of L20,1.  In wavelet analysis, a father wavelet φ and a mother wavelet ψ are chosen such that:

 

φx,ψ2 x, ψ2 x1,ψ4 x,...,ψ4 x3,ψ8 x,...

form an orthonormal basis of L20,1. In theory, φ is chosen to satisfy the conditions of a multiresolutional analysis (MRA), and then ψ  is determined from φ and the MRA. In practice, φ and ψ are assumed to satisfy the following functional equations, and the coefficients are computed according to the desired properties of the MRA.

φx=n=hnφ2xn    (1)

ψx=n=gnφ2xn    (2)

In fact, often ψ and φ cannot be determined symbolically, and are defined solely in terms of these coefficients. In such cases, the Cascades Algorithm can be used to obtain numerical approximations.

The fact that ψ and φ must be orthogonal reduces to the following numerical conditions on the hn and gn, when φ has norm 1.

n=hnhn+2h=δ0,k    (3)

n=gngn+2h=δ0,k    (4)

n=hngn+2h=0    (5)

Usually, the hn  are computed, then the gn are determined by reversing the hn and negating every term.  The hn and gn are also known as the scaling and wavelet coefficients, or the low pass and high pass filters, respectively.

Z Transforms

When the Fourier transform of equation (1) is computed,

Φw=n=∞∞hneinw2Φw2

is obtained. This is the motivation for defining:

m0x=n=∞∞hnxn2

This is known as the Z transform associated with φ. In this context, the orthogonality conditions reduce to:

 

m0xm0x+m0x+πm0xπ=2

The generation of wavelets is often phrased in this language.

 

Orthogonal and Biorthogonal Wavelets

Technically, the above discussion applies only to orthogonal wavelets. In a variant of this theory, different bases of L20,1 are used for the analysis and synthesis of signals. Such pairs of bases are generated by biorthogonal bases. Note that orthogonal wavelets can be viewed as biorthogonal wavelets for which the analysis and synthesis processes coincide. However, biorthogonal wavelets are not orthogonal; their main advantage is that the orthogonality conditions are relaxed, allowing more smoothness conditions to be imposed.

Some authors refer to wavelets that are neither orthogonal or biorthogonal, but such wavelets are not discussed here.

Wavelet Generation

This section provides examples showing how some wavelets are generated by Maple. This section is not intended as a complete guide to the generation of these wavelets, and it is not intended as a thorough discussion of their theory. It is a simplified outline of how Maple generates these wavelets.

Symlets

Symlets are a variant of the Daubechies wavelet. In fact, they are also called Daubechies least asymmetric wavelets. They have the same vanishing moments as the Daubechies wavelets, and the same size, but they have minimal phase. A complex valued function f is said to have linear phase if there are real numbers a and b such that fx=eI a x+b fx. It is known that there are no compactly supported (that is, finite length) orthogonal wavelets with linear phase. The Symlets were designed by Ingrid Daubechies to have phase as close as possible to linear.


The generation of Symlets is very similar to the generation of the Daubechies wavelets. To generate the Symlet of size 2 A, start by finding the roots of the polynomial P that is used in the generation of the Daubechies wavelet. Then transform these roots to get roots of the Laurent polynomial PZX, where ZX=2XX12. The plot demonstrates that these roots come in groups of conjugates and reciprocals, so you are able to restrict your attention to those with norm less than or equal to 1 and non zero imaginary part.

A7&colon;PaddbinomialA&plus;k1&comma;A12kXk&comma;k&equals;0..A1&semi;solsRootFinding:-AnalyticP&comma;X&comma;re&equals;100..100&comma;im&equals;100..100&colon;T:=mapt&rarr;op1t&plus;2t&plus;t212&comma;1t2t&plus;t212&comma;sols&colon;Tplotseq&real;Ti&comma;&Im;Ti&comma;i&equals;1..nopsT&colon;plotslistplotTplot&comma;style&equals;point&semi;Sselectz&rarr;z<1and0&Im;z&comma;T&colon;

P:=1&plus;72X&plus;7X2&plus;212X3&plus;1058X4&plus;23116X5&plus;23116X6

 

Now perform spectral factorization on PZX. This means that you want to find a pX such that PZX&equals;pX pX. This is done by using the Feyer Reiz algorithm. For each root λ of PX, you simply have to assign one of the roots of PZX&lambda; to px. Because of the properties of P, this is sufficient. This is where the Symlets are different from the Daubechies wavelets. In the generation of the Daubechies wavelets, you can simply pick the root in the unit circle of the complex plane. This choice of spectral factorization has maximal phase. To compute the choice of spectral factorization with minimal phase, you have to compute all spectral factorizations, and pick the one where the nonlinear part has the smallest L1 norm.  

phase procw&comma;ziftypez&comma;realconsthenarctan1&plus;ztanw21z  w2eliftypez&comma;complexconsthenarctan1 z2sinw1&plus;z2cosw2 zcosargumentzw&plus;piecewisecosw<2 zcosargumentz1&plus;z2&comma;π&comma;0end ifend proc&colon;

minphase &colon;minphaseat &colon;

 

Note that i and n should be thought of as binary numbers whose digits encode the factorization being used.

 

To save time, you can assign the first value. This means that you only range over half of all possible factorizations, but this is enough. Within the loop below, save a list of the phases, PhaseList, so that you can graph them later. Remember that minimal phase is defined by minimum L1 norm.

forifrom0to2nopsS11do   ni&semi;   &phi;phasew&comma;S1&semi;   for j from 2 to nopsSdo      &phi;&phi;&plus;2 modpn&comma;212phasew&comma;Sj&semi;      nnmodpn&comma;22&semi;     end do&semi;   PhaseListi&phi;&semi;   myIntevalfInt&phi;&comma;w&equals;0..&pi;&comma;method&equals;_Gquad&semi;  ifmyInt<minphasethen      minphasemyInt&semi;      minphaseati&semi;   end if&semi; end do&colon;

 

With this done, you can now graph all of the phases that were considered, with the minimal phase displayed with a bold line.

plotsdisplayplotseqPhaseListi&comma;i&equals;0..2nopsS11&comma;w&equals;0..&pi;&comma;plotPhaseListminphaseat&comma;w&equals;0..&pi;&comma;thickness&equals;5

 

Now construct the spectral factorization with the minimal factorization computed above, and extract the normalized scaling (low pass) coefficients.

n minphaseat&colon;if typeS1&comma;realconsthen   pxS1else   pxS1xS1&conjugate0; end if&colon; for j from 2 to nopsS do   if typeSj&comma;realcons then      pp xSj2modpn&comma;212   else      pp xSj2modpn&comma;212xSj&conjugate0;2modpn&comma;212   end if&semi;   nnmodpn&comma;22&semi; end do&colon;

pncollect&lpar;p1&plus;xA21Apx&equals;1|px&equals;1&comma;x&rpar;&colon; CListToolsReverseseq&real;coeffpn&comma;x&comma;k&comma;k&equals;0..2A1&colon;CaddCi2&comma;i&equals;1..nopsC 

0.002681814559&comma;0.001047384898&comma;0.01263630337&comma;0.03051551326&comma;0.06789269364&comma;0.04955283511&comma;0.01744125457&comma;0.5361019159&comma;0.7677643179&comma;0.2886296316&comma;0.1400472406&comma;0.1078082374&comma;0.004010244846&comma;0.01026817669

(3.1.1)

 

Verify this against the Maple function that computes wavelets.

l&comma;h DiscreteTransforms:-WaveletCoefficientsSymlet&comma;2A&colon;converth&comma;list

0.00268181456826015&comma;0.00104738488867974&comma;0.0126363034032406&comma;0.0305155131658779&comma;0.0678926935012206&comma;0.0495528349370428&comma;0.0174412550868357&comma;0.536101917090569&comma;0.767764317004883&comma;0.288629631750648&comma;0.140047240442934&comma;0.107808237703290&comma;0.00401024487152240&comma;0.0102681767084648

(3.1.2)

 

The Discrete Wavelet Transform

The wavelet transform can be accomplished for discrete signals by using an algorithm known as the (fast) discrete wavelet transform. Recall the coefficients hn and gn from equations (1) to (5). The low pass filter, w2, is the hn, and the high pass filter, w1, is the gn (in vector form). In almost all useful cases, these are finite. The size of w1 and w2 is called the filter length. If w1 has n elements, w1 is often called an n-tap wavelet.

 

The discrete wavelet transform produces two outputs, each half the size of the input. The first output is the high detail coefficients, and the second is the low detail coefficients. They are computed by convolving w1 and w2, respectively, against the input data, and then downsampling (throwing away every second term).

The low detail coefficients are then recursively processed. It is not obvious, but this in fact computes the wavelet coefficients (the coefficients of each function in the orthonormal base of wavelets).

The discrete wavelet transform is a linear transformation on the input signal. If the high and lows pass filter are:

 

w1&equals;w11&comma;w12&comma;w13&comma;w14

w2&equals;w21&comma;w22&comma;w23&comma;w24

 

then this linear transform on a signal of length 6 can be viewed as multiplication by the Matrix:

 

w11w12w13w1400w21w22w23w240000w11w12w13w1400w21w22w23w24w13w1400w11w12w23w2400w21w22

 

The result of the multiplication of this Matrix by the signal is an interlacing of the high and low pass coefficients.

The orthogonality conditions on wij (recall that w1 and w2 are the wavelet and scaling coefficients from equations (3) to (5)) are equivalent to the Matrix being orthogonal! This makes the discrete wavelet transform an orthogonal linear transformation, and hence very easy to invert.

End Conditions

In the above convolutions, the filter mask "falls off" the end of the data. To maintain the orthogonality of the discrete wavelet transform (and the resulting easy invertibility), the data must be assumed to be periodic (as shown in preceding diagrams). However, two other common alternatives exist. The data can be padded with zeros, or the data can be reflected to generate extra data. In both cases, the transform is usually not invertible with orthogonal or biorthogonal wavelets, but can sometimes be modified to maintain easy invertibility. Such modifications are not discussed here.

For example, given the signal [1,2,3,4,5,6], the following are the periodic, zeros, and reflection end conditions for extending the data.

Periodic:

1&comma;2&comma;3&comma;4&comma;5&comma;6&comma;1&comma;2&comma;3&comma;..&period;

Zeros:

1&comma;2&comma;3&comma;4&comma;5&comma;6&comma;0&comma;0&comma;0&comma;..&period;

Reflection:

1&comma;2&comma;3&comma;4&comma;5&comma;6&comma;5&comma;4&comma;3&comma;..&period;

 

Maple Functions for Wavelets

All of Maple's functions for wavelets are part of the SignalProcessing and DiscreteTransforms packages.

 

The SignalProcessing commands are:
- DWT
- InverseDWT

The DiscreteTransform package commands are:

- DiscreteWaveletTransform
- InverseDiscreteWaveletTransform

- WaveletCoefficients
- WaveletPlot

See the corresponding help pages for basic information and examples. Examples from the DiscreteTransforms package are described below. For more examples from the SignalProcessing package, see the SignalProcessing examples page.

Examples

This is a quick example to explore the Daubechies length 4 wavelet.

withDiscreteTransforms&semi;withImageTools&colon;

DiscreteWaveletTransform&comma;FourierTransform&comma;InverseDiscreteWaveletTransform&comma;InverseFourierTransform&comma;WaveletCoefficients&comma;WaveletPlot

(5.1.1)

ExW1&comma;ExW2 WaveletCoefficientsDaubechies&comma;4

ExW1&comma;ExW2:=0.1294095225512600.2241438680420130.8365163037378080.482962913144534&comma;0.4829629131445340.8365163037378080.2241438680420130.129409522551260

(5.1.2)

 

Plot the Daubechies mother and father wavelets by using the WaveletPlot command. WaveletPlot uses a numerical algorithm called Cascades to approximate and plot functions satisfying equations (1) and (2). No explicit definition exists of the Daubechies length 4 mother and father wavelets.

WaveletPlotExW1&comma;ExW2

 

The WaveletCoefficients command respects the Digits setting. In this case, if you increase the setting of Digits, you can correctly identify the symbolic expressions of the the wavelet coefficients.

Digits20&colon;ExW1&comma;ExW2WaveletCoefficientsDaubechies&comma;4&semi;mapidentify&comma;ExW1&semi;Digits10&colon;

ExW1&comma;ExW2:=0.12940952255126038117444941881202416417440.22414386804201338102597276224040035546800.83651630373780790557529378091687320345870.4829629131445341433748715998644486838167&comma;0.48296291314453414337487159986444868381670.83651630373780790557529378091687320345870.22414386804201338102597276224040035546800.1294095225512603811744494188120241641744

186&plus;182382&plus;186382&plus;186186182

(5.1.3)

 

And of course, you can use the Daubechies 4 wavelet to transform data.

ExV Vector10&comma;ii&comma;datatype&equals;float8&colon;ExTDiscreteWaveletTransformExV&comma;Daubechies&comma;4&semi;InverseDiscreteWaveletTransformExT&comma;Daubechies&comma;4

ExT:=2.2204460492503110-164.4408920985006310-164.4408920985006310-160.3.53553390593274&comma;2.310789034541155.139216159287347.9676432840335310.796070408779712.6771540786184

1.2.000000000000003.000000000000004.000000000000005.000000000000006.000000000000007.000000000000008.9.10.0000000000000

(5.1.4)

 

You can also transform an image, by using the ImageTools package.

PreImg Readcatkerneloptsdatadir&comma;/images/Lichtenstein.JPG&colon;ExImgMatrixToGrayscalePreImg&colon;ExImgHi&comma;ExImgLo DiscreteWaveletTransformExImg&comma;1&comma;Daubechies&comma;20&colon; ExImgHiHi&comma;ExImgHiLo&comma;ExImgLoHi&comma;ExImgLoLo DiscreteWaveletTransformExImgHi&comma;2&comma;Daubechies&comma;20&comma;DiscreteWaveletTransformExImgLo&comma;2&comma;Daubechies&comma;20&colon;Show FitIntensityExImgLoLo&verbar;FitIntensityExImgLoHi&comma;FitIntensityExImgHiLo&verbar;FitIntensityExImgHiHi&colon;

EmbedCreateShow

 

Sample Applications

Wavelets are of growing importance in a number of diverse fields, including seismology, underwater acoustics, computer vision, and signal processing. The largest application seems to be in image compression; below are three sample applications to illustrate of the capabilities of Maple's new wavelet functions.

Procedures and Initialization

First, define functions to transform (and inverse transform) a grayscale image. Discrete wavelet transforms of images are done by transforming first one dimension, and then the other. This process generates four outputs, the high-high, high-low, low-high, and low-low coefficients. The low-low coefficients can then be transformed recursively.

Also define a small procedure to count the zeros in a Matrix, returning the result as a percentage of the number of elements in the Matrix, and two procedures needed for signal denoising.

restart&semi;withDiscreteTransforms&colon;withImageTools&colon;

ImageDWT  procimg::OrMatrix&comma;Array&comma;w1::Vector&comma;w2::Vector&comma;iters::posint           local rows&comma;cols&comma;n&comma;final&comma;current&comma;temp&semi;   rows&comma;colsLinearAlgebra:-Dimensionsimg&semi;   n1&semi;   finalMatriximg&semi;   if modprows&comma;2iters<>0ormodpcols&comma;2iters<>0then      error img's dimensions are not divisible by %1&comma;2iters   end if&semi;   tempMatrixrtable_dimsfinal&comma;datatype&equals;float8&semi;   for n to iters do      DiscreteWaveletTransformrows2n1&comma;cols2n1&comma;final&comma;1&comma;w2&comma;w1&comma;temp&comma;storagetype&equals;singlearray&semi;      DiscreteWaveletTransformcols2n1&comma;rows2n1&comma;temp&comma;2&comma;w2&comma;w1&comma;final&comma;storagetype&equals;singlearray&semi;   end do&semi;   return final end proc&colon;

InverseImageDWT  procimg::OrMatrix&comma;Array&comma;w1::Vector&comma;w2::Vector&comma;iters::posint          local rows&comma;cols&comma;n&comma;final&comma;current&comma;temp&comma;temp2&comma;i&comma;j&semi;   rows&comma;colsLinearAlgebra:-Dimensionsimg&semi;   niters&semi;   finalMatriximg&semi;   if modprows&comma;2iters<>0ormodpcols&comma;2iters<>0then      error img's dimensions are not divisible by %1&comma;2iters   end if&semi;   tempMatrixrtable_dimsfinal&comma;datatype&equals;float8&semi;   for n from iters by 1 to 1 do      InverseDiscreteWaveletTransformcols2n1&comma;rows2n1&comma;final&comma;2&comma;w2&comma;w1&comma;temp&comma;storagetype&equals;singlearray&semi;      InverseDiscreteWaveletTransformrows2n1&comma;cols2n1&comma;temp&comma;1&comma;w2&comma;w1&comma;final&comma;storagetype&equals;singlearray&semi;   end do&semi;   return final end proc&colon;

count0 procimg&comma; $          local rows&comma;cols&comma;i&comma;j&comma;count&semi;   rows&comma;colsLinearAlgebra:-Dimensionsimg&semi;      count0&semi;   for j to cols do      for i to rows do         if imgi&comma;j&equals;0then           countcount&plus;1         end if&semi;      end do&semi;   end do&semi;   return round100 countrowscols&semi; end proc&colon;

img MatrixToGrayscaleReadcatkerneloptsdatadir&comma;/images/FingerPrint.jpg

img:= 240 x 256 MatrixData Type: float8Storage: rectangularOrder: C_order

(6.1.1)

VectorDWT procV&comma;iters&comma;w1&comma;w2&comma; $          local i&comma;s&comma;temp1&comma;temp2&semi;   sLinearAlgebra:-DimensionsV&semi;   if modps&comma;2iters<>0then      error the length of V is not divisible by enough powers of 2   end if&semi;   temp1VectorV&semi;   temp2VectorV&semi;   for i to iters do      DiscreteWaveletTransforms2i1&comma;temp1&comma;w2&comma;w1&comma;storagetype&equals;singlearray&comma;temp2&semi;      ArrayTools:-Copytemp2&comma;temp1&semi;   end do&semi;   return temp1 end proc&colon;

InverseVectorDWT procV&comma;iters&comma;w1&comma;w2&comma; $         local s&comma;i&comma;temp1&comma;temp2&semi;   sLinearAlgebra:-DimensionsV&semi;   if modps&comma;2&Hat;iters<>0then      error the length of V is not divisible by enough powers of 2   end if&semi;   temp1VectorV&semi;   temp2VectorV&semi;   for i from iters by 1 to 1 do      InverseDiscreteWaveletTransforms&sol;2&Hat;i1&comma;temp1&comma;w2&comma;w1&comma;storagetype&equals;singlearray&comma;temp2&semi;      ArrayTools:-Copytemp2&comma;temp1&semi;   end do