Wavelets and Applications - Maple Programming Help

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,  is used as an orthonormal basis of ${L}^{2}\left[0,1\right]$.  In wavelet analysis, a father wavelet $\mathrm{φ}$ and a mother wavelet $\mathrm{ψ}$ are chosen such that:

form an orthonormal basis of ${L}^{2}\left[0,1\right]$. In theory, $\mathrm{\phi }$ is chosen to satisfy the conditions of a multiresolutional analysis (MRA), and then $\mathrm{\psi }$  is determined from $\mathrm{\phi }$ and the MRA. In practice, $\mathrm{\phi }$ and $\mathrm{\psi }$ are assumed to satisfy the following functional equations, and the coefficients are computed according to the desired properties of the MRA.

$\mathrm{\phi }\left(x\right)=\sum _{n=-\infty }^{\infty }\phantom{\rule[-0.0ex]{5.0px}{0.0ex}}{h}_{n}\mathrm{\phi }\left(2x-n\right)$    (1)

$\mathrm{\psi }\left(x\right)=\sum _{n=-\infty }^{\infty }\phantom{\rule[-0.0ex]{5.0px}{0.0ex}}{g}_{n}\mathrm{\phi }\left(2x-n\right)$    (2)

In fact, often $\mathrm{\psi }$ and $\mathrm{\phi }$ 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 $\mathrm{\psi }$ and $\mathrm{\phi }$ must be orthogonal reduces to the following numerical conditions on the ${h}_{n}$ and ${g}_{n}$, when $\mathrm{\phi }$ has norm 1.

$\sum _{n=-\infty }^{\infty }\phantom{\rule[-0.0ex]{5.0px}{0.0ex}}{h}_{n}{h}_{n+2h}={\mathrm{\delta }}_{0,k}$    (3)

$\sum _{n=-\infty }^{\infty }\phantom{\rule[-0.0ex]{5.0px}{0.0ex}}{g}_{n}{g}_{n+2h}={\mathrm{\delta }}_{0,k}$    (4)

$\sum _{n=-\infty }^{\infty }\phantom{\rule[-0.0ex]{5.0px}{0.0ex}}{h}_{n}{g}_{n+2h}=0$    (5)

Usually, the ${h}_{n}$  are computed, then the ${g}_{n}$ are determined by reversing the ${h}_{n}$ and negating every term.  The ${h}_{n}$ and ${g}_{n}$ 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, $\mathrm{Φ}\left(w\right)=\left(\sum _{n=-\mathrm{∞}}^{\mathrm{∞}}\phantom{\rule[-0.0ex]{5.0px}{0.0ex}}\frac{{h}_{n}{e}^{-inw}}{2}\right)\mathrm{Φ}\left(\frac{w}{2}\right)$ is obtained. This is the motivation for defining: ${m}_{0}\left(x\right)=\sum _{n=-\mathrm{∞}}^{\mathrm{∞}}\phantom{\rule[-0.0ex]{5.0px}{0.0ex}}\frac{{h}_{n}{x}^{n}}{2}$ This is known as the Z transform associated with $\mathrm{\phi }$. In this context, the orthogonality conditions reduce to:   ${m}_{0}\left(x\right){m}_{0}\left(-x\right)+{m}_{0}\left(x+\mathrm{\pi }\right){m}_{0}\left(-x-\mathrm{\pi }\right)=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 ${L}^{2}\left[0,1\right]$ 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 . 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 , 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 $P\left(Z\left(X\right)\right)$, where $Z\left(X\right)=\frac{2-X-{X}^{-1}}{2}$. 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.

 > $A≔7:$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}P≔\mathrm{add}\left(\mathrm{binomial}\left(A+k-1,A-1\right){2}^{-k}{X}^{k},k=0..A-1\right);$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{sols}≔\left[\mathrm{RootFinding}:-\mathrm{Analytic}\left(P,X,\mathrm{re}=-100..100,\mathrm{im}=-100..100\right)\right]:\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}$$T:=\mathrm{map}\left(t→\mathrm{op}\left(\left[1-t+{\left(-2t+{t}^{2}\right)}^{\frac{1}{2}},1-t-{\left(-2t+{t}^{2}\right)}^{\frac{1}{2}}\right]\right),\mathrm{sols}\right):$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{Tplot}≔\left[\mathrm{seq}\left(\left[\mathrm{ℜ}\left({T}_{i}\right),\mathrm{ℑ}\left({T}_{i}\right)\right],i=1..\mathrm{nops}\left(T\right)\right)\right]:$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{plots}\left[\mathrm{listplot}\right]\left(\mathrm{Tplot},\mathrm{style}=\mathrm{point}\right);$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}S≔\mathrm{select}\left(z→\left|z\right|<1\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}\mathbf{and}\phantom{\rule[-0.0ex]{0.5em}{0.0ex}}0\le \mathrm{ℑ}\left(z\right),T\right):$
 ${P}{:=}{1}{+}\frac{{7}}{{2}}{}{X}{+}{7}{}{{X}}^{{2}}{+}\frac{{21}}{{2}}{}{{X}}^{{3}}{+}\frac{{105}}{{8}}{}{{X}}^{{4}}{+}\frac{{231}}{{16}}{}{{X}}^{{5}}{+}\frac{{231}}{{16}}{}{{X}}^{{6}}$

Now perform spectral factorization on $P\left(Z\left(X\right)\right)$. This means that you want to find a $p\left(X\right)$ such that . This is done by using the Feyer Reiz algorithm. For each root $\mathrm{\lambda }$ of $P\left(X\right)$, you simply have to assign one of the roots of $P\left(Z\left(X\right)\right)-\mathrm{λ}$ to $p\left(x\right)$. 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.

 >
 > 

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, $\mathrm{PhaseList}$, so that you can graph them later. Remember that minimal phase is defined by minimum L1 norm.

 >

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

 > $\mathrm{plots}\left[\mathrm{display}\right]\left(\mathrm{plot}\left(\left[\mathrm{seq}\left(\mathrm{PhaseList}\left[i\right],i=0..{2}^{\mathrm{nops}\left(S\right)-1}-1\right)\right],w=0..\mathrm{π}\right),\mathrm{plot}\left(\mathrm{PhaseList}\left[\mathrm{minphaseat}\right],w=0..\mathrm{π},\mathrm{thickness}=5\right)\right)$

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

 > 
 >
 $\left[{0.002681814559}{,}{-}{0.001047384898}{,}{-}{0.01263630337}{,}{0.03051551326}{,}{0.06789269364}{,}{-}{0.04955283511}{,}{0.01744125457}{,}{0.5361019159}{,}{0.7677643179}{,}{0.2886296316}{,}{-}{0.1400472406}{,}{-}{0.1078082374}{,}{0.004010244846}{,}{0.01026817669}\right]$ (3.1.1)

Verify this against the Maple function that computes wavelets.

 > $\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{convert}\left(h,\mathrm{list}\right)$
 $\left[{0.00268181456826015}{,}{-}{0.00104738488867974}{,}{-}{0.0126363034032406}{,}{0.0305155131658779}{,}{0.0678926935012206}{,}{-}{0.0495528349370428}{,}{0.0174412550868357}{,}{0.536101917090569}{,}{0.767764317004883}{,}{0.288629631750648}{,}{-}{0.140047240442934}{,}{-}{0.107808237703290}{,}{0.00401024487152240}{,}{0.0102681767084648}\right]$ (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 ${h}_{n}$ and ${g}_{n}$ from equations (1) to (5). The low pass filter, $\mathrm{w2}$, is the ${h}_{n}$, and the high pass filter, $\mathrm{w1}$, is the ${g}_{n}$ (in vector form). In almost all useful cases, these are finite. The size of $\mathrm{w1}$ and $\mathrm{w2}$ is called the filter length. If $\mathrm{w1}$ has $n$ elements, $\mathrm{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 $\mathrm{w1}$ and $\mathrm{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:

${w}_{1}=\left[{w}_{11},{w}_{12},{w}_{13},{w}_{14}\right]$

${w}_{2}=\left[{w}_{21},{w}_{22},{w}_{23},{w}_{24}\right]$

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

$\left[\begin{array}{cccccc}{w}_{11}& {w}_{12}& {w}_{13}& {w}_{14}& 0& 0\\ {w}_{21}& {w}_{22}& {w}_{23}& {w}_{24}& 0& 0\\ 0& 0& {w}_{11}& {w}_{12}& {w}_{13}& {w}_{14}\\ 0& 0& {w}_{21}& {w}_{22}& {w}_{23}& {w}_{24}\\ {w}_{13}& {w}_{14}& 0& 0& {w}_{11}& {w}_{12}\\ {w}_{23}& {w}_{24}& 0& 0& {w}_{21}& {w}_{22}\end{array}\right]$

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 ${w}_{ij}$ (recall that ${w}_{1}$ and ${w}_{2}$ 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: $\left[1,2,3,4,5,6,1,2,3,...\right]$ Zeros: $\left[1,2,3,4,5,6,0,0,0,...\right]$ Reflection: $\left[1,2,3,4,5,6,5,4,3,...\right]$

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:

- 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.

 > $\mathrm{with}\left(\mathrm{DiscreteTransforms}\right);$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{with}\left(\mathrm{ImageTools}\right):$
 $\left[{\mathrm{DiscreteWaveletTransform}}{,}{\mathrm{FourierTransform}}{,}{\mathrm{InverseDiscreteWaveletTransform}}{,}{\mathrm{InverseFourierTransform}}{,}{\mathrm{WaveletCoefficients}}{,}{\mathrm{WaveletPlot}}\right]$ (5.1.1)
 >
 ${\mathrm{ExW1}}{,}{\mathrm{ExW2}}{:=}\left[\begin{array}{c}{-}{0.129409522551260}\\ {-}{0.224143868042013}\\ {0.836516303737808}\\ {-}{0.482962913144534}\end{array}\right]{,}\left[\begin{array}{c}{0.482962913144534}\\ {0.836516303737808}\\ {0.224143868042013}\\ {-}{0.129409522551260}\end{array}\right]$ (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.

 > $\mathrm{WaveletPlot}\left(\mathrm{ExW1},\mathrm{ExW2}\right)$

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.

 > $\mathrm{Digits}≔20:$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{ExW1},\mathrm{ExW2}≔\mathrm{WaveletCoefficients}\left(\mathrm{Daubechies},4\right);$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{map}\left(\mathrm{identify},\mathrm{ExW1}\right);$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{Digits}≔10:$
 ${\mathrm{ExW1}}{,}{\mathrm{ExW2}}{:=}\left[\begin{array}{c}{-}{0.1294095225512603811744494188120241641744}\\ {-}{0.2241438680420133810259727622404003554680}\\ {0.8365163037378079055752937809168732034587}\\ {-}{0.4829629131445341433748715998644486838167}\end{array}\right]{,}\left[\begin{array}{c}{0.4829629131445341433748715998644486838167}\\ {0.8365163037378079055752937809168732034587}\\ {0.2241438680420133810259727622404003554680}\\ {-}{0.1294095225512603811744494188120241641744}\end{array}\right]$
 $\left[\begin{array}{c}{-}\frac{{1}}{{8}}{}\sqrt{{6}}{+}\frac{{1}}{{8}}{}\sqrt{{2}}\\ {-}\frac{{3}}{{8}}{}\sqrt{{2}}{+}\frac{{1}}{{8}}{}\sqrt{{6}}\\ \frac{{3}}{{8}}{}\sqrt{{2}}{+}\frac{{1}}{{8}}{}\sqrt{{6}}\\ {-}\frac{{1}}{{8}}{}\sqrt{{6}}{-}\frac{{1}}{{8}}{}\sqrt{{2}}\end{array}\right]$ (5.1.3)

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

 > $\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{ExT}≔\mathrm{DiscreteWaveletTransform}\left(\mathrm{ExV},\mathrm{Daubechies},4\right);$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{InverseDiscreteWaveletTransform}\left(\mathrm{ExT},\mathrm{Daubechies},4\right)$
 ${\mathrm{ExT}}{:=}\left[\begin{array}{c}{2.22044604925031}{}{{10}}^{{-16}}\\ {4.44089209850063}{}{{10}}^{{-16}}\\ {4.44089209850063}{}{{10}}^{{-16}}\\ {0.}\\ {-}{3.53553390593274}\end{array}\right]{,}\left[\begin{array}{c}{2.31078903454115}\\ {5.13921615928734}\\ {7.96764328403353}\\ {10.7960704087797}\\ {12.6771540786184}\end{array}\right]$
 $\left[\begin{array}{c}{1.}\\ {2.00000000000000}\\ {3.00000000000000}\\ {4.00000000000000}\\ {5.00000000000000}\\ {6.00000000000000}\\ {7.00000000000000}\\ {8.}\\ {9.}\\ {10.0000000000000}\end{array}\right]$ (5.1.4)

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

 > $\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{ExImg}≔\mathrm{Matrix}\left(\mathrm{ToGrayscale}\left(\mathrm{PreImg}\right)\right):$
 > $\mathrm{Embed}\left(\mathrm{Create}\left(\mathrm{Show}\right)\right)$

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.

 > $\mathrm{restart};\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{with}\left(\mathrm{DiscreteTransforms}\right):$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{with}\left(\mathrm{ImageTools}\right):\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}$
 >
 >
 >
 >
 ${\mathrm{img}}{:=}\left[\begin{array}{c}{\mathrm{240 x 256}}{\mathrm{Matrix}}\\ {\mathrm{Data Type:}}{{\mathrm{float}}}_{{8}}\\ {\mathrm{Storage:}}{\mathrm{rectangular}}\\ {\mathrm{Order:}}{\mathrm{C_order}}\end{array}\right]$ (6.1.1)
 >
 >