How to get started with wavelets and Wavelab

From CCM
Jump to: navigation, search

Wavelab is a free collection of Matlab routines for computation with wavelets.

Contents

Installation

  • Download from http://www-stat.stanford.edu/software/wavelab
  • Expand into a directory. I used /Users/jmandel/Matlab/Wavelab850. Change to this directory.
  • chmod u+w WavePath.m
  • edit the file WavePath.m and add the line to set WAVELABPATH just before while (exist(WAVELABPATH)~=7) (note the / at the end):
WAVELABPATH='/Users/jmandel/Matlab/Wavelab850/'
  • mv startup.m wavelab_startup.m
  • Edit your startup.m file (mine is /Users/jmandel/Matlab/startup.m) and add the Wavelab directory to your path:
addpath /Users/jmandel/Matlab/Wavelab850
  • Start Matlab and run wavelab_startup

Examples of wavelets

  • plot(MakeWavelet(4,8,'Coiflet',5,'Mother',512)) 4-8-Coiflet-5-Mother.png
  • plot(MakeWavelet(3,4,'Daubechies',4,'Mother',1024)) 3-4-Daubechies-4-Mother.png

Wavelet transform (discrete, orthogonal, periodicized)

Wavelet transform is the expansion of a vector in a basis of wavelets. The value of the transform are the coefficients of the expansion, just like in discrete Fourier transform. Wavelet transform is very efficient and in fact even faster than the fast Fourier transform.

First, we make a short constant vector, called a filter, which is used much like a stencil in finite differences to compute the wavelet transforms. This vector encodes all there is to know about the chosen wavelet:

qmf = MakeONFilter('Coiflet',2)

We'll use this "filter" vector for the rest of this section.

Your first wavelet transform

Make a test vector: n=64;x=abs([1:n]-(n+1)*2/3); plot(x)

Signal-x.png

and its wavelet transform: wc = FWT_PO(x,1,qmf); plot(wc)

Wavelet-transform-of-signal-x.png

Then use the inverse transform and verify you got your vector back: xc = IWT_PO(wc,1,qmf); norm(xc-x). You should see something small, like 8e-09.

Get some wavelet pictures and the wavelet matrix

When you apply inverse wavelet transform to a standard unit vector, you get the corresponding wavelet:

i=25; n=512; wc=zeros(1,n); wc(i)=1; xc=IWT_PO(wc,1,qmf); plot(xc)

Coiflet-25.png

Try a loop over i. Can you guess how are the wavelets numbered?

Now generate matrix X that has wavelets as columns, in lower dimension to see better. First a movie of the columns one by one

n=64; for i=1:n, wc=zeros(n,1); wc(i)=1; X=zeros(n); X(:,i)=IWT_PO(wc,1,qmf); mesh(X), pause(1),end

Wavelet-column-X.png

then the whole matrix,

n=64; X=zeros(n); for i=1:n, wc=zeros(n,1); wc(i)=1; X(:,i)=IWT_PO(wc,1,qmf); end; mesh(X)

Wavelet-X.png

Taking the wavelet transform apart

Now we'll generate the contribution of each wavelet to the inverse transform one at a time, and watch how the approximation gets better. Just copy and paste this file into Matlab editor, and enjoy the show:

Increment-6.png

Note that you get also larger contributions from some wavelets split across the beginning and the end of the interval because the signal and the wavelets roll over periodically, and the signal has a jump there:

Increment-5.png

Dropping small terms in the wavelet expansion is the main idea of wavelet compression.

The wavelet transform matrix

Apply the wavelet transform to the columns of the identity matrix

n=64;XI=eye(n);W=zeros(n);for i=1:n,W(:,i)=FWT_PO(XI(:,i),1,qmf);end, mesh(W)

Wavelet-transform-matrix.png

so that multiplication by this matrix is the same as applying the wavelet transform:

x=rand(n,1);w=FWT_PO(x,1,qmf);norm(w-W*x) (you should get numerical zero)

The wavelet transform matrix is orthogonal

norm(W*W'-eye(n)) (you should get numerical zero)

and of course, W=transpose(X)=inv(X), where X is the wavelet matrix we got above.

Wavelet transform in 2D

A 2D array is transformed by wavelets by applying the transform in one direction, then in the other direction:

a=rand(64);
a1=zeros(n);for i=1:n,a1(:,i)=FWT_PO(a(:,i),1,qmf);end
a2=zeros(n);for i=1:n,a2(i,:)=FWT_PO(a1(i,:),1,qmf);end

The order of the dimensions does not matter. This is immediately clear because using the wavelet transform matrix W above,

a1=W*a; a2=a1*W';

which gives the same result a2=(W*a)*W'=W*(a*W')=W*a*W', regardless of the order of the row and the column transform.

Notes

  • Wavelab does not have much argument checking and error handling. When something goes wrong, such as a bad or missing argument to some routine, you get a cryptical error someplace deep inside.

Works with

This is the environment this Howto was tested with. It should work with any reasonably recent Matlab as well.

  • Wavelab .850, Matlab 2010b, Mac OSX Snow Leopard

See also

The following references happen to be all local from my town, Boulder, Colorado

Other references and links:

Personal tools
Namespaces

Variants
Actions
Navigation
Toolbox