Example Worksheet - Wavelet Transforms

All Products Maple MapleSim

Home : Support : Online Help : Mathematics : Discrete Mathematics : DiscreteTransforms : Example Worksheet - Wavelet Transforms

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

$\{φ (x), ψ (2 x), ψ (2 x - 1), ψ (4 x), .. &period;, ψ (4 x - 3), ψ (8 x), .. &period;\}$

form an orthonormal basis of $L^{2} [0, 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) = \sum_{n = - \infty}^{\infty} h_{n} φ (2 x - n)$ (1)

$ψ (x) = \sum_{n = - \infty}^{\infty} g_{n} φ (2 x - n)$ (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 $h_{n}$ and $g_{n}$ , when $φ$ has norm 1.

$\sum_{n = - \infty}^{\infty} h_{n} h_{n + 2 h} = δ_{0, k}$ (3)

$\sum_{n = - \infty}^{\infty} g_{n} g_{n + 2 h} = δ_{0, k}$ (4)

$\sum_{n = - \infty}^{\infty} h_{n} g_{n + 2 h} = 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,

$Φ (w) = (\sum_{n = - \infty}^{\infty} \frac{h_{n} e^{- i n w}}{2}) Φ (\frac{w}{2})$

is obtained. This is the motivation for defining:

$m_{0} (x) = \sum_{n = - \infty}^{\infty} \frac{h_{n} x^{n}}{2}$

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

$m_{0} (x) m_{0} (- x) + m_{0} (x + π) m_{0} (- x - π) = 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} [0, 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 $f (x) = e^{(a x + b) I} |f (x)|$ . 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 $P (Z (X))$ , where $Z (X) = \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 nonzero imaginary part.

$A ≔ 7 &colon;$ $P ≔ add (binomial (A + k - 1, A - 1) 2^{- k} X^{k}, k = 0 .. A - 1) &semi;$ $sols ≔ [RootFinding :- Analytic (P, X, re = - 100 .. 100, im = - 100 .. 100)] &colon;$ $T := map (t \to op ([1 - t + {(- 2 t + t^{2})}^{\frac{1}{2}}, 1 - t - {(- 2 t + t^{2})}^{\frac{1}{2}}]), sols) &colon;$ $Tplot ≔ [seq ([&real; (T_{i}), &Im; (T_{i})], i = 1 .. nops (T))] &colon;$ $plots [listplot] (Tplot, style = point) &semi;$ $S ≔ select (z \to |z| < 1 and 0 \leq &Im; (z), T) &colon;$

$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 (Z (X))$ . This means that you want to find a $p (X)$ such that $P (Z (X)) = p (X) p (- X)$ . This is done by using the Feyer Reiz algorithm. For each root $λ$ of $P (X)$ , you simply have to assign one of the roots of $P (Z (X)) - λ$ to $p (x)$ . 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 ≔ proc (w, z) if type (z, realcons) then \arctan (\frac{(1 + z) \cdot \tan (\frac{w}{2})}{1 - z}) - \frac{w}{2} elif type (z, complexcons) then \arctan (\frac{(1 - {|z|}^{2}) \cdot \sin (w)}{(1 + {|z|}^{2}) \cdot \cos (w) - 2 |z| \cdot \cos (argument (z))}) - w + piecewise (\cos (w) < \frac{2 |z| \cdot \cos (argument (z))}{1 + {|z|}^{2}}, π, 0) end if end proc &colon;$

>	$minphase ≔ \infty &colon;$ $minphaseat ≔ \infty &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.

$for i from 0 to 2^{nops (S) - 1} - 1 do n ≔ i &semi; φ ≔ phase (w, S_{1}) &semi; for j from 2 to nops (S) do φ ≔ φ + 2 (modp (n, 2) - \frac{1}{2}) phase (w, S_{j}) &semi; n ≔ \frac{n - modp (n, 2)}{2} &semi; end do &semi; {PhaseList}_{i} ≔ φ &semi; myInt ≔ evalf (Int (|φ|, w = 0 .. π, method =_Gquad)) &semi; if myInt < minphase then minphase ≔ myInt &semi; minphaseat ≔ i &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.

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

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

$n ≔ minphaseat &colon;$ $if type (S_{1}, realcons) then p ≔ x - S_{1} else p ≔ (x - S_{1}) (x - \overset{&conjugate0;}{S_{1}}) end if &colon;$ $for j from 2 to nops (S) do if type (S_{j}, realcons) then p ≔ p (x - S_{j}^{(2 (modp (n, 2) - \frac{1}{2}))}) else p ≔ p (x - S_{j}^{(2 (modp (n, 2) - \frac{1}{2}))}) (x - {\overset{&conjugate0;}{S_{j}}}^{2 (modp (n, 2) - \frac{1}{2})}) end if &semi; n ≔ \frac{n - modp (n, 2)}{2} &semi; end do &colon;$

>	$pn ≔ collect (\frac{p {(1 + x)}^{A} 2^{1 - A}}{\binom{p}{} \| \binom{}{x = 1}}, x) &colon; C ≔ ListTools [Reverse] ([seq (&real; (coeff (pn, x, k)), k = 0 .. 2 A - 1)]) &colon; \frac{C}{\sqrt{add (C_{i}^{2}, i = 1 .. nops (C))}}$

$[0.002681814559, −0.001047384905, −0.01263630334, 0.03051551319, 0.06789269328, −0.04955283505, 0.01744125436, 0.5361019163, 0.7677643179, 0.2886296316, −0.1400472406, −0.1078082374, 0.004010244846, 0.01026817669]$

(3.1.1)

Verify this against the Maple function that computes wavelets.

>	$l, h ≔ DiscreteTransforms :- WaveletCoefficients (Symlet, 2 A) &colon;$ $convert (h, list)$

$[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]$

(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, $w2$ , is the $h_{n}$ , and the high pass filter, $w1$ , is the $g_{n}$ (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 down-sampling (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} = [w_{11}, w_{12}, w_{13}, w_{14}]$

$w_{2} = [w_{21}, w_{22}, w_{23}, w_{24}]$

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

$[\begin{array}{c} 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}]$

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_{i j}$ (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:	$[1, 2, 3, 4, 5, 6, 1, 2, 3, .. &period;]$
Zeros:	$[1, 2, 3, 4, 5, 6, 0, 0, 0, .. &period;]$
Reflection:	$[1, 2, 3, 4, 5, 6, 5, 4, 3, .. &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.

>	$with (DiscreteTransforms) &semi;$ $with (ImageTools) &colon;$

$[DiscreteWaveletTransform, FourierTransform, InverseDiscreteWaveletTransform, InverseFourierTransform, WaveletCoefficients, WaveletPlot]$

(5.1.1)

>	$ExW1, ExW2 ≔ WaveletCoefficients (Daubechies, 4)$

$ExW1, ExW2 ≔ [\begin{array}{c} −0.129409522551260 \\ −0.224143868042013 \\ 0.836516303737808 \\ −0.482962913144534 \end{array}], [\begin{array}{c} 0.482962913144534 \\ 0.836516303737808 \\ 0.224143868042013 \\ −0.129409522551260 \end{array}]$

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

>	$WaveletPlot (ExW1, 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.

>	$Digits ≔ 20 &colon;$ $ExW1, ExW2 ≔ WaveletCoefficients (Daubechies, 4) &semi;$ $map (identify, ExW1) &semi;$ $Digits ≔ 10 &colon;$

$ExW1, ExW2 ≔ [\begin{array}{c} −0.1294095225512603811744494188120241641744 \\ −0.2241438680420133810259727622404003554680 \\ 0.8365163037378079055752937809168732034587 \\ −0.4829629131445341433748715998644486838167 \end{array}], [\begin{array}{c} 0.4829629131445341433748715998644486838167 \\ 0.8365163037378079055752937809168732034587 \\ 0.2241438680420133810259727622404003554680 \\ −0.1294095225512603811744494188120241641744 \end{array}]$

$[\begin{array}{c} - \frac{\sqrt{3} \sqrt{2}}{8} + \frac{\sqrt{2}}{8} \\ - \frac{3 \sqrt{2}}{8} + \frac{\sqrt{3} \sqrt{2}}{8} \\ \frac{3 \sqrt{2}}{8} + \frac{\sqrt{3} \sqrt{2}}{8} \\ - \frac{\sqrt{3} \sqrt{2}}{8} - \frac{\sqrt{2}}{8} \end{array}]$

(5.1.3)

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

>	$ExV ≔ Vector (10, i \to i, datatype = float [8]) &colon;$ $ExT ≔ DiscreteWaveletTransform (ExV, Daubechies, 4) &semi;$ $InverseDiscreteWaveletTransform (ExT, Daubechies, 4)$

$ExT ≔ [\begin{array}{c} 2.22044604925031 \times 10^{−16} \\ 4.44089209850063 \times 10^{−16} \\ 4.44089209850063 \times 10^{−16} \\ 0. \\ −3.53553390593274 \end{array}], [\begin{array}{c} 2.31078903454115 \\ 5.13921615928734 \\ 7.96764328403353 \\ 10.7960704087797 \\ 12.6771540786184 \end{array}]$

$[\begin{array}{c} 1. \\ 2.00000000000000 \\ 3.00000000000000 \\ 4.00000000000000 \\ 5.00000000000000 \\ 6.00000000000000 \\ 7.00000000000000 \\ 8. \\ 9. \\ 10.0000000000000 \end{array}]$

(5.1.4)

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

$PreImg ≔ Read (cat (kernelopts (datadir), /images/lichtenstein.jpg)) &colon;$ $ExImg ≔ Matrix (ToGrayscale (PreImg)) &colon;$ $ExImgHi, ExImgLo ≔ DiscreteWaveletTransform (ExImg, 1, Daubechies, 20) &colon;$ $ExImgHiHi, ExImgHiLo, ExImgLoHi, ExImgLoLo ≔ DiscreteWaveletTransform (ExImgHi, 2, Daubechies, 20), DiscreteWaveletTransform (ExImgLo, 2, Daubechies, 20) &colon;$ $Show ≔ ⟨⟨FitIntensity (ExImgLoLo) | FitIntensity (ExImgLoHi)⟩, ⟨FitIntensity (ExImgHiLo) | FitIntensity (ExImgHiHi)⟩⟩ &colon;$

>	$Embed (Create (Show))$

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; with (DiscreteTransforms) &colon;$ $with (ImageTools) &colon;$

$ImageDWT ≔ proc (img :: (Or (Matrix, Array)), w1 :: Vector, w2 :: Vector, iters :: posint) local rows, cols, n, final, current, temp &semi; rows, cols ≔ LinearAlgebra :- Dimensions (img) &semi; n ≔ 1 &semi; final ≔ Matrix (img) &semi; if modp (rows, 2^{iters}) <> 0 or modp (cols, 2^{iters}) <> 0 then error img's dimensions are not divisible by %1, 2^{iters} end if &semi; temp ≔ Matrix (rtable_dims (final), datatype = float [8]) &semi; for n to iters do DiscreteWaveletTransform (\frac{rows}{2^{n - 1}}, \frac{cols}{2^{n - 1}}, final, 1, w2, w1, temp, storagetype = singlearray) &semi; DiscreteWaveletTransform (\frac{cols}{2^{n - 1}}, \frac{rows}{2^{n - 1}}, temp, 2, w2, w1, final, storagetype = singlearray) &semi; end do &semi; return final end proc &colon;$

$InverseImageDWT ≔ proc (img :: (Or (Matrix, Array)), w1 :: Vector, w2 :: Vector, iters :: posint) local rows, cols, n, final, current, temp, temp2, i, j &semi; rows, cols ≔ LinearAlgebra :- Dimensions (img) &semi; n ≔ iters &semi; final ≔ Matrix (img) &semi; if modp (rows, 2^{iters}) <> 0 or modp (cols, 2^{iters}) <> 0 then error img's dimensions are not divisible by %1, 2^{iters} end if &semi; temp ≔ Matrix (rtable_dims (final), datatype = float [8]) &semi; for n from iters by - 1 to 1 do InverseDiscreteWaveletTransform (\frac{cols}{2^{n - 1}}, \frac{rows}{2^{n - 1}}, final, 2, w2, w1, temp, storagetype = singlearray) &semi; InverseDiscreteWaveletTransform (\frac{rows}{2^{n - 1}}, \frac{cols}{2^{n - 1}}, temp, 1, w2, w1, final, storagetype = singlearray) &semi; end do &semi; return final end proc &colon;$

$count0 ≔ proc (img, $) local rows, cols, i, j, count &semi; rows, cols ≔ LinearAlgebra :- Dimensions (img) &semi; count ≔ 0 &semi; for j to cols do for i to rows do if img [i, j] = 0 then count ≔ count + 1 end if &semi; end do &semi; end do &semi; return round (\frac{100 \cdot count}{rows \cdot cols}) &semi; end proc &colon;$

>	$img ≔ Matrix (ToGrayscale (Read (cat (kernelopts (datadir), /images/fingerprint.jpg))))$

$VectorDWT ≔ proc (V, iters, w1, w2, $) local i, s, temp1, temp2 &semi; s ≔ LinearAlgebra :- Dimensions (V) &semi; if modp (s, 2^{iters}) <> 0 then error the length of V is not divisible by enough powers of 2 end if &semi; temp1 ≔ Vector (V) &semi; temp2 ≔ Vector (V) &semi; for i to iters do DiscreteWaveletTransform (\frac{s}{2^{i - 1}}, temp1, w2, w1, storagetype = singlearray, temp2) &semi; ArrayTools :- Copy (temp2, temp1) &semi; end do &semi; return temp1 end proc &colon;$

$InverseVectorDWT ≔ proc (V, iters, w1, w2, $) local s, i, temp1, temp2 &semi; s ≔ LinearAlgebra :- Dimensions (V) &semi; if modp (s, 2 &Hat; iters) <> 0 then error the length of V is not divisible by enough powers of 2 end if &semi; temp1 ≔ Vector (V) &semi; temp2 ≔ Vector (V) &semi; for i from iters by - 1 to 1 do InverseDiscreteWaveletTransform (s / 2 &Hat; (i - 1), temp1, w2, w1, storagetype = singlearray, temp2) &semi; ArrayTools :- Copy (temp2, temp1) &semi; end do &semi; return temp1 end proc &colon;$

Image Compression

The discrete wavelet transform is the analytical core of both the FBI fingerprint compression standard and the JPEG 2000 image compression standard. The success of wavelets in the FBI's case is nothing short of stunning: fingerprints are stored in a small fraction of the space while maintaining all of their distinguishing features.

The power of wavelets in this area comes from the zero and near zero values that appear in images transformed by using wavelets. In lossy compression, these values can be set to zero, making the data much easier to compress. The examples below demonstrate the power of this method without actually performing any compression.

You have already loaded an image above.

>	$Embed (Create (img))$

Now set the level of transform to be used, the percentage of zeros desired, and the wavelet to be used.

>	$imglevels := 2 &colon;$ $imgper ≔ 85 &colon;$ $imgWLname ≔ Daubechies &colon;$ $imgWLparam ≔ 20 &colon;$

Now you can transform and threshold the data. In the thresholding step, the smallest (in absolute value) per percent of the values in the transformed data are set to zero. As a check, output the percentage of zeros in R, the thresholded, transformed data.

$imgT ≔ ImageDWT (img, WaveletCoefficients (imgWLname, imgWLparam), imglevels) &colon;$ $imgthresh ≔ Statistics [Percentile] (ListTools :- Flatten (convert (map (abs, imgT), listlist)), imgper) &colon;$ $imgR := Matrix (map (z \to `if` (|z| \leq imgthresh, 0, z), imgT), datatype = float [8]) &colon;$

>	$count0 (imgR)$

$85$

(6.2.1)

>	$Embed (Create (imgT))$

And now you can invert the transformation, using the thresholded data.

>	$imgNew ≔ InverseImageDWT (imgR, WaveletCoefficients (imgWLname, imgWLparam), imglevels) &colon;$ $Side ≔ ⟨img \| imgNew⟩ &colon;$

>	$Embed (Create (Side))$

The quality of the new image is amazing, given that it was constructed from data that was per percent zero!

In addition to visually judging the quality of the reconstructed image, you can compute the average error (relative to the original image), view a black and white Matrix representing high error pixels, and use the ImageTools:-Quality command.

>	$\frac{LinearAlgebra :- Norm (imgNew - img, \infty)}{rtable_num_elems (img)} &semi;$ $Er := Matrix (LinearAlgebra :- Dimensions (img), (i, j) \to `if` (0.05 < \|imgNew [i, j] - img [i, j]\|, 1, 0)) &colon;$ $Quality (img, imgNew)$

$0.0002205079906$

$0.00288222335046411059$

(6.2.2)

>	$Embed (Create (Er))$

Given the ability to transform images to produce lots of zeroes without significantly affecting image quality, the applications of wavelets in image compression are obvious.

Signal Denoising

First create a signal, $S$ , and a noisy version, $NS$ .

>	$Signal := Vector (128, i \to evalf (\sin (\frac{4 i π}{128}) + \frac{1 \sin (\frac{12 i π}{128})}{3}), datatype = float [8]) &colon;$ $NS := Vector (128, i \to {Signal}_{i} + \frac{1 rand ()}{5000000000000}, datatype = float [8]) &colon;$ $plots [listplot] (NS)$

Now transform the data, threshold it, and perform the inverse transform.

>	$SDlevels ≔ 3 &colon;$ $SDthresh ≔ 0.25 &colon;$ $SDw1, SDw2 ≔ WaveletCoefficients (Coiflet, 4) &colon;$ $SDT ≔ VectorDWT (NS, SDlevels, SDw1, SDw2) &colon;$ $plots [listplot] (SDT)$

>	$SDR := Vector (map (z \to `if` (\|z\| \leq SDthresh, 0, z), SDT), datatype = float [8]) &colon;$ $plots [listplot] (SDR)$

>	$SDNew := InverseVectorDWT (SDR, SDlevels, SDw1, SDw2) &colon;$ $plots [listplot] (SDNew)$

Amazing!

Matrix-Vector Multiplication

In the situation where a fixed Matrix $M$ must be multiplied by many Vectors $V$ , the discrete wavelet transform can sometimes be used to speed up the calculations. The algorithm presented below preprocesses $M$ by using wavelets to exploit patterns. If $M$ consists of random values, there will be no speed up. By thresholding, this algorithm tries to reduce the number of floating-point multiplications required, but does not represent any decrease in algorithmic complexity.

Note that orthogonal wavelets must be used in this algorithm. This application exploits the fact that, for orthogonal wavelets, the discrete wavelet transform is an orthogonal linear transformation. (See the section Discrete Wavelet Transform for a brief comment on the orthogonality of the discrete wavelet transform.)

>	$TV ≔ Vector [column] (LinearAlgebra :- RandomVector (200), datatype = float [8]) &colon;$ $TM ≔ Matrix (200, 200, datatype = float [8], (i, j) \to evalf (\sin (i + j)), datatype = float [8]) &colon;$

Because the discrete wavelet transform is an orthogonal linear transformation, you can simultaneously transform the rows of the Matrix M and the whole Vector $V$ , while leaving their inner product unchanged. This is related to the fact that, for Vectors $x$ and $y$ and an orthogonal Matrix $A$ ,

$A x &period; A y = (A &period; A^{T} &period; x) &period; y$ = x.y

The fact that the transformation does not affect $M &period; V$ can be verified numerically.

>	$Tw1, Tw2 ≔ WaveletCoefficients (Daubechies, 4) &colon;$ $TV1 ≔ VectorDWT (TV, 1, Tw1, Tw2) &colon;$ $TM1 ≔ DiscreteWaveletTransform (TM, 2, Tw2, Tw1, storagetype = singlearray) &colon;$ $LinearAlgebra :- Norm (TM1 &period; TV1 - TM &period; TV)$

$1.59161572810262442 \times 10^{−12}$

(6.4.1)

It is also possible to transform the columns of $M$ , but then the result has to be inverse transformed.

>	$TM2 ≔ DiscreteWaveletTransform (TM1, 1, Tw1, Tw2, storagetype = singlearray) &colon;$ $TRes ≔ InverseDiscreteWaveletTransform (TM2 &period; TV1, Tw1, Tw2, storagetype = singlearray) &colon;$ $LinearAlgebra :- Norm (TRes - TM &period; TV)$

$1.70530256582424045 \times 10^{−12}$

(6.4.2)

So this means that you can effectively apply the same techniques that were used to transform and threshold images to simplify Matrix-Vector multiplication.

References

Introductory books

•	Aboufadel, Edward and Schlicker, Steven. Discovering Wavelets. Interscience, 1999.

Intermediate books; note that these books use fourier transforms extensively

•	Mallat, Stephane. A Wavelet Tour of Signal Processing. Academic Press, 1999.

•	Hernandez, Eugenio and Weiss, Guido. First course on wavelets. CRC, 1996.

The definitive book on wavelets from an expert standpoint

•	Daubechies, Ingrid. "Ten Lectures on Wavelets." SIAM, 1992.

Return to Index for Example Worksheets

Download Help Document

Maple

Maple Add-Ons

MapleSim

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim

Maple

Powerful math software that is easy to use

Maple Add-Ons

MapleSim

Advanced System Level Modeling

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim