Sampling and reconstruction

Sound in the air is continuous: the pressure has a value at every instant. A computer can only store a list of numbers. Sampling turns one into the other, by measuring the signal at regular instants. Every other page on this site works on those lists of numbers, and quietly relies on the results here: what sampling keeps, what it loses, and how the original can be recovered, in some cases exactly.

Sampling measures the signal every $T$ seconds: $f_s = 1/T$ times a second, the sample rate. The result is a list of values $x[n] = x(nT)$, each saying what the signal was at that instant. That’s why the plots on this site draw samples as lollipops: a value at each sample time, and nothing claimed in between.

Joining the dots with straight lines is tempting, and very common in plotting software. But it draws something that isn’t in the data: a guess that the signal moved in straight lines between samples. Real signals don’t. The right way to fill the gaps is quite different. And under certain conditions, which sections 4 and 5 spell out, it recovers the original continuous signal exactly, even though all we kept were values at separate instants.

Films sometimes show a fast-turning wagon wheel spinning slowly, or even backwards. The camera takes a picture, a sample, 24 times a second. If a spoke moves almost a full turn between frames, it looks as if it moved a little backwards. The camera can’t tell “turned 350°” from “turned −10°”.

Samples of a sinusoid have exactly the same problem. A sinusoid at frequency $f$ turns by $2\pi f/f_s$ radians per sample. Add a whole extra turn per sample, which is frequency $f + f_s$, and every sample is identical, since a whole turn ends where it started. So sampling can’t tell $f$ from $f + f_s$, $f + 2f_s$, … or from $f - f_s$. A high frequency masquerades as a low one: it’s aliased.

Every frequency has a stand-in in the range $-f_s/2$ to $+f_s/2$ that produces the same samples. Only frequencies already inside that range come through as themselves. Half the sample rate, $f_s/2$, is called the Nyquist frequency.

In terms of spectra, aliasing has a neat description. Since the samples can’t tell $f$ from $f + f_s$ or $f - f_s$ (or any other whole number of extra turns per sample, forwards or backwards), the spectrum of the samples contains every component of the original at every shifted position $f + k f_s$ at once, for $k = \dots, -2, -1, 0, 1, 2, \dots$ Sampling copies the whole spectrum, once centred on each multiple of the sample rate, on both sides of zero.

The figure draws spectra with negative frequencies too. As the Fourier page showed, a real sinusoid is the sum of two arrows spinning in opposite directions, one at $+f$ and one at $-f$. So a real signal’s spectrum always has a mirror image on the negative side. Both halves get copied, and it’s the copies of the negative half, shifted up by $f_s$, that reach down towards the original’s positive half. That’s the falling alias from section 2: a tone at $f$, just above $f_s/2$, has its mirror image at $-f$. Shifted up by $f_s$, the mirror image lands at $f_s - f$, just below $f_s/2$, and as $f$ rises, $f_s - f$ falls.

The picture suggests something remarkable. If the copies don’t overlap, the original spectrum is still in there, untouched, in the middle. Keep just that copy, by removing everything outside $-f_s/2 \dots f_s/2$ with an ideal low-pass filter, and you have the original signal back: exactly, between the samples too. Tick the reconstruction filter in the figure above to see which part it keeps.

That’s the sampling theorem (Nyquist, Shannon and others, 1920s–1940s):

If a signal contains no frequencies at or above half the sample rate, it is completely determined by its samples, and can be rebuilt from them exactly.

Note the condition: no frequencies at or above $f_s/2$. Such a signal is called band-limited.

What does the ideal low-pass filter do in time? On the FIR page we met its impulse response: the sinc function, $\operatorname{sinc}(t) = \sin(\pi t)/(\pi t)$. So reconstruction means: put a sinc at every sample, scaled by that sample’s value, and add them all up.

Each sinc equals 1 at its own sample and 0 at every other sample instant. So at the sample times, the sum just gives back each sample. Between them, the sincs’ wiggles combine to fill in the curve, and if the signal was band-limited they fill it in perfectly.

The theorem’s condition matters. If a signal contains frequencies above $f_s/2$ when it’s sampled, they alias down into the band, and once they’re there nothing can tell them apart from genuine low frequencies. Aliasing can’t be undone after the fact.

A classic way to fall into this is generating waveforms digitally. A square wave contains every odd harmonic forever, at 1/3, 1/5, 1/7… of the fundamental. Make one the obvious way, as $+1$ or $-1$ depending on the sign of a sine, and all the harmonics above the Nyquist frequency fold back down. Unlike the real harmonics, they land at inharmonic frequencies, unrelated to the note.

Both versions keep a phase: how far round the cycle the wave has got, in radians. Each sample, the phase moves on by $2\pi f / f_s$, the same “turn per sample” as the wheel in section 2. The difference is only in how a sample is made from the phase.

The naive version asks “which half of the cycle are we in?”:

# Naive: jump between +1 and -1. Simple, and aliased.
phase = 0
for n in range(num_samples):
    out[n] = +1 if sin(phase) >= 0 else -1
    phase += 2 * pi * f / fs

The jumps are the problem. A jump is an infinitely sharp edge, which needs harmonics at every odd multiple of $f$, forever. The samples can’t hold the ones above $f_s/2$, so they alias.

The band-limited version follows Fourier’s recipe from the Fourier page instead: add up the odd harmonics, at $1, \tfrac13, \tfrac15, \dots$ of the first, but stop before reaching half the sample rate:

# Band-limited: a sum of odd harmonics, keeping only those below fs/2.
phase = 0
for n in range(num_samples):
    total = 0
    k = 1
    while k * f < fs / 2:          # the line that prevents aliasing
        total += sin(k * phase) / k
        k += 2                     # odd harmonics only: 1, 3, 5, ...
    out[n] = (4 / pi) * total      # scale to run from -1 to +1
    phase += 2 * pi * f / fs

Nothing above $f_s/2$ is ever created, so there’s nothing to fold back. The price is a slight ripple near each edge (the Gibbs overshoot again), and work: a 100 Hz note at 44.1 kHz has 110 harmonics to add up for every sample. Real synthesisers get the same result more cheaply. Some store pre-computed band-limited cycles (wavetables) and pick one according to the note. Others take the naive waveform and smooth off just the samples next to each jump (BLEP, short for band-limited step), which removes most of the aliasing at a fraction of the cost.

Real converters wrap the sampling in two analogue low-pass filters, which are the theorem’s assumptions turned into hardware:

• Before sampling, an anti-aliasing filter removes everything above $f_s/2$, so there’s nothing left to alias. It can’t be perfectly sharp. That’s one reason CDs sample at 44.1 kHz, although we hear only up to about 20 kHz: the gap from 20 to 22.05 kHz is room for the filter’s transition band.
• After the DAC, a reconstruction filter does the job of the sincs. The DAC holds each value until the next one arrives, which produces the staircase. In spectral terms, the staircase still contains the copies of the spectrum at multiples of $f_s$ (slightly softened). The reconstruction filter removes them, leaving the smooth, original signal.

Sharp analogue filters are expensive and imperfect, so modern converters use a trick: oversampling. They sample many times faster than needed, so a gentle analogue filter is enough. Then they do the sharp filtering digitally, with an FIR filter, and reduce the sample rate afterwards. Changing the sample rate digitally (decimation and interpolation) is a topic of its own.

Sampling is half of digitising a signal. The other half is that each sample has to be stored as a number with a fixed number of bits, $B$. That allows $2^B$ possible values, so every sample gets rounded to the nearest one. The rounding error is called quantisation error.

For a busy signal, the rounding error looks like a faint, random noise. Each extra bit halves the size of the steps, so it lowers the noise by 6 dB. For a sine wave that just fills the available range:

So 16 bits (CD) gives about 98 dB, and 24 bits about 146 dB, more than any real analogue circuit can manage.

But for quiet or simple signals the error isn’t random at all. It follows the signal, repeating with it, and that shows up as distortion: harmonics that weren’t there. A quiet note fading away on 8-bit audio crunches and buzzes rather than just getting hissy. The cure sounds paradoxical: dither, adding a tiny amount of random noise (about one step’s worth) before rounding. The noise randomises the rounding, turning the distortion into a gentle, constant hiss that’s unrelated to the signal. It costs about 4.8 dB of noise, but the ear much prefers hiss to distortion, and details smaller than one step remain audible within it.

• Samples are values at instants. Nothing is claimed in between, so draw them as lollipops.
• Sampling can’t tell $f$ from $f \pm f_s$: high frequencies alias down to low ones.
• In the spectrum, sampling makes copies at every multiple of $f_s$.
• If the signal is band-limited (nothing at or above $f_s/2$), the copies don’t overlap, and a sinc per sample rebuilds it exactly.
• If it isn’t, aliasing mixes in for good. Filter before sampling.
• Quantisation adds about 6 dB of noise per bit removed. Dither trades a little more noise for no distortion.

• C. E. Shannon, “Communication in the presence of noise” (1949): the sampling theorem as it’s usually stated, in the paper that also founded information theory.
• Julius O. Smith III, Mathematics of the DFT: sampling theory.
• Dan Lavry, “Sampling Theory for Digital Audio” (2004): a clear engineer’s explanation of why “more samples per cycle” isn’t needed.
• Christopher Montgomery (Xiph.org), Digital Show & Tell: a wonderful hands-on video demonstrating all of this with real equipment.