Wavelet Toolbox in MATLAB — Part 1

7 min readNov 14, 2019

Wavelet Transform is one of the main image processing methods. In this post, simple examples are presented to demonstrate how MATLAB’s Wavelet toolbox can be used for computing two-dimensional (2-D) Discrete Wavelet Transform (DWT) of an input image and displaying its coefficients. Specifically, the main topics are as follows:

2-D Discrete Wavelet Transform (DWT) algorithm
Computing DWT using dwt2 MATLAB function
Displaying coefficients
Computing DWT using wavedec2 MATLAB function
Additional notes

Let’s dive into the article.

2-D Discrete Wavelet Transform (DWT) algorithm

Figure 1: Diagram of a single-level 2-D DWT. The functions G(n) and H(n) represent the coefficients of the low-pass and high-pass filters, respectively. The figure is adopted from here.

The DWT algorithm is shown in the above figure. Single-level DWT decomposes an input image into 4 sub-bands: one approximate component (A1), and three detail components (H1, V1, D1). This procedure can be recursively applied over approximate component of each level to perform a multi-level decomposition.

In MATLAB you can use either dwt or wavedec functions to compute DWT of a one-dimensional (1-D) signal. For images (or 2-D signals), there are functions with similar names:dwt2 andwavedec2.

The dwt2 performs single-level DWT of an input image and the wavedec2 performs a multilevel DWT decomposition. There are also subtle differences between these two functions which is shown in this article.

Computing DWT using “dwt2” MATLAB function

Example 1: single-level decomposition
Here, the famous test image “barbara” is used as input. You may want to use other images. In this example, the dwt2 analysis (or decomposes) the input image into 4 sub-band images (as shown in Fig. 1 above). Then, using these sub-band images, the inverse transform is computed and the recovery error is measured.

% Get input
im = imread('barbara.tif');
im = double(im);% Compute DWT
wname = 'db1'; % Analyzing filters used to compute the 2-D DWT
[A1,H1,V1,D1] = dwt2(im,wname);% Compute inverse DWT
im_hat = idwt2(A1,H1,V1,D1,wname);% Evaluate error
MSE=mean(mean((im(:)-im_hat(:)).^2))

The decomposition coefficients are A1 , H1 , V1 , and D1 .

In dwt2 , the wname argument specifies analyzing filters. Similarly, idwt2 gets wname as input argument to specify the synthesizing filters. Analyzing and synthesizing filters are chosen in such a way that a condition called perfect reconstruction is satisfied.

To recover (or synthesize, or reconstruct) the input image from its coefficents (A1 , H1 , V1 , and D1 ), the idwt2 function can be used.

im_hat is the recovered image and the recovery error can be measured by mean square error (MSE):

MSE =1.5936e-27

Here, we see that the error is almost zero. In the next example, a deeper decomposition is applied. Therefore, one can expect the error increases.

Example 2: Two-level decomposition

Similar to the above example, first, the DWT is performed. Then, the coefficients are used to recover the image and the error is measured.

wname = 'db1';
[A1,H1,V1,D1] = dwt2(im,wname);
[A2,H2,V2,D2] = dwt2(A1,wname);% Inverse DWT
A1_hat = idwt2(A2,H2,V2,D2,wname);
im_hat = idwt2(A1_hat,H1,V1,D1,wname);% Error
MSE=mean(mean((im(:)-im_hat(:)).^2))

After performing two-level decomposition, the coefficients are A2 , H2 , V2 , D2 , H1 , V1 , and D1. The only approximate component is A2 and the others are detail components. Therefore, performing single-level inverse transform over A2 , H2 , V2 , D2 gives us an estimation of A1 (which is indicated in the above code with A1_hat ).

The recovery error is:

MSE =5.8587e-27

Notice that the error increases compared to the previous example since we have performed a deeper decomposition.

Displaying coefficients

In MATLAB, the values of gray-scale image pixels have to be in the range 0–255 (or 0–1). However, after decomposition, the approximation/wavelet coefficients are not restricted to be in this range. For example, see the following piece of code:

>> [min(A1(:)), max(A1(:))]
ans =
38.5000  487.0000>> [min(D1(:)), max(D1(:))]
ans =
-122.5000  118.5000

Therefore, in order to display coefficients as images, you have to scale coefficients properly.

In MATLAB, the best way to scale wavelet coefficients is to use wcodemat function. However, you can still use more simpler ways as shown in the following example.

Example 3: Comparing different ways of displaying coefficients as images

Let’s display the H1 wavelet coefficient in different ways:

% The 1st method (Recommended Method)
H1img = wcodemat(H1,255,'mat',1);
figure,imshow(H1img/maxI)
title('H1 - wcodemat')% The 2nd method
figure,imshow(mat2gray(H1))
title('H1 - mat2gray')% The 3rd method
H1img_log = log(1+abs(H1));
figure,imshow(mat2gray(H1img_log))
title('H1 - log and mat2gray')

Computing DWT using “wavedec2” MATLAB function

To perform a multi-level decomposition, instead of recursively applying dwt2 (as it is done in Example 2 above), one can easily use wavedec2 function.

Example 4: two-level decomposition using “wavedec2”

im = imread('barbara.tif');
im = double(im);% DWT
N = 2; % Number of decomposition levels
wname = 'db1';[C,S] = wavedec2(im,N,wname); 
% C: coefficient vector, S: bookkeeping matrix% Inverse DWT
im_hat = waverec2(C,S,wname);% Error
MSE=mean(mean((im(:)-im_hat(:)).^2))

Output:

MSE =5.8587e-27

As it is expected, the output is the same as Example 2 above. However, comparing the above piece of code with Example 2, you will notice that the wavedec2 stores coefficients in a completely different way.

The following example shows how wavedec2 stores coefficients.

Example 5: How “wavedec2” stores coefficients?

Consider a one level-decomposition via dwt2 :

[A1, H1, V1, D1] = dwt2(im,wname);

dwt2 returns all coefficients of a single-level DWT as it is expected (Fig. 1 above). Let’s use wavedec2:

N = 1; % Number of decomposition levels
[C,S] = wavedec2(im,N,wname);

Here, all coefficients are stored in a long row vector C. Then, a bookkeeping matrix (S) is used to store the number of coefficients by level and orientations.

In this example, the barbara image has 512x512 pixels. We expect that each coefficient matrix (namely, A1, H1, V1, and D1) has 256x256 pixels. Here, C has 1x262144. This means that wavedec2 puts all coefficients (256*256*4 = 262144) one after the other in a single row.

Now, let’s see the contents of S matrix:

S(1,:) shows size of approximation coefficients at scale N (here, N = 1) and the last row (here, S(3,:)) shows size of input image.

Let’s extract A1:

start_ind = 1;
last_ind = S(1,1)*S(1,2);
A1_vector = C(1,start_ind:last_ind);A1 = reshape(A1_vector,S(1,:));

But, how we can extract detail coefficients?. As you can guess in this simple example, the middle row (or rows) in S contains size of detail coefficients at other scales. Here, we just have one scale (N = 1). So, we just have one middle row. Let’s extract detail coefficients H1, V1, andD1:

% Note that "A1" is extracted above and "last_ind" refers to index of its last element.start_ind = last_ind + 1;
last_ind = last_ind + S(2,1)*S(2,2);
H1_vector = C(1,start_ind:last_ind);start_ind = last_ind + 1;
last_ind = last_ind + S(2,1)*S(2,2);
V1_vector = C(1,start_ind:last_ind);start_ind = last_ind + 1;
last_ind = last_ind + S(2,1)*S(2,2);
D1_vector = C(1,start_ind:last_ind);

Do not forget to reshape vectors into matrices:

% Reshape the extracted coefficient vector into matrixH1 = reshape(H1_vector,S(2,:));
V1 = reshape(V1_vector,S(2,:));
D1 = reshape(D1_vector,S(2,:));

From the above discussion, we can conclude that if we have separate coefficient matrices, we can easily create a vector similar to C:

C2 = [A1(:)',H1(:)',V1(:)',D1(:)'];isequal(C2,C)
ans =
logical
1

Now, as a simple exercise, perform a 2 level DWT decomposition and see how S can be used to extract coefficients at each scale. The coefficients are stored in C in this order: A2-H2-V2-D2-H1-V1-D1

Example 6: simple way of extracting coefficients returned by wavedec2

The above procedure of extracting coefficients by using vector C and bookkeeping matrix S might be cumbersome. Fortunately, MATLAB provides two functions (namely, appcoef2 and detcoef2) which can be used to extract coefficients easily:

% Perform a 2 level DWT
N = 2
wname = 'db1';
[C,S] = wavedec2(im,N,wname);% Extract the coefficients at level 1:
level = 1;
A1 = appcoef2(C,S,wname,level);
[H1,V1,D1] = detcoef2('all',C,S,level);

The function appcoef2 returns the approximation coefficients at the specified level (level can be 0, …, N, where the 0th level is just the input image). The function detcoef2 returns the detail coefficients at the specified level (level can be 1, … , N).

As a simple exercise, perform a two-level DWT decomposition using wavedec2 and try to extract V2 without detcoef2. You will find it somewhat boring.

Additional notes

1 — Whether you are using dwt2 recursively or wavedec2, there is a limit on the number of decomposition levels. This limit can be obtained by wmaxlev function:

sz = [512 128]; % input image size
wname = 'haar';
N = wmaxlev(sz,wname)

2 — Noisy data easily degrades reconstruction quality since wavelet’s are good at representing natural images, not noise!. This property is used to design denoising methods.

%
% Get input
%
im = imread('barbara.tif');
im = double(im);%
% Add noise
%
sigma = 30;
n = randn(size(im)) * sigma;
imn = im + n;%
% Compute DWT and recover the clean image ("im")
%
wname = 'dmey'; % Discrete Meyer wavelet filters are used.
[A1,H1,V1,D1] = dwt2(im,wname);
im_hat = idwt2(A1,H1,V1,D1,wname);MSE_im=mean(mean((im(:)-im_hat(:)).^2))%
% Compute DWT and recover the noisy image ("imn")
%
[A1n,H1n,V1n,D1n] = dwt2(imn,wname);
imn_hat = idwt2(A1n,H1n,V1n,D1n,wname);
MSE_imn=mean(mean((imn(:)-imn_hat(:)).^2))

The output is:

MSE_im =2.6043e-08MSE_imn =3.8584e-08

It is clear that the DWT is better at recovering clean image.