4. Wavelets

Chapter PDF · Complete book

The reference for this chapter is [23].

4.1 Multiresolution Approximation Spaces¶

A multiresolution approximation of $L^2(\mathbb{R})$ consists of nested closed subspaces $(V_j)_j$

related by dyadic scaling and invariant under translations on the corresponding dyadic grid:

Large $j$ corresponds to coarse approximation spaces, and $2^j$ is often called the “scale”.

At the fine-scale limit on the left of (4.1), $\cup_j V_j$ is dense in $L^2(\mathbb{R})$. Equivalently, $P_{V_j}(f) \rightarrow f$ as $j \rightarrow -\infty$, where $P_V$ denotes orthogonal projection onto $V$:

The limit on the right of (4.1) means that $\cap_j V_j = \{0\}$, or equivalently that $P_{V_j}(f) \rightarrow 0$ as $j \rightarrow +\infty$.

The first example consists of piecewise constant functions on dyadic intervals

A second example is the space used for Shannon interpolation of bandlimited signals

whose elements are recovered exactly from the samples $f(2^j n)$. As for piecewise constant signals, the sampling grid has spacing $2^j$. Here, however, orthogonal projection is a bandlimiting operation:

Scaling functions.

We also require a scaling function $\varphi\in L^2(\mathbb{R})$ such that

By the dilation property, this implies that

The normalization ensures $|\!| \varphi_{j,n} |\!| = 1$.

Note that one then has

Figure 4.1 illustrates the effects of translation and scaling.

For the case of piecewise constant signals (4.2), one can use

For the case of Shannon multiresolution (4.3), one can use $\varphi(t) = \sin(\pi t)/(\pi t)$ and one verifies

Spectral orthogonalization.

In many applications, $V_0$ is generated by translates of a function that are not orthogonal: $V_0=\overline{\mathop{\mathrm{Span}}}\{\theta(\cdot-n)\}_{n\in\mathbb{Z}}$. The next proposition orthogonalizes this family in the Fourier domain.

Figure 4.3 shows the orthogonal cardinal spline functions obtained by this construction.

Spline interpolation provides a typical application. For example, cubic splines are generated by translates of a box-spline function $\theta$, which is a compactly supported piecewise polynomial.

4.2 Multiresolution Detail Spaces¶

The detail spaces are the orthogonal complements associated with the inclusions $V_j \subset V_{j-1}$. Since these subspaces are closed,

The following diagram shows these nested approximation spaces and their orthogonal complements:

Suppose that $W_0$ admits an orthonormal basis of translates $\{\psi(\cdot-n)\}_n$. Scaling then gives

Orthogonal decomposition across all scales yields

For every $f \in L^2(\mathbb{R})$, the corresponding expansion converges in $L^2(\mathbb{R})$:

This decomposition shows that

is an orthonormal basis of $L^2(\mathbb{R})$, which is called a wavelet basis.

Stopping at a coarsest scale gives the orthonormal basis

The forward wavelet transform computes all inner products of a function $f$ with the elements of this basis.

Haar wavelets.

For the Haar multiresolution (4.2), one has

One choice of mother wavelet is

as shown in Figure 4.2.

Figure 4.5 illustrates how subtracting consecutive approximation-space projections isolates the details at each scale.

Shannon and splines.

Figure 4.6 compares the sharp frequency partitions of Shannon wavelets with the smooth transitions of spline wavelets. Shannon wavelets can be viewed as a limiting case as the spline degree increases.

4.3 On Bounded Domains¶

On the periodic domain $\mathbb{T}=\mathbb{R}/\mathbb{Z}$, we use period 1 rather than the $2\pi$-period convention of the Fourier chapter. To obtain an orthonormal wavelet basis of $L^2(\mathbb{T})$, periodize the wavelets and restrict their translation points $2^j n$ to $[0,1]$, choosing $0 \leqslant n < 2^{-j}$. As in (1.6), the periodization of $f \in L^1(\mathbb{R})$ is

For an integer coarsest level $j_0\leqslant 0$, this gives the family

which is an orthonormal basis of $L^2(\mathbb{T})$, illustrated in Figure 4.7.

One can also construct wavelet bases using Neumann (mirror) boundary conditions, but this is more involved.

4.4 Fast Wavelet Transform¶

4.4.1 Discretization¶

We now work over $\mathbb{R}/\mathbb{Z}$.

We assume access to a discrete signal $a_J \in \mathbb{R}^N$ with $N = 2^{-J}$ at a fixed scale $2^J$, whose entries are inner products with the scaling functions, i.e.

for the underlying continuous function $f$. The approximation uses the normalization $\int\varphi=1$. The equality assumes that acquisition gives exact access to $P_{V_J}f$, much as Shannon sampling assumes bandlimiting. This model is reasonable when the scaling functions $(\varphi_{J,n})_n$ approximate the acquisition device’s point-spread function. Preprocessing the measurements can sometimes improve agreement with (4.5).

Starting from $a_J$, the discrete wavelet transform computes the coefficients

in order of increasing $j$. After an approximation vector $a_j$ has been used to compute the next scale, it can be discarded; the detail vectors $d_j$ are retained.

The forward discrete wavelet transform on a bounded domain is thus the orthogonal finite-dimensional map

By orthogonality, the inverse transform is the adjoint.

Figure 4.8 shows examples of wavelet coefficients. For each scale $2^j$, there are $2^{-j}$ coefficients.

4.4.2 Forward Fast Wavelet Transform (FWT)¶

The algorithm applies a sequence of elementary operators

There are $|J|=\log_2(N)$ steps. Each $\mathcal{W}_j$ changes coordinates between orthonormal bases and is therefore orthogonal.

To describe the algorithm for $\mathcal{W}_j$, we introduce the filter coefficients $h,g \in \mathbb{R}^\mathbb{Z}$

Figure 4.9 illustrates the computation of these weights for the Haar system.

We denote by $\downarrow_2 : \mathbb{R}^K \rightarrow \mathbb{R}^{K/2}$ the subsampling operator by a factor of 2, i.e.

Throughout the following derivation, the scaling function, wavelet, and filters are real, and the filters decay sufficiently fast. The notation $\bar h_n=h_{-n}$ and $\bar g_n=g_{-n}$ denotes reflection. Complex filters would require conjugation as well.

Figure 4.10 shows two successive decomposition steps, each separating a coarse approximation from its details.

The FWT thus operates as follows. The same real filters also act on complex-valued input signals:

• Input: signal $f \in \mathbb{C}^N$.

• Initialization: $a_J = f$.

• For $j=J,\ldots,j_0-1$.

  $$a_{j+1} = (a_j \star \bar h) \downarrow 2 \quad \text{and} \quad d_{j+1} = (a_j \star \bar g) \downarrow 2$$

• Output: the coefficients $\{d_j\}_{J<j\leqslant j_0} \cup \{ a_{j_0} \}$.

If the supports of $h$ and $g$ each contain at most $C$ indices, then evaluating each $\mathcal{W}_j$ requires $(2C) 2^{-j}$ operations, so that the complexity of the whole wavelet transform is

This shows that the fast wavelet transform is a linear-time algorithm.

Figure 4.11 shows the iterative extraction of the wavelet coefficients.

Figure 4.12 illustrates the storage layout: at each step, the new vectors $a_j$ and $d_j$ occupy the left portion of the output array.

Fast Haar transform.

For the Haar wavelets, one has

This corresponds to the filters

The Haar transform iterates normalized sums and differences:

• Input: signal $f \in \mathbb{C}^N$.

• Initialization: $a_J = f$.

• For $j=J,\ldots,j_0-1$.

  $$a_{j+1,n} = \frac{1}{\sqrt{2}} (a_{j,2n}+a_{j,2n+1}) \quad \text{and} \quad d_{j+1,n} = \frac{1}{\sqrt{2}} (a_{j,2n}-a_{j,2n+1}).$$

• Output: the coefficients $\{d_j\}_{J<j\leqslant j_0} \cup \{ a_{j_0} \}$.

4.4.3 Inverse Fast Transform (iFWT)¶

The inverse algorithm proceeds by inverting each step (4.6)

where $\mathcal{W}_j^{*}$ is the adjoint for the canonical inner product on $\mathbb{R}^{2^{-j+1}}$; in matrix form, it is the transpose.

Let $\uparrow_2 : \mathbb{R}^{K/2} \rightarrow \mathbb{R}^K$ denote up-sampling by inserting zeros:

Figure 4.14 makes this transpose relation concrete for one Haar step.

The inverse fast wavelet transform iteratively applies this elementary step

• Input: $\{d_j\}_{J<j\leqslant j_0} \cup \{ a_{j_0} \}$.

• For $j=j_0,j_0-1,\ldots,J+1$.

  $$a_{j-1} = (a_j \uparrow 2) \star h + (d_j \uparrow 2) \star g.$$

• Output: $f = a_J$.

The block diagram in Figure 4.15 reverses the operations of the forward transform in Figure 4.10.

4.5 2-D Wavelets¶

4.5.1 Anisotropic Wavelets¶

The anisotropic basis on $\mathbb{T}^2$ is the tensor product of the complete one-dimensional periodized basis $\mathcal B$ (including its coarsest scaling functions):

Its wavelet–wavelet elements have the form

The basis also includes wavelet–scaling, scaling–wavelet, and scaling–scaling products. The anisotropic 2-D transform follows the same separable construction as the 2-D FFT in Section 2.5.2. Treat the image $a_J \in \mathbb{R}^{2^{-J} \times 2^{-J}}$ as a matrix, apply the 1-D FWT to each row, and then apply it to each column. The total cost is $O(N)$ for $N=2^{-2J}$ pixels.

4.5.2 Isotropic Wavelets¶

The anisotropic atoms (4.12) have independent scales in the two coordinate directions. For compactly supported wavelets, their supports lie in axis-aligned rectangles of size proportional to $2^{j_1}\times2^{j_2}$. This directional bias can produce visible artifacts in denoising and compression when it does not match the image geometry.

An alternative uses a common scale in both directions. Tensor products of one-dimensional approximation spaces give a two-dimensional multiresolution

We write

for this isotropic 2-D approximation space.

Recall that the tensor product of two closed subspaces $V_1,V_2\subset L^2(\mathbb{R})$ is

If $(\varphi_k^s)_k$ are orthonormal bases of $V_s$, their tensor products $(\varphi_{k_1}^1(x_1)\varphi_{k_2}^2(x_2))_{k_1,k_2}$ form an orthonormal basis of $V_1 \otimes V_2$.

The tensor product distributes over these orthogonal sums

The labels $\{V,H,D\}$ stand for Vertical, Horizontal, Diagonal detail spaces.

The following diagram summarizes the nested spaces and their detail components:

For $j \in \mathbb{Z}$, translates of a mother wavelet span each detail space. We index them by $n=(n_1,n_2)\in\mathbb{Z}^2$, or by indices in $\{0,\ldots,2^{-j}-1\}^2$ on the torus $\mathbb{T}^2$:

where

and where the three mother wavelets are

Here the span is closed in $L^2$ on an unbounded domain; on the torus, the wavelets are periodized. Figure 4.18 displays examples of these wavelets.

Haar 2-D multiresolution.

The 2-D Haar construction gives piecewise constant approximations: functions in $V_j \otimes V_j$ are constant on squares of size $2^j \times 2^j$. Figure 4.19 shows the corresponding projections of an image.

Discrete 2-D wavelet coefficients.

As in (4.5), assume that acquisition provides inner products of the continuous signal $f$ with scaling functions at scale $2^J$, with $N=2^{-J}$ samples per direction:

Discrete wavelet coefficients are defined as

Here the wavelets are periodized.

Approximation coefficients are defined as

Figure 4.20 shows how the wavelet coefficients are arranged in an image of $N^2$ pixels.

Figure 4.21 shows other examples of wavelet decompositions.

One step of the forward 2-D transform.

Each step separates fine-scale approximation coefficients into three detail arrays and a coarse approximation:

As in one dimension, this map is orthogonal. It is implemented by applying the filtering and sub-sampling formula (4.8) in each coordinate direction.

One first applies 1-D horizontal filtering and sub-sampling

where $\star^H$ denotes convolution along the first coordinate (called horizontal here), applied to each column of the array

Here $a\in\mathbb{C}^{P\times P}$ is a matrix and $b\in\mathbb{C}^P$ is a filter. The notation $\downarrow^H 2$ denotes sub-sampling in the horizontal direction

One then applies 1-D vertical filtering and sub-sampling to $\tilde a_j$ and $\tilde d_j$ to obtain

where the vertical operators are defined analogously to the horizontal operators, acting on rows.

These two forward steps are shown in the block diagram in Figure 4.22. The updates can be performed in place, storing all coefficients in an image of $N^2$ pixels, as shown in Figure 4.23. This produces the familiar multiscale arrangement used in Figure 4.21.

Fast 2-D wavelet transform.

The 2-D FWT algorithm iterates these steps through the scales:

• Input: image $f \in \mathbb{C}^{N\times N}$.

• Initialization: $a_J = f$.

• For $j=J+1,\ldots,j_0$.

  $$\begin{array}{ll} \tilde a_j &= ( a_{j-1} \star^H \bar h ) \downarrow^H 2, \\ \tilde d_j &= (a_{j-1}\star^H\bar g) \downarrow^H 2, \\ a_j &= ( \tilde a_j \star^V \bar h ) \downarrow^V 2, \end{array}\qquad \begin{array}{ll} d_j^V &= ( \tilde a_j \star^V \bar g ) \downarrow^V 2,\\ d_j^H &= ( \tilde d_j \star^V \bar h ) \downarrow^V 2,\\ d_j^D &= ( \tilde d_j \star^V \bar g ) \downarrow^V 2. \end{array}$$

• Output: the coefficients $\{d_j^\omega\}_{J<j\leqslant j_0,\omega} \cup \{ a_{j_0} \}$.

Fast 2-D inverse wavelet transform.

The inverse transform undoes the horizontal and vertical filtering steps. The first step computes

where the vertical up-sampling is

The second inverse step computes

Figure 4.24 gives the inverse filter-bank diagram corresponding to Figure 4.22.

The inverse fast wavelet transform iteratively applies these elementary steps

• Input: $\{d_j^\omega\}_{J<j\leqslant j_0,\omega} \cup \{ a_{j_0} \}$.

• For $j=j_0,j_0-1,\ldots,J+1$.

  $$\begin{aligned} \tilde a_j &= (a_j\uparrow^V2)\star^Vh + (d_j^V\uparrow^V2)\star^Vg,\\ \tilde d_j &= (d_j^H\uparrow^V2)\star^Vh + (d_j^D\uparrow^V2)\star^Vg,\\ a_{j-1} &= (\tilde a_j\uparrow^H2)\star^Hh + (\tilde d_j\uparrow^H2)\star^Hg. \end{aligned}$$

• Output: $f = a_J$.

4.6 Wavelet Design¶

Implementing the FWT requires the filters $h$ and $g$, rather than explicit formulas for the scaling function $\varphi$ and wavelet $\psi$. Most useful wavelets have no elementary closed form; their refinement equations define them implicitly, and the cascade algorithm approximates them.

This section derives constraints on the filters and gives examples. Under the usual quadrature-mirror construction, the low-pass filter $h$ determines the high-pass filter $g$.

4.6.1 Low-pass Filter Constraints¶

We introduce the following three conditions on a filter $h$

These conditions use the Fourier series of the filter $h \in \mathbb{R}^\mathbb{Z}$:

If $h \in \ell^1(\mathbb{Z})$, this defines a continuous periodic function $\hat h \in \mathcal{C}^0(\mathbb{R}/2\pi\mathbb{Z})$, and the definition extends to $h \in \ell^2(\mathbb{Z})$ to give $\hat h \in L^2(\mathbb{R}/2\pi\mathbb{Z})$.

Condition ($C^*$) prevents zeros of $\hat h$ on the central half-period. It is a sufficient nondegeneracy condition for the converse construction used here.

4.6.2 High-pass Filter Constraints¶

We now introduce the following two conditions on a pair of filters $(g,h)$

Figure 4.25 illustrates these filter constraints.

Quadrature mirror filters.

Quadrature mirror filters (QMF) define $g$ as a function of $h$ so that the conditions of Theorem 4.5 are automatically satisfied. Most wavelet constructions implicitly use this choice (other phase choices can yield different wavelets spanning the same detail spaces).

4.6.3 Wavelet Design Constraints¶

One can construct a multiresolution wavelet pair by designing a finite filter $h$ satisfying ($C_1$)–($C^*$). The scaling function is obtained through the infinite product (4.16). Except in special cases such as Haar wavelets, it has no elementary closed form. The QMF relation (4.19) then defines $g$, and (4.17) defines $\psi$.

Fourier atoms are fixed exponentials, whereas mother wavelets can be designed to meet several requirements:

• Size of the support.

• Cancellation of polynomials, measured by the number $p$ of vanishing moments.

• Symmetry (for compactly supported scalar orthogonal wavelets, the Haar system is the exceptional symmetric case).

• Smoothness (number of derivatives).

We now explain how these design requirements interact with conditions ($C_1$)–($C_4$).

Vanishing moments.

A wavelet $\psi$ has $p$ vanishing moments if

Vanishing moments make $\langle f,\,\psi_{j,n}\rangle$ small when $f$ has regularity $\text{C}^\alpha$ with $\alpha<p$ on Supp$(\psi_{j,n})$. A Taylor expansion of $f$ about $2^j n$ on the wavelet’s support explains this cancellation.

The vanishing-moment condition has the following equivalent Fourier characterization; see Figure 4.26.

Conditions ($C_1$) and ($C_2$) give $\hat h(\pi)=0$, so every admissible wavelet has at least 1 vanishing moment: $\int \psi=0$. By (4.22), additional vanishing moments require the Fourier transform $\hat h$ to vanish to higher order at $\omega=\pi$.

Support.

Figure 4.27 shows the wavelet coefficients of a piecewise smooth signal. The cancellation provided by $\psi$ suppresses coefficients in smooth regions, leaving the largest coefficients near singularities.

Choosing $\psi$ with short support limits how many wavelets intersect each singularity.

One can show that the size of the support of $\psi$ is proportional to the size of the support of $h$.

Short support competes with the vanishing-moment requirement (4.21). Indeed, one can prove that for an orthogonal wavelet basis with $p$ vanishing moments

where the left-hand side denotes the length of the smallest closed interval containing the support of a compactly supported orthogonal mother wavelet.

Chapter 5 studies in detail the tradeoff between support size and vanishing moments in nonlinear approximation of piecewise smooth signals.

Smoothness.

In compression or denoising applications, an approximate signal is recovered from a partial set $I_M$ of coefficients,

This finite sum is at least as smooth as $\psi$.

A smooth wavelet $\psi$ helps avoid visible reconstruction artifacts. Smoothness controls the regularity of the reconstructed signal, while vanishing moments and localization primarily determine the approximation rates studied below. These properties are often linked: in many families, including the Daubechies wavelets introduced next, more vanishing moments also bring greater smoothness.

4.6.4 Daubechies Wavelets¶

To build a wavelet $\psi$ with a fixed number $p$ of vanishing moments, one designs the filter $h$, and uses the quadrature mirror filter relation (4.19) to compute $g$. One therefore looks for $h$ such that

These conditions give algebraic relations between the coefficients of $h$, which can be solved explicitly using the Euclidean division algorithm for polynomials.

The resulting Daubechies wavelets have $p$ vanishing moments and attain the minimum support length $2p-1$ among compactly supported orthogonal wavelets.

The following filter coefficients are rounded to four decimal places. For $p=1$, the construction gives the Haar system (up to the overall sign of the wavelet), with

For $p=2$, one obtains the celebrated Daubechies 4 filter

and for $p=3$,

Wavelet display.

Figure 4.28 compares Daubechies mother wavelets with increasing numbers of vanishing moments. Each continuous wavelet is approximated by a discrete wavelet $\bar \psi_{j,n}$ computed on a fine grid of $N$ samples. Applying the inverse transform to the following coefficient vector produces that discrete wavelet: $d_{j',n'} = \delta_{j-j'} \delta_{n-n'}$.