Wavelets

Here a signal with three features at three scales, a slow wave, a step and a burst of fast alternation, is split into averages and differences, and the averages are split again and again. Watch where each feature lands: the burst all in the first level, the slow wave in the few averages left at the end.

Averages and differences

In Filter banks (23.1) a two-channel bank split a signal into a low band and a high band. It kept every second sample of each, and it still rebuilt the input exactly. On this page I run that bank again on its own low band, and again on the low band of that. The result is the discrete wavelet transform, the DWT.

Let’s start with the simplest pair of filters, 23.1’s Haar pair, $(1,1)/\sqrt2$ and $(1,-1)/\sqrt2$. Filter with each, and keep every second output, as in “Keep every M-th sample” of Downsampling and decimation (22.1). Each pair of samples becomes two numbers:

I call $a_1[k]$ an approximation coefficient: the average of the pair, times $\sqrt2$. And $d_1[k]$ is a detail coefficient: half the difference of the pair, times $\sqrt2$. The subscript 1 is the level. From $N$ samples you get $N/2$ of each, so there are as many numbers as before.

Depending on which outputs you keep, 23.1’s bank gives these numbers, or $-d_1$, or pairs the samples one place later. I follow PyWavelets, whose wavedec gives exactly these two formulas.

Why divide by $\sqrt2$ and not by 2? Because then the energy of How big is a signal (1.3) is shared out exactly. For one pair of samples $p$ and $q$:

The two rows $(1,1)/\sqrt2$ and $(1,-1)/\sqrt2$ are orthogonal and have length 1. They are orthonormal, and “Rows that cancel” in The DFT as a matrix (13.5) showed that such rows keep energy.

Now do the same to $a_1$. Its pairs give $a_2[k]$ and $d_2[k]$, half as many again. Then $a_2$ gives $a_3$ and $d_3$, and so on. Each level keeps its details and splits only the approximation. For three levels the energy adds up like this:

With more levels it is the same: the last approximation plus every level’s details.

What does each level hold? Up to the factor $\sqrt2$, $a_1$ is a 2-point moving average, a gentle low-pass from Simple smoothing filters (18.2). The pair difference is a gentle high-pass. So $d_1$ holds roughly the top half of the band, $\pi/2$ to $\pi$, and $a_1$ the bottom half.

Keeping every second sample stretches $a_1$‘s half band back out to fill 0 to $\pi$ (22.1). So the next split takes half of that half: $d_2$ holds roughly $\pi/4$ to $\pi/2$ of the original band. Then $d_3$ holds $\pi/8$ to $\pi/4$, and $a_3$ keeps what is below $\pi/8$. Haar’s split is gentle, so these edges are soft, as 23.1 found.

Think of zooming out of a city map. The first zoom loses the alleys, the next the streets, and the motorways are left. Each level of the DWT is one zoom: the details are what that zoom lost, and the approximation is the map that is left.

Split the low band again and again

The picture at the top of the page runs three Haar levels on a signal with three features, each at its own scale. There is a slow wave, one period in 64 samples, and a step at $n=41$, the $u[n]$ of Impulse, step and ramp (3.1). There is also a burst of fast alternation, $0.5(-1)^n$, for $n$ from 16 to 23 and 0 at every other $n$:

Its 64 samples run from $-0.707$ to $1.500$, and its energy is $\sum x[n]^2=23.087$.

Notice the burst: all of its energy is in level 1. It starts at an even sample, so each pair holds $+0.5$ and $-0.5$ on top of the wave. The average loses it, and the difference gets $1/\sqrt2=0.7071$, plus a little of the wave. The burst’s energy, $8\cdot0.25=2$, is exactly the $4\cdot0.7071^2=2$ that it adds to $d_1$.

So $a_1$ has no burst left in it already: the burst’s pairs cancelled at the first split, which is why the level-2 caption has no burst left to split. The step lands in $d_1$ because one pair, samples 40 and 41, straddles it. That pair’s difference is $-1/\sqrt2=-0.7071$, plus $+0.0466$ from the wave, which makes $-0.6605$.

Everywhere else in $d_1$ the wave changes very little from one sample to the next, so its details are at most 0.0693. At coarser levels the halves being compared lie further apart, and the wave has changed more between them. That is why $d_3$ holds more energy than $d_2$: 5.04 % against 1.91 %.

The largest $d_3$ values, $-0.5164$ and $0.5164$, sit at $k=0$ and $k=4$. Those blocks cover samples 0 to 7 and 32 to 39, where the wave is steepest. The step still shows at every level, because some pair always straddles it: in $d_2$ it gives $-0.5$, partly offset by the wave to $-0.3834$.

At level 3 the approximation is only 8 numbers, yet it keeps 81.62 % of the energy. Each $a_3[k]$ is $\sqrt8$ times the average of 8 samples. Those 8 averages trace the slow wave and the step’s height, and nothing finer.

Two ways to tile time and frequency

Each coefficient covers some samples and some band. A level-1 coefficient covers 2 samples and the top half of the band. A level-2 coefficient covers 4 samples and the next quarter. Each level down, a coefficient lasts twice as long and its band is half as wide.

Draw each coefficient as a tile, with time across and frequency up. In “Short window or long window” of Spectrograms & the STFT (15.5), one window length set the tile for every frequency at once. The figure puts that grid next to the wavelet tiles, for 16 samples and three levels.

Fast events, like the burst and the step, are made of high frequencies, and a short tile pins them down in time. A slow wave needs a long look before its frequency is clear, and it gets one. Each scale gets a tile to suit it, and this is what multiresolution means.

It suits signals whose fast parts are brief and whose slow parts last, which is common. It does not suit two close tones high in the band, the case 15.5 needed a long window for: there the wavelet tiles are tall.

Smoother wavelets: Daubechies

Nothing ties the tree to the Haar pair. 23.1’s Daubechies pair db2 has four taps, $h_0[n]$ = 0.48296, 0.83652, 0.22414 and −0.12941, and the same tree runs with it. Each filter now reaches past the end of the block, so I take the samples past the end from its start again: the block is treated as periodic.

Haar’s details are 0 wherever the signal is constant. db2’s details are 0 on straight lines too, except where the periodic wrap joins the line’s end to its start. A smooth signal looks almost straight over four samples, so it gives even fewer large details.

Take the slow wave alone. Over three levels, Haar puts 4.96 % of its energy in details, and db2 only 0.45 %, about a tenth.

Where does the name come from? Set one detail coefficient to 1 and every other coefficient to 0, and run the synthesis steps back to samples. Out comes one shape, and the deeper its level, the more samples it covers. Haar’s is a block up and then a block down; db2’s is jagged but has no jumps.

As the level grows, that shape settles to one form that is only stretched: the wavelet, a “small wave”. Doing the same with an approximation coefficient gives the scaling function. Every detail coefficient measures how much of one shifted, stretched copy of the wavelet the signal holds.

Keep the big coefficients, drop the rest

Now add noise of RMS $x_\text{noise,rms}$ to a signal, with RMS from “Four everyday sizes” in 1.3. The transform’s rows are orthonormal, so noise that is the same all along the signal spreads evenly over the coefficients. Each of the $N$ coefficients gets noise of about that same RMS.

A signal made of flat blocks is different. Its details are 0 inside every block, so it puts its energy into a few large coefficients near the edges. Most details then hold noise alone, and small.

So keep every detail larger than a threshold $\eta$, set the rest to 0, and rebuild. The approximation coefficients are always kept. This is hard thresholding. To rebuild, run each level backwards:

and the same from $a_2$ and $d_2$ back to $a_1$.

How large should $\eta$ be? Donoho and Johnstone proposed the universal threshold:

They chose it so that, as $N$ grows, the chance that any pure-noise coefficient passes it goes to 0. For $N=256$ and RMS 0.15, $\sqrt{2\ln256}=3.3302$, so $\eta=0.4995$.

The instrument’s signal is a row of blocks: 0, then 1 from $n=50$, $-0.5$ from 90, 0.7 from 150, and 0 again from 205. The noise comes from the site’s seeded generator, drawn at RMS 0.15; this draw of 256 samples happens to measure 0.1366. Its edges are not at multiples of 2, 4 or 8, so the clean signal has 17 nonzero details, not 5.

Watch the error as $\eta$ rises: 0.1366, then 0.0886, then 0.0573, less than half of where it started. Of the 16 details kept at 0.50, 15 sit where the clean signal has a detail. The other one is pure noise, 0.5076 at level 1, just above the threshold.

For comparison, a centred 8-point moving average, not drawn in the picture, leaves 0.1235 and blurs every edge. The moving average of Simple smoothing filters (18.2) has no such choice. It turns every edge into an 8-sample ramp, as in “A longer average: less noise, more delay”, and near the edges its error reaches 0.82. Thresholding keeps the large edge coefficients whole, so the rebuilt edges stay where they were.

After the clip, drag the threshold line to set $\eta$ yourself. Now try 1.2. Only 8 details survive, and the error, 0.1392, is worse than the noisy signal’s 0.1366. Too high a threshold removes the signal too: the smaller edges’ coefficients go with the noise.

Which threshold?

The universal threshold is a principled guess, not the best value. On this trace the error is smallest, 0.0477, for $\eta$ from 0.51 to 0.52, just above it. There the one noise-only coefficient, 0.5076, has gone, and all 15 signal coefficients stay. At 0.53 a signal coefficient, 0.5279, goes too.

On real data the noise RMS is not known. Donoho and Johnstone estimate it from the level-1 details, which hold mostly noise.

Soft thresholding sets the small details to 0 as well, and also shrinks every kept detail towards 0 by $\eta$. Its output has no jumps as $\eta$ changes, and it tends to look smoother. But it shrinks the large edge coefficients too: on this trace, at the universal threshold, its error is 0.1212 against hard thresholding’s 0.0573.

Worked example

1. Haar by hand. Take $x$ = 1, 3, 2, 2. The pairs give $a_1$ = $4/\sqrt2$, $4/\sqrt2$ = 2.8284, 2.8284 and $d_1$ = $-2/\sqrt2$, 0 = $-1.4142$, 0. One more level gives $a_2=(2.8284+2.8284)/\sqrt2=4$ and $d_2=0$.

The energy checks: $1+9+4+4=18$, and $4^2+0^2+(-1.4142)^2+0^2=16+2=18$. Going back, $x[0]=(2.8284-1.4142)/\sqrt2=1$ and $x[1]=(2.8284+1.4142)/\sqrt2=3$.

2. The tree. The 64-sample signal has energy 23.087. After three Haar levels, $d_1$ holds 11.44 %, $d_2$ 1.91 %, $d_3$ 5.04 % and $a_3$ 81.62 %. These add to 100 % before rounding; rounded, they make 100.01.

3. Denoising. With $N=256$ and RMS 0.15, $\ln256=5.5452$, so $\eta=0.15\sqrt{11.0904}=0.4995$. Hard thresholding keeps 16 of the 248 details, and the error falls from 0.1366 to 0.0573.

Where you’ll meet this

JPEG 2000 compresses images with a wavelet transform whose analysis and synthesis filters differ, a biorthogonal pair, which lets both be symmetric. Its two-dimensional version splits rows and columns level by level; Image compression (28.4) comes back to it. The FBI’s fingerprint archive uses a wavelet standard too, WSQ.

Wavelet denoising cleans ECGs and seismic traces, where sharp events sit in slow backgrounds; ECG processing (31.1) uses it. In machine learning, the energies of the levels make multiscale features for classifying signals.

PyWavelets reproduces this page’s wavelet numbers with pywt.wavedec(x, 'haar', mode='periodization'), pywt.threshold and pywt.waverec. For the ideas in depth, see S. Mallat, A Wavelet Tour of Signal Processing (3rd ed., 2009), and G. Strang and T. Nguyen, Wavelets and Filter Banks (1996).

Reference card