Fast Wavelet Transform (FWT) Algorithm

In 1988, Mallat produced a fast wavelet decomposition and reconstruction algorithm [Mal89]. The Mallat algorithm for discrete wavelet transform (DWT) is, in fact, a classical scheme in the signal processing community, known as a two-channel subband coder using conjugate quadrature filters or quadrature mirror filters (QMFs).

• The decomposition algorithm starts with signal s, next calculates the coordinates of A₁ and D₁, and then those of A₂ and D₂, and so on.
• The reconstruction algorithm called the inverse discrete wavelet transform (IDWT) starts from the coordinates of A_J and D_J, next calculates the coordinates of A_J-1, and then using the coordinates of A_J-1 and D_J-1 calculates those of A_J-2, and so on.

This section addresses the following topics:

• Filters Used to Calculate the DWT and IDWT
• Algorithms
• Why Does Such an Algorithm Exist?
• One-Dimensional Wavelet Capabilities
• Two-Dimensional Wavelet Capabilities

Filters Used to Calculate the DWT and IDWT

For an orthogonal wavelet, in the multiresolution framework (see [Dau92] in Using Wavelet Packets), we start with the scaling function and the wavelet function . One of the fundamental relations is the twin-scale relation (dilation equation or refinement equation):

• 

All the filters used in DWT and IDWT are intimately related to the sequence

• 

Clearly if is compactly supported, the sequence (w_n) is finite and can be viewed as a filter. The filter W, which is called the scaling filter (nonnormalized), is

• Finite Impulse Response (FIR)
• Of length 2N
• Of sum 1
• Of norm
  
  
  
• A low-pass filter

For example, for the db3 scaling filter,

• load db3 
  db3
      db3 =
      0.2352   0.5706   0.3252  -0.0955  -0.0604   0.0249
  
  sum(db3)
      ans =
           1.0000
  
  norm(db3)
      ans =
           0.7071
  

From filter W, we define four FIR filters, of length 2N and of norm 1, organized as follows.

Filters	Low-Pass	High-Pass
Decomposition	Lo_D	Hi_D
Reconstruction	Lo_R	Hi_R

The four filters are computed using the following scheme.

where qmf is such that Hi_R and Lo_R are quadrature mirror filters (i.e., Hi_R(k) = (-1) k Lo_R(2N + 1 - k)) for k = 1, 2, ..., 2N.

Note that wrev flips the filter coefficients. So Hi_D and Lo_D are also quadrature mirror filters. The computation of these filters is performed using orthfilt. Next, we illustrate these properties with the db6 wavelet. The plots associated with the following commands are shown in Figure 6-8.

• % Load scaling filter.
  load db6; w = db6; 
  subplot(421); stem(w); title('Original scaling filter');
  
  % Compute the four filters.
  [Lo_D,Hi_D,Lo_R,Hi_R] = orthfilt(w); 
  subplot(423); stem(Lo_D); 
  title('Decomposition low-pass filter Lo{\_}D'); 
  subplot(424); stem(Hi_D); 
  title('Decomposition high-pass filter Hi{\_}D'); 
  subplot(425); stem(Lo_R); 
  title('Reconstruction low-pass filter Lo{\_}R'); 
  subplot(426); stem(Hi_R); 
  title('Reconstruction high-pass filter Hi{\_}R');
  
  % High and low frequency illustration.
  n = length(Hi_D);
  freqfft = (0:n-1)/n;
  nn = 1:n;
  N = 10*n;
  for k=1:N
      lambda(k) = (k-1)/N;
      XLo_D(k) = exp(-2*pi*j*lambda(k)*(nn-1))*Lo_D';
      XHi_D(k) = exp(-2*pi*j*lambda(k)*(nn-1))*Hi_D';
  end
  fftld = fft(Lo_D);
  ffthd = fft(Hi_D);
  subplot(427); plot(lambda,abs(XLo_D),freqfft,abs(fftld),'o'); 
  title('Transfer modulus: lowpass (Lo{\_}D or Lo{\_}R') 
  subplot(428); plot(lambda,abs(XHi_D),freqfft,abs(ffthd),'o'); 
  title('Transfer modulus: highpass (Hi{\_}D or Hi{\_}R') 
  
  
  

Figure 6-8: Four Wavelet Filters for db6

Algorithms

Given a signal s of length N, the DWT consists of log₂N stages at most. Starting from s, the first step produces two sets of coefficients: approximation coefficients cA₁, and detail coefficients cD₁. These vectors are obtained by convolving s with the low-pass filter Lo_D for approximation, and with the high-pass filter Hi_D for detail, followed by dyadic decimation.

More precisely, the first step is

The length of each filter is equal to 2n. If N= length (s), the signals F and G are of length N+ 2n - 1, and then the coefficients cA₁ and cD₁ are of length

• 

The next step splits the approximation coefficients cA₁ in two parts using the same scheme, replacing s by cA₁ and producing cA₂ and cD₂, and so on.

So the wavelet decomposition of the signal s analyzed at level j has the following structure: [cA_j, cD_j, ..., cD₁].

This structure contains for J = 3 the terminal nodes of the following tree.

• Conversely, starting from cA_j and cD_j, the IDWT reconstructs cA_j-1, inverting the decomposition step by inserting zeros and convolving the results with the reconstruction filters.

  

• For images, a similar algorithm is possible for two-dimensional wavelets and scaling functions obtained from one-dimensional wavelets by tensorial product.

• This kind of two-dimensional DWT leads to a decomposition of approximation coefficients at level j in four components: the approximation at level j + 1 and the details in three orientations (horizontal, vertical, and diagonal).

  The following charts describe the basic decomposition and reconstruction steps for images.

  

  

So, for J = 2, the two-dimensional wavelet tree has the following form.

Finally, let us mention that, for biorthogonal wavelets, the same algorithms hold but the decomposition filters on one hand and the reconstruction filters on the other hand are obtained from two distinct scaling functions associated with two multiresolution analyses in duality.

In this case, the filters for decomposition and reconstruction are, in general, of different odd lengths. This situation occurs, for example, for "splines" biorthogonal wavelets used in the toolbox. By zero-padding, the four filters can be extended in such a way that they will have the same even length.

Why Does Such an Algorithm Exist?

The previous paragraph describes algorithms designed for finite-length signals or images. To understand the rationale, we must consider infinite-length signals. The methods for the extension of a given finite-length signal are described in Dealing with Border Distortion.

Let us denote h = Lo_R and g = Hi_R and focus on the one-dimensional case.

We first justify how to go from level j to level j+1, for the approximation vector. This is the main step of the decomposition algorithm for the computation of the approximations. The details are calculated in the same way using the filter g instead of filter h.

Let be the coordinates of the vector A_j:

• 

and the coordinates of the vector A_j+1:

• 

is calculated using the formula

• 

This formula resembles a convolution formula.

The computation is very simple.

Let us define

• 

The sequence is the filtered output of the sequence by the filter .

We obtain

• 

We have to take the even index values of F. This is downsampling.

The sequence is the downsampled version of the sequence .

The initialization is carried out using , where s(k) is the signal value at time k.

There are several reasons for this surprising result, all of which are linked to the multiresolution situation and to a few of the properties of the functions _j,k and _j,k.

Let us now describe some of them.

1. The family is formed of orthonormal functions. As a consequence for any j, the family is orthonormal.
2. The double indexed family is orthonormal.
3. For any j, the are orthogonal to .
4. Between two successive scales, we have a fundamental relation, called the twin-scale relation.
   
   

   Twin-Scale Relation for

This relation introduces the algorithm's h filter (). For more information, see Filters Used to Calculate the DWT and IDWT.

5. We check that:
1. The coordinate of on _j,k is and does not depend on j.
2. The coordinate of on _j,k is equal to .
6. These relations supply the ingredients for the algorithm.
7. Up to now we used the filter h. The high-pass filter g is used in the twin scales relation linking the and functions. Between two successive scales, we have the following twin-scale fundamental relation.
   
   

   Twin-Scale Relation Between and

8. After the decomposition step, we justify now the reconstruction algorithm by building it. Let us simplify the notation. Starting from A₁ and D₁, let us study A₀ = A₁ + D₁. The procedure is the same to calculate A_j = A_j+1 + D_j+1.

1. Let us define _n, _n, by

• 

1. Let us assess the coordinates as

• 

• 

• 

We will focus our study on the first sum ; the second sum

is handled in a similar manner.

The calculations are easily organized if we note that (taking k = 0 in the previous formulas, makes things simpler)

• 

• 

If we transform the (_n) sequence into a new sequence defined by

      ..., _-1, 0, ₀, 0, ₁, 0, ₂, 0, ... that is precisely

• 

Then

• 

and by extension

• 

Since

• 

the reconstruction steps are:

1. Replace the and sequences by upsampled versions and inserting zeros.
2. Filter by h and g respectively.
3. Sum the obtained sequences.

One-Dimensional Wavelet Capabilities

Basic One-Dimensional Objects.   

	Objects	Description
Signal in original time	s Ak, 0 k j Dk, 1 k j	Original signal Approximation at level k Detail at level k
Coefficients in scale-related time	cAk, 1 k j cDk, 1 k j [cAj, cDj, ..., cD1]	Approximation coefficients at level k Detail coefficients at level k Wavelet decomposition at level j, j 1

Analysis-Decomposition Capabilities.   

Purpose	Input	Output	M-File
Single-level decomposition	s	cA1, cD1	dwt
Single-level decomposition	cAj	cAj+1, cDj+1	dwt
Decomposition	s	[cAj, cDj, ..., cD1]	wavedec

Synthesis-Reconstruction Capabilities.   

Purpose	Input	Output	M-File
Single-level reconstruction	cA1, cD1	s or A0	idwt
Single-level reconstruction	cAj+1, cDj+1	cAj	idwt
Full reconstruction	[cAj, cDj, ..., cD1]	s or A0	waverec
Selective reconstruction	[cAj, cDj, ..., cD1]	Al, Dm	wrcoef

Decomposition Structure Utilities.   .

Purpose	Input	Output	M-File
Extraction of detail coefficients	[cAj, cDj, ..., cD1]	cDk, 1 k j	detcoef
Extraction of approximation coefficients	[cAj, cDj, ..., cD1]	cAk, 0 k j	appcoef
Recomposition of the decomposition structure	[cAj, cDj, ..., cD1]	[cAk, cDk, ..., cD1] 1 k j	upwlev

To illustrate command-line mode for one-dimensional capabilities, see One-Dimensional Analysis Using the Command Line.

Two-Dimensional Wavelet Capabilities

Basic Two-Dimensional Objects.   

	Objects	Description
Image in original resolution	s	Original image
	A0	Approximation at level 0
	Ak, 1 k j	Approximation at level k
	Dk, 1 k j	Details at level k
Coefficients in scale-related resolution	cAk, 1 k j	Approximation coefficients at level k
	cDk, 1 k j	Detail coefficients at level k
	[cAj, cDj, ..., cD1]	Wavelet decomposition at level j

D_k stands for , the horizontal, vertical, and diagonal details at level k.

The same holds for cD_k, which stands for .

The two-dimensional M-files are the same as those for the one-dimensional case, but with a 2 appended on the end of the command.

For example, idwt becomes idwt2. For more information, see One-Dimensional Wavelet Capabilities.

To illustrate command-line mode for two-dimensional capabilities, see Two-Dimensional Analysis Using the Command Line.

Provide feedback about this page

General Concepts

Dealing with Border Distortion

© 1984-2009- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS