A wave known only at its ticks

Every wave that reaches a computer arrives as a list of numbers: its value at one instant, then at the next, at equal ticks of a clock. Between the ticks nothing is recorded. It is natural to suppose that this costs a little accuracy everywhere and that faster ticks cost less. The truth is more abrupt in both directions. Below a certain rate of ticking, whole waves become indistinguishable from other waves, and no amount of care recovers them. Above it, nothing is lost at all — the wave between the ticks is determined exactly by the values at them.

The figure below is the first half of that statement. A sine wave of nine cycles a unit of time is read ten times a unit, and through exactly the same ten readings runs a sine of one cycle a unit, upside down. The readings do not prefer either wave.

Two waves that agree at every tick #

The coincidence is one line of trigonometry. At the tick $t = k/10$ the fast wave reads $\sin(2\pi \cdot 9k/10)$, and since $9/10 = 1 - 1/10$ this is $\sin(2\pi k - 2\pi k/10) = -\sin(2\pi k/10)$, which is the slow wave’s reading, sign reversed. A whole number of extra turns per tick is invisible, because a sine sampled once a turn sees the same value every time.

The general rule is that a sine of frequency $f$, read $f_s$ times a unit, agrees at every tick with a sine of frequency $f - m f_s$ for any whole number $m$, and with its own mirror image $-f$. The readings therefore see only one representative of the whole family: the member closest to nought, which lies between nought and $f_s/2$. Every faster member is an alias of it, a different wave answering to the same name.

That is the same arithmetic that makes a sequence of complex numbers on the unit circle repeat. Sampling a sine of frequency $f$ at rate $f_s$ is stepping round a circle by the angle $2\pi f/f_s$ each tick, and a step of $2\pi(f/f_s + m)$ lands in the same place. Frequencies differing by a multiple of the sampling rate are the same rotation, and frequencies of opposite sign are the same rotation run backwards — which, for a real wave, is the same wave.

The frequency axis folded #

Draw the frequency a sampled wave seems to have against the frequency it really has, and the graph is a zig-zag.

Up to half the sampling rate — the Nyquist frequency, after Harry Nyquist, whose 1928 study of telegraph signalling counted how many independent pulses a channel of given bandwidth could carry — the samples report the true frequency. Beyond it the axis folds back, and folds again at the full rate, and again at one and a half times, like a strip of paper concertinaed at regular intervals. Every frequency on the strip lies directly above or below some frequency in the first panel, and the samples cannot see the difference between them.

The fold is not a defect of a particular method of reading samples. It is a fact about the samples themselves: a list of numbers taken ten times a unit is consistent with infinitely many waves, and the folding diagram lists them. Any procedure that turns samples back into a wave has to choose one, and the natural choice — the slowest — is right only when the wave had nothing faster than half the rate to begin with.

Orthogonal only up to an alias #

The fold has a second, algebraic face, and it is the one a computer meets. A Fourier coefficient is found by multiplying by one harmonic and taking the area, and the method works because two different harmonics multiply to something with zero area: over a whole period, $\int \cos(2\pi m t)\cos(2\pi n t)\,dt$ is nought unless $m = n$. With only $N$ samples a period, the area becomes a sum over the samples, and the sum is

$$\sum_{k=0}^{N-1} e^{2\pi i (m - n) k/N} = \begin{cases} N & \text{if } N \text{ divides } m - n, \\ 0 & \text{otherwise.} \end{cases}$$

The second line is the fact that the $N$-th roots of unity add up to nothing. The first line is aliasing. Two harmonics whose frequencies differ by a multiple of $N$ are not orthogonal on the samples — they are identical there — so the sampled coefficient of harmonic $m$ is the sum of the true coefficients of $m$, $m \pm N$, $m \pm 2N$ and so on. A discrete Fourier transform of $N$ samples therefore reports $N$ numbers, each of which is a whole family of true coefficients added together, and the fast Fourier transform that computes them so efficiently inherits the same bookkeeping. The coefficients for frequencies above $N/2$ are printed as negative frequencies, because that is the member of the family closest to nought.

When the wave is band-limited and $N$ is large enough, each family has only one member that is not zero, and the sampled coefficients are the true ones exactly. That is the sampling theorem again, read in the language of series rather than transforms.

One pulse for every sample #

The other half of the statement is the surprise. Suppose the wave contains no frequency above some limit $B$ — it is band-limited — and is sampled faster than $2B$. Then the folding diagram’s first panel holds every frequency the wave has, no two of them alias, and the samples determine the wave exactly. The reconstruction is explicit:

$$x(t) = \sum_{k} x\!\left(\frac{k}{f_s}\right) \operatorname{sinc}(f_s t - k), \qquad \operatorname{sinc}(u) = \frac{\sin \pi u}{\pi u}.$$

Each sample contributes one pulse, centred on its own tick, equal to one there and to nought at every other tick, scaled by the sample’s value.

Edmund Whittaker found this series in 1915 while asking which smooth function best interpolates a table of values, and called it the cardinal function; Vladimir Kotelnikov proved the sampling theorem for engineers in 1933, and Claude Shannon made it the foundation of digital communication in 1949. With the 801 samples between $t = -40$ and $40$, the series rebuilds the four-sine signal to within 0.0007 everywhere on the interval drawn, and what remains is entirely the samples left out: the sinc pulse decays only like $1/t$, so distant samples still matter a little.

The pulse is not new. It is the continuous relative of the Dirichlet kernel, whose ripples make Fourier partial sums run away: both are what a sharp cut-off in frequency looks like in time. And the reason the theorem holds is the passage from series to transform. A wave’s spectrum, as its period grows, becomes a continuous curve; sampling that wave copies the curve at every multiple of $f_s$, and if the copies do not overlap — if $f_s > 2B$ — cutting out the middle copy recovers the original spectrum, and with it the wave. If they overlap, the overlapping tails are added together and no cut can separate them. Aliasing is overlapping copies, and the folding diagram is their bookkeeping.

A cliff, not a slope #

The theorem predicts more than good reconstruction above the rate $2B$. It predicts that the quality does not degrade gracefully below it.

Below 7.4 samples a unit the fastest sine folds onto a slower frequency and the reconstruction is simply wrong, by about the size of that sine, at every rate; the error barely changes as the rate climbs towards the threshold. Just above it the error drops by two orders of magnitude, to the level set by truncating the series. There is no middle ground. A rate 2 per cent too low is as bad as one 40 per cent too low, because in both the fast component has been replaced by a different wave, and the replacement is equally wrong whether it lands near the right frequency or far from it.

This is why every device that digitises a signal filters it first. Compact-disc audio samples at 44,100 times a second, a little over twice the twenty thousand cycles a second that human hearing reaches, and before sampling the sound is passed through a filter that removes everything above that limit. Without the filter, a sound at 30,000 cycles — inaudible in the room — would fold to 14,100 and become audible in the recording. The filter does not make the sampling more accurate; it makes the signal band-limited, which is the theorem’s hypothesis.

The fold can also be used on purpose. A radio signal occupying a narrow band of frequencies — say between 100 and 101 million cycles a second — does not need two hundred million samples a second. Sampled at a little over two million a second, its band folds down intact to near nought, because no other part of the signal occupies the frequencies that fold onto the same place. The rate the theorem really demands is twice the width of the band a signal occupies, not twice its highest frequency, provided the band sits where its copies do not collide; radio receivers that digitise directly at high frequency rely on exactly this deliberate aliasing.

A square wave cannot be sampled cleanly #

The square wave is the standard example of a signal that is not band-limited. It is built from round waves at every odd multiple of its own frequency, with amplitudes falling only like one over the multiple, so it has content at every frequency however high. Sample it, at any rate, and every harmonic above half the rate folds back onto some frequency below it.

The result is audible. An early digital synthesiser that produced a square wave by switching between two values at the sample instants produced, besides the harmonics a square wave should have, a haze of aliased tones at frequencies with no simple relation to the note being played — the third harmonic of a high note folding to an inharmonic pitch, the fifth folding somewhere else. For a note at 1,000 cycles a second sampled at 44,100, the harmonic at 23,000 folds to 21,100 and the one at 25,000 to 19,100, which is within hearing. Modern synthesisers avoid it by generating only the harmonics that fit below half the rate, or by rounding off the wave’s corners so that its high harmonics are already small before the sampler sees them.

The plucked string is the physical counterpart. A string released from a sharp pluck keeps its corner because no harmonic decays, and its spectrum, like the square wave’s, falls off only as a power. A recording of it is band-limited only because the microphone, the air and the filter before the sampler round the corner off — and whatever was above half the rate before them has either been removed or folded into the recording.

A wheel under a strobe #

Film samples time at twenty-four frames a second, and a rotating wheel is a periodic signal in disguise.

A wheel with eight identical spokes looks the same after any turn of 45°, so its appearance is periodic in the angle with period 45°, and the camera samples that angle once a frame. At three turns a second the wheel turns exactly 45° per frame and the film shows it standing still; at 2.9 turns a second it turns a little less than a spoke gap per frame, which the eye reads as a small turn backwards; at 4.2 it turns 63° per frame, which looks like 18° forwards, a stately 1.2 turns a second. The folding diagram applies unchanged, with the spoke gap in place of the period and the frame rate in place of the sampling rate.

A sampled rotation by an irrational fraction of a turn never repeats, and the gaps between its points come in at most three sizes; a sampled rotation by a rational fraction repeats, and the wagon wheel is the case where the fraction is close to a whole number of spokes. The strobe lamps that mechanics use to read the speed of an engine work the same way: they flash at a known rate and are adjusted until the rotating part seems to stand still.

Rings the grid invents #

In two dimensions the fold happens in two directions at once, and a fine pattern read on a coarse grid acquires features it never had.

The pattern $\cos(260(x^2 + y^2))$ has rings whose local frequency grows in proportion to the distance from the centre, so it sweeps through every frequency, and the grid’s Nyquist frequency is crossed on a circle. Inside that circle the coarse picture is a fair copy. Outside, each local frequency folds, and the grid shows new sets of rings centred where the folded frequency passes through nought — at points that have no special meaning in the original pattern at all.

These are moiré patterns, and they are a practical nuisance: a striped shirt photographed by a digital camera can acquire swirling bands that were never in the fabric. Many cameras put a thin optical filter in front of the sensor, blurring the image just enough to remove detail finer than the pixels can carry — the two-dimensional version of the audio filter. The same folding is the whole mechanism of error in integration by lattice rules, where averaging a function over a lattice of points is exact for every frequency except the ones that alias onto nought, and a good lattice is one on which few important frequencies do.

What the samples cannot settle #

The sampling theorem assumes a signal with no frequency above $B$, and that assumption is never exactly true of anything that starts and stops. A wave confined to a finite stretch of time always has some content at arbitrarily high frequencies — a function and its transform cannot both vanish outside bounded regions, the same tension that makes a heat profile’s sharp corners vanish at once — so every real signal is only approximately band-limited, and every real reconstruction carries a small aliased residue. The filter before the sampler decides how small.

The figures here also show only the ideal case. Every sample is exact, and the reconstruction uses hundreds of them. With noise in the samples, the cardinal series amplifies some of it, and with a finite record the slow decay of the sinc pulse leaves errors at the ends. Neither is visible in the plots, which draw the middle of a long, clean record.

The theorem’s rate is also a worst case. It guarantees recovery of every signal band-limited to $B$, and for signals with more structure — a handful of frequencies, unknown but few — fewer samples taken at random instants suffice, a discovery of Emmanuel Candès, Justin Romberg, Terence Tao and David Donoho around 2006 that goes under the name of compressed sensing. It does not contradict the sampling theorem, because it gives up on the signals that are not sparse.

And no figure can show the general statement that a band-limited signal is uniquely determined by its samples. The plots show one signal rebuilt; the uniqueness is the theorem’s content, and it rests on the spectral copies not overlapping, which is a fact about the transform rather than anything a reconstruction can display.

Still open: the best grid in many dimensions #

In two dimensions, a square grid is not the most economical way to sample a pattern whose frequencies fill a disc. The spectral copies sit on the grid’s dual lattice, and the sampling is free of aliasing when copies of the disc placed on that lattice do not overlap; the most economical grid is the one whose dual packs discs most densely, and that is the hexagonal grid, which needs about 13 per cent fewer samples than the square one. Daniel Petersen and David Middleton worked this out in 1962 for any number of dimensions: the best sampling lattice for signals band-limited to a ball is the dual of the densest lattice packing of balls.

That turns a question about sampling into one of the oldest questions in geometry. The densest lattice packings are known in dimensions one to eight and in twenty-four, where the eight-dimensional lattice of the octonions and its twenty-four-dimensional relative give the answers. In every other dimension the densest lattice packing of balls is unknown, and so is the most economical grid for sampling a band-limited signal of that many variables. The problem is genuinely practical — medical scanners sample in three dimensions plus time — and its answer waits on a packing problem that has resisted a century of work.

Folding as the whole story #

A list of values taken at equal ticks sees every frequency only up to a multiple of the sampling rate and up to sign, so it folds the frequency axis at half the rate. Below the fold the samples are faithful: a signal with nothing faster than half the rate is determined exactly, and rebuilt by one sinc pulse per sample. Above it they are not merely approximate but wrong, assigning fast motion to slower waves — a nine-cycle sine to a one-cycle one, a fast wheel to a slow backward one, fine rings to rings that were never there. The threshold is a cliff, because the error comes from a fold rather than from a blur.