Wavelet Transform

The wavelet transform rewrites an image as a coarse, low-resolution approximation plus a set of detail coefficients, each of which records how the image changes at one scale, in one orientation, near one position. The scales are spaced by octaves, as in an image pyramid, but the discrete wavelet transform (DWT) in its orthogonal form keeps exactly as many coefficients as there are pixels and loses nothing: the image is recovered by the inverse transform. Because natural images are smooth almost everywhere and change abruptly only at edges and in texture, most detail coefficients are close to zero, and a small fraction of them describes most of the image. That sparsity is why wavelets are used to remove noise, to compress images and to analyze texture. The construction is also called the multiresolution decomposition or wavelet representation (Mallat’s term), the Mallat decomposition (PyWavelets’ documentation), or a subband decomposition.

Definition

A wavelet is a function $\psi$ of finite energy whose dilations and translations, $2^{-j/2}\psi(2^{-j}x - k)$ in one common indexing, form a basis in which signals are expanded. The discrete wavelet transform of a sampled signal computes the coefficients of such an expansion with a pyramid of filters: a low-pass filter and a high-pass filter split the signal into an approximation and a detail, each subsampled by two, and the split is repeated on the approximation. After $J$ levels the signal is represented by one approximation at the coarsest scale and $J$ sets of details. In two dimensions the filters are applied along rows and columns, and each level produces three detail images, one per orientation.

Intuition

Seen through a filter bank, an image is the sum of two parts: what survives averaging over neighbouring pairs of pixels, and what the averaging removed. With the simplest wavelet, Haar’s, the approximation of two pixels is their scaled sum and the detail is their scaled difference. In a flat region the difference is zero; at an edge it is large. Repeating the split on the averages gives differences between pairs of averages, which respond to edges and blobs twice as large, and so on up the octaves.

The image pyramid does something similar: its Laplacian bands also store, level by level, what smoothing removed. Mallat (1989) set out two shortcomings of that representation: the data at separate levels are correlated, so it is hard to tell whether a similarity across scales belongs to the image or to the redundancy of the representation, and the bands have no orientation selectivity. The wavelet representation removes both. Each band holds only what the next coarser level cannot represent, in an orthogonal complement, so the total number of coefficients does not grow; and in two dimensions the details are split into horizontal, vertical and diagonal components.

The other comparison is with the Fourier transform. A Fourier coefficient measures a grating that extends over the whole image, so a single sharp edge spreads into coefficients at many frequencies. A wavelet coefficient measures a band of frequencies over a small region, so an edge affects only the coefficients whose support it crosses, a few per scale. Mallat describes the wavelet decomposition as intermediate between spatial and Fourier representations.

Formal Definition

Multiresolution approximation

Mallat (1989) builds the transform from the operator $A_{2^j}$ that approximates a function $f \in L^2(\mathbb{R})$ at resolution $2^j$. He asks that it be a linear orthogonal projection onto a space $V_{2^j}$; that a finer approximation contain all the information of a coarser one, $V_{2^j} \subset V_{2^{j+1}}$; that the spaces be dilations of one another, $f(x) \in V_{2^j} \Leftrightarrow f(2x) \in V_{2^{j+1}}$; that the approximation be characterized by $2^j$ samples per unit length and translate with the signal; and that the approximations converge to $f$ as the resolution grows and to zero as it falls. A sequence of spaces with these properties is a multiresolution approximation. His first example is the space of functions constant on each interval $[k2^{-j}, (k+1)2^{-j})$, which leads to the Haar wavelet.

Mallat’s Theorem 1 states that each multiresolution approximation has a unique scaling function $\phi$ such that the functions $\sqrt{2^{-j}}\,\phi_{2^j}(x - 2^{-j}n)$, with $\phi_{2^j}(x) = 2^j\phi(2^j x)$, form an orthonormal basis of $V_{2^j}$. The scaling function is a low-pass filter. It is tied to a discrete filter $h[n] = \langle \phi_{2^{-1}}(u), \phi(u - n) \rangle$ whose transfer function $H(\omega) = \sum_n h[n]e^{-in\omega}$ satisfies (Theorem 2, eqs. 17a–17b)

$$|H(0)| = 1, \qquad |H(\omega)|^2 + |H(\omega + \pi)|^2 = 1 .$$

Mallat calls filters with the second property conjugate filters.

The wavelet

The detail at resolution $2^j$ is what $V_{2^{j+1}}$ contains and $V_{2^j}$ does not: the projection of $f$ onto the orthogonal complement $O_{2^j}$ of $V_{2^j}$ in $V_{2^{j+1}}$. Mallat’s Theorem 3 builds an orthonormal basis of $O_{2^j}$ from one function, the orthogonal wavelet, defined through its Fourier transform

$$\hat{\psi}(\omega) = G\!\left(\tfrac{\omega}{2}\right)\hat{\phi}\!\left(\tfrac{\omega}{2}\right), \qquad G(\omega) = e^{-i\omega}\,\overline{H(\omega + \pi)} ,$$

so that $\sqrt{2^{-j}}\,\psi_{2^j}(x - 2^{-j}n)$, over all $n$, is an orthonormal basis of $O_{2^j}$, and over all $j$ and $n$ an orthonormal basis of $L^2(\mathbb{R})$. In the time domain (eq. 29) the high-pass filter is the mirror of the low-pass one,

$$g[n] = (-1)^{1-n}\, h[1 - n],$$

and $H$ and $G$ form a pair of quadrature mirror filters. For the piecewise-constant example the wavelet is Haar’s: $\psi(x) = 1$ for $0 \le x < \tfrac12$, $-1$ for $\tfrac12 \le x < 1$, and 0 otherwise.

The pyramid algorithm

Since $\phi$ and $\psi$ at one resolution are combinations of $\phi$ at the next finer one, the inner products at a coarse level can be computed from those at the finer level by discrete filtering. Write $a_\ell$ for the approximation coefficients after $\ell$ levels ($a_0$ is the measured signal) and $d_\ell$ for the details; Mallat indexes by resolution $2^j$ with $j = -\ell$, while this article, like PyWavelets, counts levels. In Mallat’s normalization, decomposition is a correlation with $h$ or $g$ followed by keeping every other sample (eqs. 16 and 28),

$$a_{\ell+1}[n] = \sum_k h[k - 2n]\, a_\ell[k], \qquad d_{\ell+1}[n] = \sum_k g[k - 2n]\, a_\ell[k],$$

and reconstruction inserts a zero between samples, filters and adds (eq. 32),

$$a_\ell[n] = 2\sum_k h[n - 2k]\, a_{\ell+1}[k] + 2\sum_k g[n - 2k]\, d_{\ell+1}[k] .$$

Each level halves the number of samples in both outputs, so $J$ levels of an $N$-sample signal give $N/2 + N/4 + \dots + N/2^J$ detail samples and $N/2^J$ approximation samples: $N$ in all. Libraries such as PyWavelets scale the filters so that $\sum_n h[n] = \sqrt{2}$ (the Haar low-pass filter is $[1/\sqrt2, 1/\sqrt2]$), which for orthogonal wavelets makes each analysis step an orthonormal transform and, by our arithmetic, absorbs the factor 2 in the synthesis formula.

Two dimensions

For images Mallat uses a separable multiresolution approximation, $\Phi(x, y) = \phi(x)\phi(y)$, and his Theorem 4 gives three wavelets,

$$\Psi^1(x, y) = \phi(x)\psi(y), \qquad \Psi^2(x, y) = \psi(x)\phi(y), \qquad \Psi^3(x, y) = \psi(x)\psi(y),$$

whose dilations and translations form an orthonormal basis of the detail space at each resolution. One level filters the rows with $h$ and $g$, keeps every other row, filters the columns of both results and keeps every other column, giving an approximation and three detail images. With ideal filters, the first detail image holds the vertical high frequencies (horizontal edges), the second the horizontal high frequencies (vertical edges), and the third high frequencies in both directions, which Mallat calls corners. After $J$ levels the image is represented by $3J + 1$ images; if the image has $N$ pixels, each of the four images at level $\ell$ has $N/4^\ell$, and Mallat notes that the total equals the number of pixels. He contrasts this with the Laplacian pyramid, which increases the number of samples by a factor of 2 in one dimension and $4/3$ in two.

Only the approximation is split again at each level. PyWavelets’ documentation describes two alternatives: splitting every axis fully before moving to the next (the fully separable transform, at twice the cost, which it says represents anisotropic data more compactly), and splitting the detail bands too (the wavelet packet transform).

Biorthogonal wavelets

An orthogonal filter bank uses the same filters, reversed, for analysis and synthesis. A biorthogonal one uses different decomposition and reconstruction filters, chosen so that reconstruction is still exact although neither set of basis functions is orthogonal. PyWavelets lists its Daubechies wavelets as orthogonal and asymmetric and its biorthogonal family as non-orthogonal and symmetric. The price is that the transform no longer preserves energy, which the example below measures.

Properties

• Perfect reconstruction. For both orthogonal and biorthogonal filter banks the inverse transform returns the input up to floating-point rounding, whatever the image.
• Critical sampling, and borders. A finite image needs values beyond its edges for the filters to read. PyWavelets offers several extensions (zero, constant, symmetric, reflect, periodic, smooth, antisymmetric, antireflect) and a periodization mode; its documentation says the extension modes make the transform “slightly redundant” but ensure perfect reconstruction, while periodization gives the smallest possible number of coefficients. Its length rule is $\lfloor (n + L - 1)/2 \rfloor$ coefficients per band for a length-$n$ axis and a filter of length $L$, or $\lceil n/2 \rceil$ under periodization. The count of exactly $N$ therefore holds for periodization (and for the Haar filter on even lengths), not for every implementation. Mallat avoided border problems by assuming the signal symmetric about its first and last samples.
• Orthogonality and energy. An orthonormal wavelet basis preserves energy: the sum of squared coefficients equals the sum of squared pixels (Parseval’s relation for this basis). Mallat uses this to split a coding error among the bands: the mean squared error of the reconstruction is a weighted sum of the errors in each band. White noise in the image stays white, with the same variance, in every band of an orthonormal transform, which scikit-image’s documentation gives as the reason thresholding methods assume orthogonal wavelets.
• Vanishing moments. A wavelet has $p$ vanishing moments when $\int x^m \psi(x)\,dx = 0$ for $m = 0, \dots, p - 1$; for the filter, $\sum_n n^m g[n] = 0$ for the same $m$. It follows that every detail coefficient computed over a polynomial piece of degree below $p$ is zero: writing the samples under the filter as a polynomial in the filter index, each power is cancelled by one of the vanishing sums. A smooth region, which is locally close to a low-degree polynomial, produces small details; large details occur only where the image is not smooth. Mallat notes that the local regularity of a function can be read from how fast its wavelet coefficients decrease across scales. PyWavelets records the Haar wavelet as having one vanishing moment, Daubechies’ dbN as having $N$, and bior4.4 four.
• Localization in space and scale. Each detail coefficient depends on a neighbourhood of about one filter length at its scale. In frequency, Mallat describes the wavelet as a band-pass filter whose energy is essentially concentrated in $\pi \le |\omega| \le 2\pi$, so that each resolution’s detail covers about one octave, while noting that the channels of real filters overlap.
• Orientation. The three detail bands separate horizontal, vertical and diagonal structure, a coarse orientation selectivity. Mallat notes that non-separable wavelet bases can have as many orientation tunings as desired.
• Shift variance. Because each level subsamples, shifting the image by one pixel does not, in our reading, simply shift the coefficients but changes them, as the image pyramid article notes for pyramids. PyWavelets’ documentation calls the stationary wavelet transform, which does not decimate, a translation-invariance modification of the DWT.
• Coefficient statistics. Mallat found that the histograms of detail images of natural images are, at every resolution and orientation, symmetric peaks centred on zero, and modelled them as $h(u) = K e^{-(|u|/\alpha)^\beta}$, with $\alpha$ and $\beta$ fitted from the first two moments. In our reading, peaked, heavy-tailed histograms are another way of saying that most coefficients are small and a few are large.

Wavelet families

Family	PyWavelets name	Orthogonal	Symmetric	Vanishing moments of $\psi$	Filter length
Haar	haar (db1)	yes	no	1	2
Daubechies	dbN	yes	no	$N$	$2N$
Biorthogonal	biorNr.Nd	no	yes	first number of the name	depends on the pair

The table follows PyWavelets’ source (wavelets.c), which sets these attributes. Mallat (1989) reports the result of Daubechies (1988) that for any $n$ there is a filter whose wavelet has compact support and is $n$ times continuously differentiable. The smooth wavelets of Mallat’s own paper, built from splines and studied by Lemarié and Battle, are exponentially decreasing. In our reading, longer filters buy vanishing moments and smoothness at the price of wider support, so each coefficient sees a larger neighbourhood and an edge touches more of them.

Variants

• Stationary (undecimated) wavelet transform. Omitting the subsampling makes the transform translation-invariant and redundant: every band keeps the size of the image. PyWavelets’ documentation equates it with the algorithme à trous, whose filters are upsampled by inserting zeros; the dilated convolution article traces how that algorithm gave atrous convolution its name.
• Wavelet packets and the fully separable transform. These split more of the frequency plane, as described above.
• Continuous wavelet transform. Grossmann and Morlet, whom Mallat credits with introducing wavelets, expanded functions in translations and dilations over continuous parameters. PyWavelets implements a one-dimensional continuous transform alongside its discrete ones.
• Overcomplete oriented transforms. Steerable pyramids split each band into more orientations and are nearly free of aliasing, at a redundancy of $4k/3$ for $k$ orientations; the image pyramid article compares them with critically sampled wavelets.

Two families are called wavelets but are not the DWT described here. Gabor wavelets are complex sinusoids under Gaussian envelopes at many scales and orientations; they form a non-orthogonal family and are covered in Gabor filter. Haar-like features, the rectangle features of Viola and Jones, are named after the Haar wavelet but come in any size, position and aspect ratio and include three-rectangle bars; the Haar-like features article explains the difference. SURF’s “Haar wavelet responses” are box-filter differences of that kind.

Examples

Two examples, with NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2 and PyWavelets 1.8.0, run on scikit-image’s bundled camera image (512 × 512, grey levels 0–255, converted to floating point) and on synthetic images. Arrays are indexed [row, column], with rows running down the image. PyWavelets’ dwt2 returns cA, (cH, cV, cD), where cH is the detail along axis 0 (down the rows), cV along axis 1 and cD along both; its documentation calls axis 0 “horizontal” in the mathematical convention, which in image terms means that cH responds to horizontal edges, intensity changes from one row to the next. All transforms use four levels.

The first example checks the exact properties: perfect reconstruction and coefficient counts in two boundary modes, a hand-written filter bank against the library, energy preservation, which band an edge lands in, and vanishing moments on polynomial images. The polynomial images were sized at $64 \times 64$ when the example was planned so that the fourth-order term, which only bior4.4 sees, stands well above rounding error; that part illustrates the derivation above rather than testing it. Before running, we expected every check to come out as the derivations say, within planned tolerances: reconstruction errors below $10^{-8}$, relative interior details below $10^{-10}$ where the degree is below the number of vanishing moments and above $10^{-9}$ otherwise, and differences from the library below $10^{-10}$. The script runs in well under a second and printed identical output on two runs.

import os
os.environ["OMP_NUM_THREADS"] = "1"          # one BLAS thread, set before NumPy loads
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
import pywt
from skimage import data

WAVELETS = ["haar", "db2", "bior4.4"]
MODES = ["periodization", "symmetric"]
LEVELS = 4
img = data.camera().astype(np.float64)       # 512 x 512, grey levels 0-255
N = img.shape[0]


def count(coeffs):
    return coeffs[0].size + sum(band.size for level in coeffs[1:] for band in level)


def predicted_count(n, L, mode):
    """Coefficient count from PyWavelets' length rule (dwt_coeff_len)."""
    total = 0
    for _ in range(LEVELS):
        n = (n + 1) // 2 if mode == "periodization" else (n + L - 1) // 2
        total += 3 * n * n                   # cH, cV, cD at this level
    return total + n * n                     # plus the coarsest approximation


print("1. Reconstruction error and coefficient count, 4 levels on camera (512 x 512)")
for w in WAVELETS:
    L = pywt.Wavelet(w).dec_len
    for mode in MODES:
        c = pywt.wavedec2(img, w, mode=mode, level=LEVELS)
        err = np.abs(pywt.waverec2(c, w, mode=mode) - img).max()
        print(f"   {w:8s} {mode:13s} max error {err:.1e}   coefficients {count(c):6d} "
              f"(rule {predicted_count(N, L, mode):6d}) = {count(c) / img.size:.4f} x pixels")
lap = sum((N >> l) ** 2 for l in range(LEVELS + 1))  # Laplacian bands 512..64 plus a 32 x 32 top
print(f"   Laplacian pyramid, same depth: {lap} = {lap / img.size:.4f} x pixels")


def analyze(x, lo, hi, axis):
    """One periodized analysis step along one axis: filter, keep every second sample."""
    x = np.moveaxis(x, axis, -1)
    n, L = x.shape[-1], len(lo)
    k = np.arange(n // 2)
    a = sum(lo[j] * x[..., (2 * k + L // 2 - j) % n] for j in range(L))
    d = sum(hi[j] * x[..., (2 * k + L // 2 - j) % n] for j in range(L))
    return np.moveaxis(a, -1, axis), np.moveaxis(d, -1, axis)


print("2. Hand-written filter bank against pywt.dwt2 (one level, periodization)")
for w in WAVELETS:
    W = pywt.Wavelet(w)
    lo0, hi0 = analyze(img, W.dec_lo, W.dec_hi, axis=0)     # along axis 0 (down the rows)
    aa, ad = analyze(lo0, W.dec_lo, W.dec_hi, axis=1)       # then along axis 1
    da, dd = analyze(hi0, W.dec_lo, W.dec_hi, axis=1)
    cA, (cH, cV, cD) = pywt.dwt2(img, w, mode="periodization")
    diff = max(np.abs(p - q).max() for p, q in [(aa, cA), (da, cH), (ad, cV), (dd, cD)])
    print(f"   {w:8s} max difference {diff:.1e}")

print("3. Energy of the coefficients / energy of the zero-mean image, minus 1")
z = img - img.mean()
for w in WAVELETS:
    arr, _ = pywt.coeffs_to_array(pywt.wavedec2(z, w, mode="periodization", level=LEVELS))
    print(f"   {w:8s} {np.sum(arr ** 2) / np.sum(z ** 2) - 1:.1e}")

print("4. A horizontal edge (rows 101 and below set to 255): largest |coefficient| per band")
step = np.zeros((N, N))
step[101:, :] = 255.0
for w in WAVELETS:
    cA, (cH, cV, cD) = pywt.dwt2(step, w, mode="periodization")
    print(f"   {w:8s} cH {np.abs(cH).max():.3f}   cV {np.abs(cV).max():.1e}   cD {np.abs(cD).max():.1e}")

cp, _ = pywt.coeffs_to_array(pywt.wavedec2(img, "haar", mode="periodization", level=LEVELS))
cs, _ = pywt.coeffs_to_array(pywt.wavedec2(img, "haar", mode="symmetric", level=LEVELS))
print(f"5. haar, symmetric against periodization: max difference {np.abs(cp - cs).max():.1e}")

print("6. Vanishing moments: |sum_n n^j g[n]| for j = 0..5, and interior detail on")
print("   f = 100 (u + 0.5 v + 0.3)^p (64 x 64) as max|detail| / max|f| for p = 0..4")
P = 64
u = np.arange(P)[:, None] / P
v = np.arange(P)[None, :] / P
for w in WAVELETS:
    W = pywt.Wavelet(w)
    n = np.arange(W.dec_len, dtype=float)
    moments = [abs(np.sum(n ** j * np.array(W.dec_hi))) for j in range(6)]
    inner = slice(W.dec_len, P // 2 - W.dec_len)              # away from the periodic wrap
    rel = []
    for p in range(5):
        f = 100.0 * (u + 0.5 * v + 0.3) ** p
        _, bands = pywt.dwt2(f, w, mode="periodization")
        rel.append(max(np.abs(b[inner, inner]).max() for b in bands) / np.abs(f).max())
    print(f"   {w:8s} ({W.vanishing_moments_psi} vanishing moments) moments "
          + " ".join(f"{m:.1e}" for m in moments))
    print(f"   {'':8s} detail for p = 0..4: " + " ".join(f"{r:.1e}" for r in rel))

Output:

1. Reconstruction error and coefficient count, 4 levels on camera (512 x 512)
   haar     periodization max error 4.5e-13   coefficients 262144 (rule 262144) = 1.0000 x pixels
   haar     symmetric     max error 4.5e-13   coefficients 262144 (rule 262144) = 1.0000 x pixels
   db2      periodization max error 5.1e-13   coefficients 262144 (rule 262144) = 1.0000 x pixels
   db2      symmetric     max error 4.5e-13   coefficients 266539 (rule 266539) = 1.0168 x pixels
   bior4.4  periodization max error 1.1e-09   coefficients 262144 (rule 262144) = 1.0000 x pixels
   bior4.4  symmetric     max error 1.0e-09   coefficients 278191 (rule 278191) = 1.0612 x pixels
   Laplacian pyramid, same depth: 349184 = 1.3320 x pixels
2. Hand-written filter bank against pywt.dwt2 (one level, periodization)
   haar     max difference 1.1e-13
   db2      max difference 1.7e-13
   bior4.4  max difference 2.3e-13
3. Energy of the coefficients / energy of the zero-mean image, minus 1
   haar     1.1e-15
   db2      8.9e-16
   bior4.4  -4.7e-02
4. A horizontal edge (rows 101 and below set to 255): largest |coefficient| per band
   haar     cH 255.000   cV 1.6e-15   cD 7.8e-16
   db2      cH 174.168   cV 2.7e-14   cD 9.7e-15
   bior4.4  cH 142.174   cV 5.4e-10   cD 1.4e-10
5. haar, symmetric against periodization: max difference 0.0e+00
6. Vanishing moments: |sum_n n^j g[n]| for j = 0..5, and interior detail on
   f = 100 (u + 0.5 v + 0.3)^p (64 x 64) as max|detail| / max|f| for p = 0..4
   haar     (1 vanishing moments) moments 0.0e+00 7.1e-01 7.1e-01 7.1e-01 7.1e-01 7.1e-01
            detail for p = 0..4: 3.7e-17 8.8e-03 1.7e-02 2.3e-02 2.9e-02
   db2      (2 vanishing moments) moments 2.8e-17 5.6e-17 1.2e+00 4.5e+00 1.3e+01 3.8e+01
            detail for p = 0..4: 1.3e-16 1.4e-16 1.3e-04 3.6e-04 6.3e-04
   bior4.4  (4 vanishing moments) moments 1.4e-12 5.7e-12 2.3e-11 9.1e-11 8.3e+00 1.7e+02
            detail for p = 0..4: 2.0e-12 1.5e-12 1.1e-12 7.8e-13 7.0e-08

Every reconstruction is exact to rounding, within the planned tolerance of $10^{-8}$: about $5 \times 10^{-13}$ grey levels for the orthogonal wavelets and $10^{-9}$ for bior4.4. In our reading the larger error comes from bior4.4‘s stored filters, whose printed moments, exactly zero below order 4 in theory, are of order $10^{-12}$ to $10^{-11}$. The counts agree with the length rule in every case. Periodization keeps exactly one coefficient per pixel for all three wavelets; the symmetric mode, PyWavelets’ default, adds 1.7% for db2 and 6.1% for bior4.4, whose filters (length 4 and 10 in PyWavelets’ source) reach further past the border, and nothing for Haar, whose two taps never leave the image on even lengths; Haar’s coefficients are in fact identical in both modes. A Laplacian pyramid of the same depth holds 1.33 times as many samples as the image.

The NumPy filter bank, a dozen lines that filter and keep every second sample, reproduces pywt.dwt2 to $10^{-13}$, which pins down the alignment PyWavelets uses under periodization (coefficient $k$ reads samples $2k + L/2 - j$, modulo the length). The energy check separates the two kinds of wavelet: Haar and db2 preserve energy to $10^{-15}$, while bior4.4’s coefficients hold 4.7% less energy than the image, so for a biorthogonal transform, coefficient energy is not image energy. The horizontal edge appears only in cH; cV and cD are zero to rounding, as they must be for an image that does not change along its rows.

The vanishing-moment rows show the derivation at work. The filter moments vanish up to order 0 for Haar, 1 for db2 and 3 for bior4.4, and the interior details of the polynomial images vanish (to between $10^{-17}$ and $2 \times 10^{-12}$) exactly for the degrees below each wavelet’s count: a plane already produces Haar details of almost 1% of the image’s maximum, but db2 details only from the quadratic on, and bior4.4 details only for the quartic, at $7 \times 10^{-8}$.

The second example measures sparsity on the zero-mean camera image: for each representation, the smallest number $k$ of largest-magnitude coefficients whose reconstruction has a relative squared error of at most 1% or 0.1%. The pixel basis and a global orthonormal 2-D discrete cosine transform (type II) serve as baselines. For orthonormal transforms the error is simply the energy of the dropped coefficients; for bior4.4 and for the symmetric mode, which are not orthogonal, the code searches over $k$ on the actual reconstruction, checks that $k$ meets the target and $k - 1$ does not, and for Haar and db2 under periodization compares the search with the energy count. The boundary mode is a modelling choice that the sources leave open, so both modes are printed.

The planned hypotheses were that every wavelet needs fewer coefficients than the pixel basis (H1); that under periodization db2 and bior4.4 each need fewer than Haar (H2); that bior4.4 needs fewer than the global DCT in each mode (H3); and that db2 and bior4.4 need fewer under the symmetric mode than under periodization, which wraps each border onto the opposite one and so creates an artificial edge (H4). All at both targets.

The last part is one brief thresholding run; image denoising covers wavelet shrinkage and its successors. Gaussian noise with $\sigma = 20$ grey levels is added (seeds 0–4, no clipping), the noise level is estimated from the finest diagonal band by the median absolute value of its nonzero coefficients divided by $\Phi^{-1}(0.75) \approx 0.6745$, and every detail coefficient is thresholded at the universal threshold $\hat\sigma\sqrt{2\ln n}$, with $n$ the number of pixels, about $5\hat\sigma$ here, with soft thresholding (shrink every kept coefficient toward zero by the threshold) and hard thresholding (keep it unchanged). This is VisuShrink as scikit-image implements it; scikit-image’s source credits the threshold and the noise estimate (its section 4.2) to Donoho and Johnstone (1994), and the code checks the hand version against skimage.restoration.denoise_wavelet. Quality is PSNR. We expected thresholding to raise the PSNR above the noisy image’s on every seed (H5a) and hard thresholding to beat soft on every seed (H5b), since soft thresholding at about 100 grey levels shifts every kept coefficient by that much. The script runs in about 2 seconds and printed identical output on two runs.

import os
os.environ["OMP_NUM_THREADS"] = "1"          # one BLAS thread, set before NumPy loads
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
import warnings
import numpy as np
import pywt
from scipy import fft, stats
from skimage import data
from skimage.restoration import denoise_wavelet

WAVELETS = ["haar", "db2", "bior4.4"]
MODES = ["periodization", "symmetric"]
LEVELS = 4
TARGETS = [1e-2, 1e-3]                       # allowed relative squared error
img = data.camera().astype(np.float64)
z = img - img.mean()
E = np.sum(z ** 2)


def k_orthonormal(coeffs, target):
    """Orthonormal transform: the error of keeping the k largest is the energy dropped."""
    s = np.sort(coeffs.ravel() ** 2)[::-1]
    return int(np.searchsorted(np.cumsum(s) >= (1 - target) * E, True) + 1)


def k_wavelet(w, mode, target):
    """Smallest k whose reconstruction from the k largest |coefficients| meets the target."""
    arr, slices = pywt.coeffs_to_array(pywt.wavedec2(z, w, mode=mode, level=LEVELS))
    order = np.argsort(np.abs(arr), axis=None)[::-1]

    def err(k):
        kept = np.zeros(arr.size)
        kept[order[:k]] = arr.ravel()[order[:k]]
        c = pywt.array_to_coeffs(kept.reshape(arr.shape), slices, output_format="wavedec2")
        r = pywt.waverec2(c, w, mode=mode)[: z.shape[0], : z.shape[1]]
        return np.sum((r - z) ** 2) / E

    lo, hi = 1, arr.size
    while lo < hi:                                 # binary search on the real error
        mid = (lo + hi) // 2
        lo, hi = (lo, mid) if err(mid) <= target else (mid + 1, hi)
    local = err(lo) <= target < err(lo - 1)
    exact = None
    if mode == "periodization" and pywt.Wavelet(w).orthogonal:
        exact = lo == k_orthonormal(arr, target)
    return lo, local, exact


print("A. Coefficients needed (k, share of 262,144) for 1% and 0.1% relative squared error")
for t in TARGETS:
    print(f"   target {t:.0e}")
    for name, coeffs in [("pixels", z), ("global DCT", fft.dctn(z, norm="ortho"))]:
        k = k_orthonormal(coeffs, t)
        print(f"   {name:24s} k = {k:6d} ({100 * k / z.size:7.3f}%)")
    for mode in MODES:
        for w in WAVELETS:
            k, local, exact = k_wavelet(w, mode, t)
            extra = "" if exact is None else f", equals cumulative-energy count: {exact}"
            print(f"   {w + ' ' + mode:24s} k = {k:6d} ({100 * k / z.size:7.3f}%)"
                  f"  local check: {local}{extra}")

SIGMA = 20.0


def visushrink(noisy, w, how):
    c = pywt.wavedec2(noisy, w, mode="symmetric", level=LEVELS)
    finest = c[-1][2]                                    # cD at level 1
    finest = finest[finest != 0]
    sigma = np.median(np.abs(finest)) / stats.norm.ppf(0.75)   # MAD estimate
    t = sigma * np.sqrt(2 * np.log(noisy.size))          # universal threshold
    c = [c[0]] + [tuple(pywt.threshold(b, t, mode=how) for b in lev) for lev in c[1:]]
    return pywt.waverec2(c, w, mode="symmetric")[: noisy.shape[0], : noisy.shape[1]], sigma


def psnr(x):
    return 10 * np.log10(255.0 ** 2 / np.mean((x - img) ** 2))


print(f"B. VisuShrink on camera + Gaussian noise (sigma = {SIGMA:.0f}), seeds 0-4, PSNR in dB")
noisy_psnr, sigmas, results, worst = [], [], {}, 0.0
for seed in range(5):
    rng = np.random.default_rng(seed)
    noisy = img + SIGMA * rng.standard_normal(img.shape)
    noisy_psnr.append(psnr(noisy))
    for w in WAVELETS:
        for how in ["soft", "hard"]:
            out, sigma = visushrink(noisy, w, how)
            results.setdefault((w, how), []).append(psnr(out))
            with warnings.catch_warnings():           # bior4.4 is not orthogonal: skimage warns
                warnings.simplefilter("ignore")
                lib = denoise_wavelet(noisy, wavelet=w, mode=how, wavelet_levels=LEVELS,
                                      method="VisuShrink")
            worst = max(worst, np.abs(out - lib).max())
            if how == "soft":                         # one estimate per wavelet and seed
                sigmas.append(sigma)
print(f"   noisy            {np.mean(noisy_psnr):.2f} ({min(noisy_psnr):.2f}-{max(noisy_psnr):.2f})")
for (w, how), vals in results.items():
    print(f"   {w:8s} {how:4s}    {np.mean(vals):.2f} ({min(vals):.2f}-{max(vals):.2f})")
print(f"   MAD estimate of sigma (3 wavelets x 5 seeds): {min(sigmas):.3f}-{max(sigmas):.3f}")
print(f"   max difference from skimage denoise_wavelet: {worst:.1e}")

Output:

A. Coefficients needed (k, share of 262,144) for 1% and 0.1% relative squared error
   target 1e-02
   pixels                   k = 208084 ( 79.378%)
   global DCT               k =  28895 ( 11.023%)
   haar periodization       k =  12759 (  4.867%)  local check: True, equals cumulative-energy count: True
   db2 periodization        k =  12630 (  4.818%)  local check: True, equals cumulative-energy count: True
   bior4.4 periodization    k =  10985 (  4.190%)  local check: True
   haar symmetric           k =  12759 (  4.867%)  local check: True
   db2 symmetric            k =  11732 (  4.475%)  local check: True
   bior4.4 symmetric        k =  11438 (  4.363%)  local check: True
   target 1e-03
   pixels                   k = 239371 ( 91.313%)
   global DCT               k = 130921 ( 49.942%)
   haar periodization       k =  57601 ( 21.973%)  local check: True, equals cumulative-energy count: True
   db2 periodization        k =  56800 ( 21.667%)  local check: True, equals cumulative-energy count: True
   bior4.4 periodization    k =  54443 ( 20.768%)  local check: True
   haar symmetric           k =  57601 ( 21.973%)  local check: True
   db2 symmetric            k =  55931 ( 21.336%)  local check: True
   bior4.4 symmetric        k =  56614 ( 21.597%)  local check: True
B. VisuShrink on camera + Gaussian noise (sigma = 20), seeds 0-4, PSNR in dB
   noisy            22.11 (22.10-22.12)
   haar     soft    24.33 (24.31-24.36)
   haar     hard    25.91 (25.86-25.94)
   db2      soft    24.77 (24.75-24.78)
   db2      hard    26.27 (26.24-26.31)
   bior4.4  soft    24.88 (24.88-24.90)
   bior4.4  hard    26.57 (26.56-26.58)
   MAD estimate of sigma (3 wavelets x 5 seeds): 20.255-21.134
   max difference from skimage denoise_wavelet: 0.0e+00

The pixel basis needs 79% of its values to reach 1% error, the global DCT 11%, and the wavelets between 4.2% and 4.9%; at 0.1% error the pixel basis needs 91%, the DCT 50% and the wavelets 20.8–22.0%. H1 held under both modes at both targets, and so did H3: bior4.4 needs fewer than half the DCT’s coefficients in either mode. In our reading, the difference is the edges: a DCT coefficient spans the whole image, so each edge costs coefficients at many frequencies, while a wavelet coefficient spans a neighbourhood, so an edge costs a few coefficients per scale and the flat regions cost almost none.

H2 held: under periodization, db2 needs fewer coefficients than Haar and bior4.4 fewer still, at both targets. The margin of db2 over Haar is small, 129 coefficients at 1% and 801 at 0.1%, 1.0% and 1.4% of Haar’s count; bior4.4’s margin is larger (1,774 and 3,158). H4 held for db2, which needs fewer coefficients under the symmetric mode than under periodization at both targets, and failed for bior4.4, which needs more under the symmetric mode at both targets (11,438 against 10,985, and 56,614 against 54,443). The example does not show why. Haar, as the first example showed, is unaffected by the mode. The searches satisfy their local check everywhere, and for Haar and db2 they agree with the energy count.

Thresholding raised the PSNR by about 2.2–4.5 dB over the noisy 22.1 dB for every wavelet and both rules; the lowest denoised value on any seed (24.31 dB) is above the highest noisy one (22.12 dB), so H5a held. Hard thresholding beat soft by 1.5–1.7 dB in the means, and for every wavelet the lowest hard value is above the highest soft one, so H5b held too. The noise estimate came out between 20.255 and 21.134 against a true 20, inside the planned 10%; in our reading, the image’s own fine detail in the diagonal band adds to the noise there. The hand implementation matches scikit-image exactly, so every planned check of both examples passed. scikit-image warns that bior4.4 is not orthogonal and that its thresholds may then be chosen poorly; its documentation also says VisuShrink tends to produce images that look overly smooth, and offers BayesShrink, a threshold per band, as the default.

Common Misconceptions

• “A wavelet transform always has as many coefficients as pixels.” Only the critically sampled transform with periodization (or Haar on even lengths) does. Common boundary modes add coefficients, 1.7% and 6.1% in the example, and the stationary transform keeps a full-size image per band.
• “The wavelet coefficients carry the image’s energy.” Only for orthonormal wavelets. A biorthogonal transform reconstructs exactly but does not preserve energy (4.7% less in the example), so it follows that the energy of the dropped coefficients is no longer the error of a truncation, and the example has to measure that error on the reconstruction.
• “Haar-like features and Gabor wavelets are wavelet transforms.” Neither is the DWT. Haar-like features are an overcomplete set of rectangle differences evaluated with an integral image, and Gabor wavelets an overcomplete family of oriented complex filters.
• “cH holds the horizontal frequencies.” In PyWavelets, cH holds the detail along axis 0, which in an image is the vertical direction: it responds to horizontal edges. Mallat’s first detail image likewise holds the vertical high frequencies, which he describes as horizontal edges. Naming conventions differ, so check which axis a library filters.
• “Wavelets are better than the DCT.” For whole-image approximation of camera they need far fewer coefficients, as printed. Block DCTs, as in JPEG, are a different representation; this example does not compare them.

Where It Is Used

• Denoising. Thresholding wavelet coefficients removes noise because the signal is concentrated in a few large coefficients and white noise is spread over all of them; the image denoising article covers Donoho and Johnstone’s shrinkage, its Gibbs-like artifacts near edges, and the translation-invariant variants that reduce them. Wavelet detail coefficients are also used to estimate the noise level itself, as the image noise article notes.
• Compression. Mallat (1989) proposed adapting the coding error of each band to the sensitivity of the visual system for that frequency and orientation, and reported that, with such representations, images can be coded at less than 1.5 bits per pixel with few visible distortions, citing his own work and two papers on pyramid and subband image coding. The image processing article names compression among the main uses of wavelets.
• Reconstruction from incomplete data. Compressed-sensing MRI needs a transform in which the image is sparse, such as wavelets or finite differences; the medical image reconstruction article describes it.
• Texture and image statistics. Mallat interpreted wavelet coefficients as Julesz’s textons and used their statistics for texture discrimination and fractal analysis. The texture analysis and image quality assessment articles describe later methods built on wavelet-coefficient statistics, such as Portilla and Simoncelli’s texture model and the VIF and RRED quality indices.
• Networks. The algorithme à trous for the undecimated transform lives on as atrous or dilated convolution.

History

Mallat’s introduction places the transform at the meeting of several lines of work. Grossmann and Morlet introduced wavelets as functions whose translations and dilations could expand square-integrable functions; Meyer showed that there are wavelets whose dyadic dilations and translations form orthonormal bases of $L^2(\mathbb{R})$, generalizing the Haar basis, where before, Mallat writes, it had been believed that no construction could give simple orthonormal bases with good localization in both space and frequency; and Burt, with the Laplacian pyramid, and Crowley had introduced pyramidal implementations for computing image details at several resolutions. Mallat’s contribution, presented at an IEEE computer vision workshop in December 1987 (his reference 27) and in the 1989 journal paper, was to derive wavelet bases from multiresolution approximations and to compute them with a pyramid of quadrature mirror filters, filters his references trace to split-band voice coding and subband image coding. Daubechies (1988), as Mallat reports, constructed orthonormal wavelets with compact support, which have finite filters. Donoho and Johnstone (1994) proposed thresholding wavelet coefficients for denoising, as the image denoising article describes.