Fast Wavelet Transform (FWT) Algorithm

In 1988, Mallat produced a fast wavelet decomposition and reconstruction algorithm [1]. 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

For an orthogonal wavelet, in the multiresolution framework, 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):

$\frac{1}{2} ϕ (\frac{x}{2}) = \sum_{n \in Z} w_{n} ϕ (x - n)$

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

(w_n)_n∊Z

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 $\frac{1}{\sqrt{2}}$
Of norm 1
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) ^kLo_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.

Load the Daubechies’ extremal phase scaling filter and plot the coefficients.

load db6;
subplot(421); stem(db6,'markerfacecolor',[0 0 1]);
title('Original scaling filter');

Use orthfilt to return the analysis (decomposition) and synthesis (reconstruction) filters.

Obtain the discrete Fourier transforms (DFT) of the lowpass and highpass analysis filters. Plot the modulus of the DFT.

LoDFT = fft(Lo_D,64);
HiDFT = fft(Hi_D,64);
freq = -pi+(2*pi)/64:(2*pi)/64:pi;
subplot(427); plot(freq,fftshift(abs(LoDFT)));
set(gca,'xlim',[-pi,pi]); xlabel('Radians/sample');
title('DFT Modulus - Lowpass Filter')
subplot(428); plot(freq,fftshift(abs(HiDFT)));
set(gca,'xlim',[-pi,pi]); xlabel('Radians/sample');
title('Highpass Filter');

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 2L. The result of convolving a length N signal with a length 2L filter is N+2L–1. Therefore, the signals F and G are of length N + 2L – 1. After downsampling by 2, the coefficient vectors cA₁ and cD₁ are of length

$⌊ \frac{N - 1}{2} + L ⌋ .$

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 1-D wavelets by tensorial product.
This kind of 2-D 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 2-D 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 Border Effects.

Let us denote h = Lo_R and g = Hi_R and focus on the 1-D 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 (A_k^(j))k∊Z be the coordinates of the vector A_j:

$A_{j} = \sum_{k} A_{k}^{(j)} ϕ_{j, k}$

and A_k^(j+1) the coordinates of the vector A_j₊₁:

$A_{j + 1} = \sum_{k} A_{k}^{(j + 1)} ϕ_{j + 1, k}$

A_k^(j+1) is calculated using the formula

$A_{k}^{(j + 1)} = \sum_{n} h_{n - 2 k} A_{n}^{(j)}$

This formula resembles a convolution formula.

The computation is very simple.

Let us define

$\tilde{h} (k) = h (- k), and F_{k}^{(j + 1)} = \sum_{n} {\tilde{h}}_{k - n} A n^{(j)} .$

The sequence F^(j+1) is the filtered output of the sequence A^(j) by the filter $\tilde{h}$ .

We obtain

A_k^(j+1) = F_2k^(j+1)

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

The sequence A^(j+1) is the downsampled version of the sequence F^(j+1).

The initialization is carried out using A_k⁽⁰⁾ = s(k), 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.

The family $(ϕ_{0, k}, k \in Z)$ is formed of orthonormal functions. As a consequence for any j, the family $(ϕ_{j, k}, k \in Z)$ is orthonormal.
The double indexed family

$(ψ_{j, k}, j \in Z, k \in Z)$
is orthonormal.
For any j, the $(ϕ_{j, k}, k \in Z)$ are orthogonal to $(ψ_{j^{'}, k}, j^{'} \leq j, k \in Z)$ .
Between two successive scales, we have a fundamental relation, called the twin-scale relation.

Twin-Scale Relation for $ϕ$

$ϕ_{1, 0} = \sum_{k \in Z} h_{k} ϕ_{0, k}$

$ϕ_{j + 1, 0} = \sum_{k \in Z} h_{k} ϕ_{j, k}$

This relation introduces the algorithm's h filter $(h_{n} = \sqrt{2 ω_{n}})$ . For more information, see Filters Used to Calculate the DWT and IDWT.
We check that:
1. The coordinate of φ_j+1,0 on φ_j,kis h_k and does not depend on j.
2. The coordinate of φ_j+1,n on φ_j,k is equal to $〈 ϕ_{j + 1, n}, ϕ_{j, k} 〉 = h_{k - 2 n}$ .
These relations supply the ingredients for the algorithm.
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 $ϕ$

$ψ_{1, 0} = \sum_{k \in Z} g_{k} ϕ_{0, k}$

$ψ_{j + 1, 0} = \sum_{k \in Z} g_{k} ϕ_{j, k}$
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_j₁. The procedure is the same to calculate A = A_j₊₁ + D_j₊₁.
Let us define α_n, δ_n, $α_{k}^{0}$ by
$A_{1} = \sum_{n} a_{n} ϕ_{1, n} D_{1} = \sum_{n} δ_{n} ψ_{1, n} A_{0} = \sum_{k} a_{k}^{0} ϕ_{0, k}$
Let us assess the $α_{k}^{0}$ coordinates as
$\begin{array}{l} a_{k}^{0} = 〈 A_{0}, ϕ_{0, k} 〉 = 〈 A_{1} + D_{1}, ϕ_{0, k} 〉 = 〈 A_{1}, ϕ_{0, k} 〉 + 〈 D_{1}, ϕ_{0, k} 〉 \\ = \sum_{n} a_{n} 〈 ϕ_{1, n}, ϕ_{0, k} 〉 + \sum_{n} δ_{n} 〈 ψ_{1, n}, ϕ_{0, k} 〉 \\ = \sum_{n} a_{n} h_{k - 2 n} + \sum_{n} δ_{n} g_{k - 2 n} \end{array}$

Twin-Scale Relation for $ϕ$
$ϕ_{1, 0} = \sum_{k \in Z} h_{k} ϕ_{0, k}$	$ϕ_{j + 1, 0} = \sum_{k \in Z} h_{k} ϕ_{j, k}$

Twin-Scale Relation Between $ψ$ and $ϕ$
$ψ_{1, 0} = \sum_{k \in Z} g_{k} ϕ_{0, k}$	$ψ_{j + 1, 0} = \sum_{k \in Z} g_{k} ϕ_{j, k}$

We will focus our study on the first sum $\sum_{n} a_{n} h_{k - 2 n};$ the second sum $\sum_{n} δ_{n} g_{k - 2 n}$ 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)

$\begin{array}{l} \sum_{n} a_{n} h_{- 2 n} = \dots + a_{- 1} h_{2} + a_{0} h_{0} + a_{1} h_{- 2} + a_{2} h_{- 4} + \dots \\ = \dots + a_{- 1} h_{2} + 0 h_{1} + a_{0} h_{0} + 0 h_{- 1} + a_{1} h_{- 2} + 0 h_{- 3} + a_{2} h_{- 4} + \dots \end{array}$

If we transform the $(α_{n})$ sequence into a new sequence $({\tilde{α}}_{n})$ defined by

..., α_–1, 0, α₀, 0, α₁, 0, α₂, 0, ... that is precisely

${\tilde{a}}_{2 n} = a_{n}, {\tilde{a}}_{2 n} + 1 = 0$

Then

$\sum_{n} a_{n} h_{- 2 n} = \sum_{n} {\tilde{a}}_{n} h_{- n}$

and by extension

$\sum_{n} a_{n} h_{k - 2 n} = \sum_{n} {\tilde{a}}_{n} h_{k - n}$

Since

$a_{k}^{0} = \sum_{n} {\tilde{a}}_{n} h_{k - n} + \sum_{n} {\tilde{δ}}_{n} g_{k - n}$

the reconstruction steps are:

Replace the α and δ sequences by upsampled versions α˜ and $\tilde{δ}$ inserting zeros.
Filter by h and g respectively.
Sum the obtained sequences.

1-D Wavelet Capabilities

Basic 1-D Objects

	Objects	Description
Signal in original time	s A_k, 0 ≤ k ≤ j D_k, 1 ≤ k ≤ j	Original signal Approximation at level k Detail at level k
Coefficients in scale-related time	cA_k, 1 ≤ k ≤ j cD_k, 1 ≤ k ≤ j [cA_j, cD_j, ..., cD₁]	Approximation coefficients at level k Detail coefficients at level k Wavelet decomposition at level j, j ≥ 1

Objects

Description

Signal in original time

A_k, 0 ≤ k ≤ j

D_k, 1 ≤ k ≤ j

Original signal

Approximation at level k

Detail at level k

Coefficients in scale-related time

cA_k, 1 ≤ k ≤ j

cD_k, 1 ≤ k ≤ j

[cA_j, cD_j, ..., cD₁]

Approximation coefficients at level k

Detail coefficients at level k

Wavelet decomposition at level j, j ≥ 1

Analysis-Decomposition Capabilities

Purpose	Input	Output	File
Single-level decomposition	s	cA₁, cD₁	`dwt`
Single-level decomposition	cA_j	cA_j₊₁, cD_j₊₁	`dwt`
Decomposition	s	[cA_j, cD_j, ..., cD₁]	`wavedec`

Synthesis-Reconstruction Capabilities

Purpose	Input	Output	File
Single-level reconstruction	cA₁, cD₁	s or A₀	`idwt`
Single-level reconstruction	cA_j₊₁, cD_j₊₁	cA_j	`idwt`
Full reconstruction	[cA_j, cD_j, ..., cD₁]	s or A₀	`waverec`
Selective reconstruction	[cA_j, cD_j, ..., cD₁]	A_l, D_m	`wrcoef`

Decomposition Structure Utilities

Purpose	Input	Output	File
Extraction of detail coefficients	[cA_j, cD_j, ..., cD₁]	cD_k, 1 ≤ k ≤ j	`detcoef`
Extraction of approximation coefficients	[cA_j, cD_j, ..., cD₁]	cA_k, 0≤ k ≤ j	`appcoef`
Recomposition of the decomposition structure	[cA_j, cD_j, ..., cD₁]	[cA_k, cD_k, ..., cD₁] 1 ≤ k ≤ j	`upwlev`

To illustrate command-line mode for 1-D capabilities, see 1-D Decimated Wavelet Analysis.

2-D Wavelet Capabilities

Basic 2-D Objects

	Objects	Description
Image in original resolution	s	Original image
	A₀	Approximation at level 0
	A_k, 1 ≤ k ≤ j	Approximation at level k
	D_k, 1 ≤ k ≤ j	Details at level k
Coefficients in scale-related resolution	cA_k, 1 ≤ k ≤ j	Approximation coefficients at level k
	cD_k, 1 ≤ k ≤ j	Detail coefficients at level k
	[cA_j, cD_j, ..., cD₁]	Wavelet decomposition at level j

D_k stands for $[D_{k}^{(h)}, D_{k}^{(v)}, D_{k}^{(d)}]$ , the horizontal, vertical, and diagonal details at level k.

The same holds for cD_k, which stands for $[c D_{k}^{(h)}, c D_{k}^{(v)}, c D_{k}^{(d)}]$ .

The 2-D files are the same as those for the 1-D case, but with a 2 appended on the end of the command.

For example, idwt becomes idwt2. For more information, see 1-D Wavelet Capabilities.

To illustrate command-line mode for 2-D capabilities, see Wavelet Image Analysis and Compression..

References

[1] Mallat, S. G. “A Theory for Multiresolution Signal Decomposition: The Wavelet Representation,” IEEE Transactions on Pattern Analysis and Machine Intelligence. Vol. 11, Issue 7, July 1989, pp. 674–693.