How to get started with wavelets and Wavelab
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
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)
and its wavelet transform: wc = FWT_PO(x,1,qmf); plot(wc)
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)
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
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)
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:
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:
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)
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 local from my town, Boulder, Colorado
- Data assimilation seminar: Aimé Fournier
- A Practical Guide to Wavelet Analysis by C. Torrence and G. P. Compo
Other references and links: