The Haar wavelet transform

Kristian Sandberg
Dept. of Applied Mathematics
University of Colorado at Boulder

A fractal along with a two level Haar (standard) decomposition. Image source: Carlson's Fractal Gallery.

1 Goal

The goals with this lab are:

1.

Understand the properties of the Haar basis.

2.

Understand the averaging and differencing technique and how it relates to a wavelet decomposition.

3.

Use the Haar transform to compress and de-noise one dimensional signals and images.

2 Introduction

This lab will illustrate some fundamental concepts about wavelets and multiresolution theory. Most examples and definitions will be made with respect to the Haar scaling function and the Haar wavelet. However, the concepts covered in this lab are true for more ``advanced'' wavelets such as the Daubechies wavelets.

3 The Haar basis

Definition 1 (The Haar scaling function)   Let $\phi:\ \mathbb{R}\rightarrow \mathbb{R}$ be defined by

$$\phi(t)=\left\{ \begin{array}{ll} 1 & t\in [0,1) \\ 0 & t\not\in [0,1) \end{array} \right.$$ (1)

Define $\phi^j_i:\ \mathbb{R}\rightarrow \mathbb{R}$ as

$$\phi^j_i(t)=\sqrt{2^j}\phi(2^jt-i),\ j=0,1,\ldots\mbox{ and }i=0,1,\ldots,2^j-1\ .$$ (2)

Define the vector space V^j as

$$V^j=\mbox{\em {sp}}\{\phi^j_i\}_{i=0,\ldots,2^j-1}$$ (3)

where $\mbox{\em {sp}}$ denotes the linear span.

The index j refers to dilation and i refers to translation.

Warning! This definition may be different in different texts! Sometimes one defines the scaling function for negative integers as well and uses a different sign convention. Whenever you use a text on wavelets, make sure you know what sign convention is being used.

The constant $\sqrt{2^j}$ is chosen such that $<\phi^j_i,\phi^j_i>=\int_0^1 \phi^j_i(t)^2\ dt=1$. If one considers the scaling function on other intervals than [0,1], the normalization constant will change.

Exercise 1: The Haar scaling function

a) Using, e.g., the Matlab functions haarscale.m and haargrid.m, plot the basis vectors for V^j for j=0,1 and 2. Describe the significance of the index j in V^j and the significance of the index i in $\phi^j_i$ when j is fixed.

b) Claim: The following property holds for V^j:

$$V^j\subseteq V^{j+1}\ .$$ (4)

Based on your plots in a) above, convince yourself (and the reader of your lab report) that this property holds for the Haar scaling function. (If you see this relation written in ``the opposite direction'' in some text, they used a different sign convention. The relation given above is consistent with Definition 1.)

Definition 2 (The Haar wavelet function)   Let $\psi:\ \mathbb{R}\rightarrow \mathbb{R}$ be defined by

$$\psi(t)=\left\{ \begin{array}{ll} 1 & t\in [0,1/2) \\ -1 & t\in [1/2,1) \\ 0 & t\not\in [0,1) \end{array} \right.$$ (5)

Define $\psi^j_i:\ \mathbb{R}\rightarrow \mathbb{R}$ as

$$\psi^j_i(t)=\sqrt{2^j}\psi(2^jt-i),\ j=0,1,\ldots\mbox{ and }i=0,1,\ldots,2^j-1\ .$$ (6)

Define the vector space W^j

$$W^j=\mbox{\em {sp}}\{\psi^j_i\}_{i=0,\ldots,2^j-1}$$ (7)

where $\mbox{\em {sp}}$ denotes the linear span.

Warning! This definition may be different in different texts! Sometimes one defines the wavelet function for negative integers as well and uses a different sign convention. Whenever you use a text on wavelets, make sure you know what sign convention is being used.

Again, the constant $\sqrt{2^j}$ is chosen such that $<\psi^j_i,\psi^j_i>=\int_0^1 \psi^j_i(t)^2\ dt=1$. If one considers the wavelet function on other intervals than [0,1], the normalization constant will change.

Exercise 2: The Haar wavelet

a) Using, e.g., the Matlab functions haargrid.m and haarwavelet.m, plot the basis vectors for W^j for j=0,1 and 2. Describe the significance of the index j in W^j and the significance of the index i in $\psi^j_i$ when j is fixed.

b) Claim: The following property holds for W^j:

$$W^j\subseteq V^{j+1}\ .$$ (8)

Based on your plots in a) above, convince yourself (and the reader of your lab report) that this property holds for the Haar wavelet function. (If you see this relation written in ``the opposite direction'' in some text, they used a different sign convention. The relation given above is consistent with Definition 2.)

c) A standard argument for using wavelets rather than the Fourier basis for signal processing is that the Haar wavelet is localized. By looking at your plots from a) above, can you see why the Haar wavelet is called a ``localized'' function? Is the Fourier basis localized?

4 The structure of the Haar basis

The Haar basis (and other wavelet basis) has the very important property that

$$V^{j+1}=V^j\oplus W^j\ .$$ (9)

This property is explored in the following exercise.

Exercise 3: The relation between V^j and W^j

Using, e.g., the Matlab functions haargrid.m and haarwavelet.m, produce a few plots that illustrate the following relation:

$$V^{j+1}=V^j\oplus W^j\ .$$ (10)

(The formulation of this exercise is very vague. However, the purpose is to visualize the important relation $V^{j+1}=V^j\oplus W^j$. Use your imagination to do this.)

5 Expansion of a signal in an ON-basis

Let $f\in\mathbb{R} ^N$ denote a finite signal with N samples (cf. an N-dimensional vector) given by

$$f=(f[0],f[1],\ldots,f[N-1])\ .$$ (11)

Let $\{v_k\}_{k=0,1,\ldots,N-1}$ be an ON-basis in $\mathbb{R} ^N$.Then

$$f=<f,v_0>v_0+<f,v_1>v_1+\ldots+<f,v_{N-1}>v_{N-1}\ .$$ (12)

Let us first review how this was used for Fourier analyzing a signal and then see how this can be used for wavelet expansion of a signal.

5.1 Review of Fourier analyzing a signal

Let

$$v_k=\frac{1}{\sqrt{N}}(1,e^{2\pi i\frac{k}{N}},e^{2\pi i\frac... ...s,e^{2\pi i\frac{nk}{N}},\ldots,e^{2\pi i\frac{(N-1)k}{N}})\ .$$ (13)

Denote the (n+1):th element of the v_k vector as v_k[n]. Then f can be expanded as in (12) with

$$<f,v_k>=\sum_{n=0}^{N-1}f[n]v_k[n]^*=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}f[n]e^{-2\pi i\frac{nk}{N}}$$ (14)

where ^* denotes complex conjugation. Numerically, we compute <f,v_k> by using the Fast Fourier Transform.

  
5.2 Wavelet transformation of a signal

In this example we consider a signal with 8 samples for clarity. However, the procedure can be generalized to any finite signal. The analysis simplifies if the length (dimension) of the signal is 2^k for some positive integer k.

Let $f\in\mathbb{R} ^8$ and expand this signal into $V^0\oplus W^0\oplus W^1\oplus W^2$ spanned by the basis vectors according to Definitions 1 and 2 (discretized). Then

$$\begin{split} f &= <f,\phi^0_0>\phi^0_0+<f,\psi^0_0>\psi^0_0... ...si^2_1+<f,\psi^2_2>\psi^2_2+<f,\psi^2_3>\psi^2_3\ . \end{split}$$ (15)

The correspondence to Fourier expanding a signal is now to compute the inner products $<f,\phi^j_i>$ and $<f,\psi^j_i>$. In principle, we could do this by evaluating an integral (or a sum in the discrete case). In practice, computation of the inner products for a wavelet basis is done as illustrated by the following example.

  
5.2.1 Example: Averaging and differencing of a finite signal

Consider f=(2,5,8,9,7,4,-1,1). We wish to expand this signal in the Haar basis according to (15). In practice, we do this by performing the following steps:

Step 1:
$f=(2+5,8+9,7+4,-1+1,2-5,8-9,7-4,-1-1)/\sqrt{2}=$
$=(7,17,11,0,-3,-1,3,-2)/\sqrt{2}$

Step 2:
$f=(\frac{7+17}{\sqrt{2}},\frac{11+0}{\sqrt{2}},\frac{7-17}{\sqrt{2}},\frac{11-0}{\sqrt{2}},-3,-1,3,-2)/\sqrt{2}=$
$=(\frac{24}{\sqrt{2}},\frac{11}{\sqrt{2}},\frac{-10}{\sqrt{2}},\frac{11}{\sqrt{2}},-3,-1,3,-2)/\sqrt{2}$

Step 3:
$f=(\frac{24+11}{(\sqrt{2})^2},\frac{24-11}{(\sqrt{2})^2},\frac{-10}{\sqrt{2}},\frac{11}{\sqrt{2}},-3,-1,3,-2)/\sqrt{2}=$
$=(\frac{35}{2},\frac{13}{2},\frac{-10}{\sqrt{2}},\frac{11}{\sqrt{2}},-3,-1,3,-2)/\sqrt{2}\simeq$
$\simeq(12.4,4.60,-5.00,5.50,-2.12,-0.707,2.12,-1.41)$

The numbers in the final vector are the coefficients in the expansion (15). You probably notice the pattern, we add/subtract a pair of numbers, divide by the normalization factor $\sqrt{2}$ and between every step we don't touch the last elements of the vector. This procedure is called ``averaging and differencing''. In the next section, we will study this example in more depth and see what we really did mathematically.

5.2.2 The mathematics of averaging and differencing

The example in the previous section shows the steps that we want to implement when we code up a one dimensional Haar wavelet transform. When we consider two dimensional signals (images) we will see that all we need to do is to act with the one dimensional transform on each row of the image followed by a transformation of each column of the image. Hence, the above example should serve as a guide for programming your own Haar wavelet transform later on.

Although the example above may be a good illustration of how the Haar transform is coded, it is not that clear how the ``averaging and differencing'' technique actually relates to the more mathematical description in (15). The next paragraph will go through the same steps as in the example above, but try to explain what really was done mathematically.

We still consider the signal f=(2,5,8,9,7,4,-1,1). The coefficients of the original signal are the coefficients of the vector f expanded in V³, i.e., $<f,\phi^3_0>=2,\ <f,\phi^3_1>=5,\ <f,\phi^3_2>=8$, etc. Recall that the basis vectors of V³ are translation of the Haar mother scaling function which has been dilated so that each basis function has a support equal to 1/8 (=1 pixel) of the support of f which is 8 pixels.

Our goal is to expand f into $V^0\oplus W^0\oplus W^1\oplus W^2$. We will do this in three steps.

Step 1: Expand f into $V^2\oplus W^2$.
The first step in the example above can be described as the matrix-vector multiplication f₁=W₁f where

$$W_1=\frac{1}{\sqrt{2}}\left( \begin{array}{cccccccc} 1 & 1 &... ... 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \\ \end{array}\right)$$ (16)

Notice that the first four rows correspond to the basis vectors $\phi^2_0,\ \phi^2_1,\ \phi^2_2$ and $\phi^2_3$ respectively which span V².

The last four rows correspond to the basis vectors $\psi^2_0,\ \psi^2_1,\ \psi^2_2$ and $\psi^2_3$ which span W².

Hence, the final vector in Step 1 in the example above is nothing more than the coefficients of f in the expansion

$$f=\frac{1}{\sqrt{2}}( 7\phi^2_0+17\phi^2_1+11\phi^2_2+0\phi^2_3 -3\psi^2_0-1\psi^2_1+3\psi^2_2-2\psi^2_3 )\ ,$$ (17)

i.e., we have expanded our signal into $V^2\oplus W^2$.

Step 2: Expand f into $V^1\oplus W^1 \oplus W^2$.
In this step we notice that we already have the coefficients for the basis vectors of W². Hence we don't have to worry about these coefficients in this step, we simply just keep these coefficients (the last four entries of our vector in the previous step) and only work with the first four entries.

The second step can be described as the matrix-vector multiplication f₂=W₂f₁ where

$$W_2=\left( \begin{array}{cccccccc} \frac{1}{\sqrt{2}} & \fra... ... \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array}\right)\ .$$ (18)

We can combine the first and the second step as f₂=W₂W₁f where

$$W_2W_1=\left( \begin{array}{cccccccc} \frac{1}{2} & \frac{1}... ...{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \\ \end{array}\right)\ .$$ (19)

Notice that the first two rows correspond to the basis vectors $\phi^1_0$ and $\phi^1_1$ which span V¹.

The third and fourth rows correspond to the basis vectors $\psi^1_0$ and $\psi^1_1$ which span W¹. The last four rows correspond to the basis vectors of W².

Hence, the final vector in Step 2 in the example above is nothing else than the coefficients of f in the expansion

$$f=\frac{1}{\sqrt{2}}( \frac{24}{\sqrt{2}}\phi^1_0+\frac{11}{\... ...qrt{2}}\psi^1_1 -3\psi^2_0+-1\psi^2_1+3\psi^2_2-2\psi^2_3 )\ .$$ (20)

Step 3: Expand f into $V^0\oplus W^0\oplus W^1\oplus W^2$.
In this step we notice that we already have the coefficients for the basis vectors of $W^1\oplus W^2$. Hence we don't have to worry about these coefficients in this step, we simply just keep these coefficients (the last six entries of our vector in the previous step) and only work with the first two entries.

The third step can be described as the matrix-vector multiplication f₃=W₃f₂ where

$$W_3=\left( \begin{array}{cccccccc} \frac{1}{\sqrt{2}} & \fra... ... \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array}\right)\ .$$ (21)

We can combine the first, second and third step as f₃=W₃W₂W₁f where

$$W_3W_2W_1= \left( \begin{array}{cccccccc} \frac{1}{2\sqrt{2... ...rac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \end{array} \right)\ .$$ (22)

Notice that the first row corresponds to the basis vector $\phi^0_0$ which spans V⁰, the second row corresponds to the basis vector $\psi^0_0$ which spans W⁰ and so on.

Hence, the final vector in Step 3 in the example above is nothing else than the coefficients of f in the expansion

$$\begin{array}{ll} f & =\frac{1}{\sqrt{2}}( \frac{35}{2}\phi^0... ...2\psi^2_0-0.707\psi^2_1+2.12\psi^2_2-1.41\psi^2_3 ) \end{array}$$ (23)

5.3 The wavelet transform as a unitary operator

The procedure in the example above can be expressed as

$$Wf=\frac{1}{\sqrt{2}}( \frac{35}{2}\phi^0_0+\frac{13}{2}\psi^... ...}{\sqrt{2}}\psi^1_1 -3\psi^2_0-1\psi^2_1+3\psi^2_2-2\psi^2_3 )$$ (24)

where W is the wavelet transform.

Exercise 4: The wavelet transform as a unitary operator

a) Verify (using, for example, Matlab or Mathematica to reproduce the example above) that the matrix in (22) describes the wavelet transform for an 8 dimensional signal.

b) Verify that the wavelet operator in (22) is a unitary operator.

c) Mention a few properties that characterize unitary operators.

5.4 Wavelets and filters

The method of averaging and differencing explained above can also be expressed as filtering the data. Averaging corresponds to a low pass filtering. It removes high frequencies of the data. Since details (sharp changes in the data) correspond to high frequencies, the averaging procedure tends to smooth the data. The low pass filter can be expressed as $\frac{1}{\sqrt{2}}(1,1)$ in the Haar case and when we average the data, we move this filter along our input data.

The differencing corresponds to high pass filtering. It removes low frequencies and responds to details of an image since details correspond to high frequencies. It also responds to noise in an image, since noise usually is located in the high frequencies. The high pass filter can be expressed as $\frac{1}{\sqrt{2}}(1,-1)$ in the Haar case and when we difference the data, we simply move this filter along our input data. The low pass and high pass filters make up what in signal processing language is referred to as a filter bank. The method of averaging and differencing is referred to as analysis. The reverse procedure (going the opposite way in the example above) is called synthesis.

Hence, the wavelet transform separates low and high frequencies, just as the Fourier transform. Since different features of a signal (background, details, noise, edges, etc.) correspond to different frequencies, this is a key to use wavelets in signal processing. The nice thing is that wavelets are localized since they only live on part of the interval of the data, as opposed to the trigonometric functions used in Fourier analysis which live on the entire interval of the data.

6 Exercises

In the previous lab on Fourier techniques in image processing, we learned that when we want to manipulate an image (compress, de-noise, etc.), we first transformed the image into the Fourier basis, performed our manipulations in this basis, and then transformed back. The same principle holds when using wavelets techniques for image processing! Now we will transform our image into a wavelet basis, perform our manipulations in this basis, and then transform back. The foundation of signal processing can be described as to find a basis which is suitable for the manipulations we need. This means that image processing relies on concepts from linear algebra such as orthogonality, norms, inner products, etc.

When doing the exercises below, there are a few things you should keep in mind. Just as you probably has developed a feeling for what a Fourier spectrum looks like and how it can be interpreted, try to get a feeling for how a signal looks like in the Haar wavelet basis. Where can you find the high/low frequencies? Where is the magnitude of the coefficients large/small?

When you wavelet decompose images, there is a trick to make it easier to study the decomposition. Usually a major part of the decomposition (see the image in the beginning of this instruction) is very dark. If you use Matlab, use the command brighten() to find a suitable brightness of the picture.

When you run your wavelet transform in Matlab, you may notice that the transform is quite slow. Even though you may have implemented your wavelet transform as a ``fast algorithm'', the architecture of the Matlab language may slow down your algorithm considerably. This has to do with the fact that loops and recursive functions calls (which you may need) are time consuming in Matlab. Don't worry about that, but keep in mind that your algorithm probably is much faster in C or Fortran which can handle loops and subroutine calls much more efficient than Matlab.

Exercise 5: The one dimensional Haar transform

While it is very important to keep in mind that the wavelet transform can be described by a unitary matrix, it is not efficient to perform the transformation by multiplying the matrix to a vector. The way to wavelet transform numerically, is to proceed as in the example in section 5.2.1 by moving the filters $(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$ and $(\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}})$ along the signal recursively.

Program a function with the following in-/output:
Input:

• A vector of length 2^k where k is a positive integer.
• The ``level'' of decomposition. (E.g., the first step in the example in section 5.2.1 corresponds to a first level decomposition, the second step corresponds to a second level decomposition and the third step corresponds to a third level decomposition.)

Output:

• A vector which is the Haar transformation of the vector you gave as input corresponding to the level you gave as an input.

The purpose of this exercise is to build a Haar transform to use for the following exercises. If you want to use slightly different input/outputs, feel free to do so as long as you can use your Haar transform efficiently.

a) Verify that your transform can reproduce the example in section 5.2.1.

b) Verify that the l²-norm is preserved in each of the three steps in the example.

c) Build an inverse Haar transform, i.e., a code that takes a Haar transformed vector and returns a non-transformed vector. (A ``reverse'' of the code in a).)

Exercise 6: Compression of a one dimensional signal

Generate a signal that samples the function

$$f(t)=\left\{ \begin{array}{ll} e^{-t}\cos(32\pi t) & t\in [0... ... t) & t\in [3/4,7/8) \\ 0 & t\in [7/8,1) \end{array} \right.$$ (25)

at at least 1024 points. Plot the wavelet decomposition of the signal at a few different levels.

By using thresholding (like in the previous lab) on the Haar transformed signal, investigate the compression performance of the Haar transform. Compare the compression performance to a compression using the FFT of the same signal.

Bonus question: If this signal represents a sound signal, what does the signal sound like?

Exercise 7: The two dimensional Haar transform

A standard decomposition of a two dimensional signal (image) is easily done by first performing a one dimensional transformation on each row followed by a one dimensional transformation of each column.

Construct a function that Haar transforms an image for different levels. Test your function on the same images as used for the previous labs (the basket and the image with campus and the Rockies). Compress the image and investigate the compression performance of the Haar transform compared to using the FFT.

(Hint. The image in the beginning of this lab instruction illustrates a level 2 Haar transformation of an image. The structure of your second level Haar transform should look like this one.)

Exercise 8: De-noising images

Using the Haar transform, try to de-noise the following image:

Fractal from Carlson's Fractal Gallery with random noise added.

(Do not expect your de-noising to do miracles. There are many approaches to de-noise. The simplest is probably to do thresholding on your wavelet decomposed image but you may come up with other methods that work better.)

Exercise 9: Horizontal, vertical and diagonal features

Down load the following image:

Do a Haar wavelet decomposition of this image at a few different levels. What are your observations? (You may have to adjust the brightness of the wavelet decomposition using brighten() in Matlab.)

Exercise 10: ``De-scratching''

Remove (or reduce) the ``scratch'' from the following image!

An MRI (Magnetic Resonance Imaging) image from First Radiology Clinic with "a scratch" added. Medical imaging is a very important application of digital image processing.

7 Format of your lab report

Write down your comments to all of the exercises above (word processed). Include a representative collection of the code(s) that you used as an appendix. Make sure to include comments in your code explaining what it does. Include any plots or images you need to justify your conclusions.

8 Questions?

In case you have questions regarding the material in this lab, do not hesitate to contact me at kristian.sandberg@colorado.edu or visit me during my lab hours Mondays 4:30-6 pm and Thursdays 4-5:30 pm in ECCR 143.

Lycka till!

About this document ...

The Haar wavelet transform

This document was generated using the LaTeX2HTML translator Version 98.1 release (February 19th, 1998)

Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.

The command line arguments were:
latex2html -split 0 -no_navigation haar.tex.

The translation was initiated by Kristian Sandberg on 2000-04-01