Wavelets & Time-Frequency Trade-offs

"The Fourier transform knows every note in the symphony but refuses to say when any of them was played. I know the what and the where, within limits my lawyer describes as 'Heisenberg-compliant'."

A Well-Localized Wavelet

The pyramids of the previous section decomposed images by scale, but they are slightly wasteful (a Laplacian pyramid carries about a third more numbers than the image it encodes) and their levels are not independent. Wavelets fix both at once, and they answer a complaint you may have been nursing since Section 4.1: a spectrum tells you a strong diagonal texture exists somewhere in the image, but not where. For inspection, compression, and denoising, where matters. This section explains why no transform can answer both questions perfectly, then builds the family that comes closest.

1. What Fourier Forgets, and the Price of Remembering Intermediate

Each Fourier basis wave spans the entire image: a single coefficient $F(u,v)$ mixes contributions from every pixel. That is why one bright dead pixel perturbs the whole spectrum, and why the spectrum of "stripes on the left half" and "stripes on the right half" can look nearly identical. The pixel basis sits at the opposite extreme: perfect location, zero frequency knowledge. The obvious compromise is to window the analysis, chopping the image into patches and transforming each, the short-time Fourier transform (STFT), whose 2D incarnation with Gaussian windows gives the Gabor filters that reappear as texture features later in this book and, remarkably, as good models of the receptive fields in your own primary visual cortex.

But windowing has a built-in flaw: one window size must serve all frequencies. A window wide enough to estimate low frequencies precisely smears the location of fine details; a window narrow enough to localize an edge cannot tell a 30-cycle wave from a 32-cycle one. The constraint is fundamental, not an engineering failure. For any analysis function with spatial spread $\sigma_x$ and frequency spread $\sigma_u$,

$$\sigma_x \, \sigma_u \;\ge\; \frac{1}{4\pi}$$

the Gabor uncertainty principle, mathematically the same inequality as Heisenberg's. You cannot know both precisely; you can only choose how to spend your uncertainty budget. The wavelet insight is that the right spending plan is frequency-dependent: analyze low frequencies with wide windows (good frequency resolution, who cares exactly where the sky gradient sits) and high frequencies with narrow windows (good localization, which is the whole point of an edge). Figure 4.6.1 compares the four budgets.

2. The Haar Wavelet From Scratch Intermediate

The simplest wavelet transform, due to Alfréd Haar in 1910, predates the word "wavelet" by seven decades and needs only addition and subtraction. Take a row of pixels two at a time and record, for each pair, its scaled average $a = (p_1 + p_2)/\sqrt{2}$ and its scaled difference $d = (p_1 - p_2)/\sqrt{2}$. Averages summarize; differences remember what summarizing lost. (The $\sqrt{2}$ keeps total energy unchanged, making the transform orthogonal.) In 2D, apply the pairing first along rows, then along columns of both results, yielding four quarter-size subbands:

import numpy as np
from skimage import data

img = data.camera().astype(np.float64)            # 512 x 512

def haar2d_level(a):
    """One level of the 2D Haar transform: returns LL, (LH, HL, HH)."""
    lo = (a[:, 0::2] + a[:, 1::2]) / np.sqrt(2)   # pair averages along x
    hi = (a[:, 0::2] - a[:, 1::2]) / np.sqrt(2)   # pair differences along x
    LL = (lo[0::2, :] + lo[1::2, :]) / np.sqrt(2) # then the same along y
    LH = (lo[0::2, :] - lo[1::2, :]) / np.sqrt(2)
    HL = (hi[0::2, :] + hi[1::2, :]) / np.sqrt(2)
    HH = (hi[0::2, :] - hi[1::2, :]) / np.sqrt(2)
    return LL, (LH, HL, HH)

LL, (LH, HL, HH) = haar2d_level(img)
print(LL.shape)    # (256, 256): four quarter-size bands, zero redundancy
# LL is a clean half-resolution image; LH lights up on horizontal edges,
# HL on vertical edges, HH on diagonals and fine texture.

Display the four bands and you will recognize old friends: LL is a properly anti-aliased half-size image (it is a Gaussian-pyramid level, by humbler means), while LH, HL, and HH split the lost detail by orientation, the same band that a Laplacian level stores all mixed together. The multi-level transform now writes itself: the LL band is just a smaller image, so transform it again, and again, until the approximation is a handful of pixels. The result, the discrete wavelet transform (DWT), is the critically sampled cousin of Section 4.5's pyramid, with the nested layout shown in Figure 4.6.2.

3. Multi-Level DWT with PyWavelets Beginner

Haar is the clearest wavelet but not the best for photographs: its blocky basis leaves staircase artifacts under heavy compression. Practical systems use smoother families, Daubechies (db2, db4), symlets, biorthogonal pairs (JPEG 2000 uses the bior-family CDF 9/7 and 5/3 wavelets), all available through PyWavelets:

import pywt

coeffs = pywt.wavedec2(img, wavelet="db2", level=3)
print(coeffs[0].shape)                        # (66, 66) coarse approximation
print([lvl[0].shape for lvl in coeffs[1:]])   # [(66, 66), (130, 130), (257, 257)]

arr, slices = pywt.coeffs_to_array(coeffs)    # pack into one Figure-4.6.2 mosaic

# Energy compaction: what fraction of total energy do the largest 1% hold?
flat = np.sort(np.abs(arr).ravel())[::-1]
k = int(0.01 * flat.size)
print((flat[:k] ** 2).sum() / (flat ** 2).sum())
# Prints a number close to 1.0 (typically 0.95 to 0.99 for natural photos):
# a few coefficients carry almost everything.

That last printed number is the economic engine of everything wavelets are used for. An orthogonal transform cannot create or destroy energy, only reorganize it, and the DWT reorganizes a photograph's energy into very few large coefficients (edges, structures) and a vast desert of near-zeros (smooth regions, at every scale). Compression keeps the few and cheaply encodes the desert; denoising exploits the same sparsity from the opposite direction.

4. Wavelet Shrinkage: Denoising by Thresholding Intermediate

White noise is the one signal that refuses to compact: orthogonal transforms spread it evenly, so every wavelet coefficient receives the same small noise contribution. Real image content, by contrast, concentrates into few large coefficients. The two populations barely overlap, and that gap is an algorithm: transform, shrink every coefficient toward zero by a threshold calibrated to the noise level, reconstruct. Donoho and Johnstone's classic recipe estimates the noise from the finest diagonal band (which is almost pure noise in most photos) and applies the "universal" threshold:

rng = np.random.default_rng(0)
noisy = img + 20 * rng.standard_normal(img.shape)     # sigma = 20 noise

c = pywt.wavedec2(noisy, "db2", level=4)
sigma = np.median(np.abs(c[-1][-1])) / 0.6745         # noise std from finest HH
thresh = sigma * np.sqrt(2 * np.log(noisy.size))      # universal threshold

c_den = [c[0]] + [
    tuple(pywt.threshold(band, thresh, mode="soft") for band in level)
    for level in c[1:]
]
restored = pywt.waverec2(c_den, "db2")[: img.shape[0], : img.shape[1]]
print(round(sigma, 1))      # close to 20.0: the estimate found the noise

Compare this with the Gaussian smoothing of Chapter 3: blurring suppresses noise and edges alike, because in pixel space they occupy the same neighborhoods. In wavelet space they occupy different amplitudes, so the threshold can kill one and spare the other. This idea, denoise where the signal is sparse, is the seed of the entire restoration toolbox of Chapter 7, and its philosophical descendants reach all the way to the learned denoisers at the heart of modern generative models.

from skimage.restoration import denoise_wavelet
clean = denoise_wavelet(noisy / 255.0, method="BayesShrink",
                        mode="soft", wavelet="db2", rescale_sigma=True)

5. Wavelets at Work: JPEG 2000 and the Compression Lineage Intermediate

The JPEG format you met in Chapter 1 transforms fixed 8x8 blocks with the discrete cosine transform, a Fourier relative, and at strong compression its block boundaries surface as the familiar mosaic of artifacts. JPEG 2000 replaced the blocks with a whole-image DWT in exactly Figure 4.6.2's layout, and the change buys three structural advantages: no block boundaries (artifacts degrade into soft blur instead of mosaics), an embedded progressive bitstream (truncate the file anywhere and you still have the best possible image at that byte budget), and free multi-resolution access (the LL bands are ready-made zoom levels, so a viewer can fetch a thumbnail or a region at reduced scale without decoding the full image). Those properties made it the standard in digital cinema distribution, medical imaging archives, and map servers, even though baseline JPEG kept the consumer web on simplicity and decoder ubiquity.