Sampling and Reconstruction

3.1  The bridge: two different approximations, one signal chain

Sensors such as microphones, strain gauges, antennas and thermocouples produce a continuous-time signal \(x_c(t)\). A computer cannot store a value at every instant or represent every possible amplitude exactly. The signal must therefore be converted into a sequence of numbers.

This conversion has two separate stages. Understanding their difference is central to the lecture:

Students routinely blur the two into “digitising”. Keeping them apart explains almost every practical decision an engineer makes here: you choose \(f_s\) to satisfy a hard constraint, and you choose \(R\) to hit a noise budget.

3.2  Ideal sampling: multiplication by an impulse train

To analyse sampling we need a model of it that the Lecture 2 toolkit can act on. The trick, and it is only a trick, is to write sampling as a multiplication. Define the impulse train (the “Dirac comb”)

and form the product with the signal. The sifting property of Lecture 1 (\(x(t)\delta(t-t_0)=x(t_0)\delta(t-t_0)\)) then does the rest, turning a product of two awkward objects into a tidy weighted train:

Each impulse carries an area equal to the signal's value at that instant. Two cautions, both worth stating because the notation invites confusion.

• \(x_s(t)\) is still continuous-time. It is a function of \(t\) that happens to be zero except at multiples of \(T\). This is exactly what makes it useful: the CTFT of Lecture 2 can be applied to it, whereas it cannot be applied to a sequence.
• \(x[n]\) is not \(x_s(t)\). The sequence \(x[n]=x_c(nT)\) is indexed on an integer and has been time-normalised: it has forgotten \(T\) entirely. The spacing between consecutive values of \(x[n]\) is always 1, never \(T\) seconds. Section 3.6 shows that this forgetting is exactly what rescales the frequency axis.

This is an ideal mathematical model, not a literal circuit. Real samplers use a switch with a short but finite aperture and hold the measured voltage on a capacitor; they do not generate Dirac impulses. The impulse-train model is still valuable because it predicts the correct frequency-domain effects in a simple form.

3.3  What sampling does to a spectrum: replication

Now cash in the model. Eq. (3.2) is a product in time, and Lecture 2's Eq. (2.30) says a product in time is a convolution in frequency. We need one ingredient: the CTFT of the impulse train, which is itself an impulse train, spaced by the sampling frequency \(\Omega_s=2\pi/T\):

Because the sampling impulse train is periodic with period \(T\), its spectrum contains impulses at multiples of \(\Omega_s=2\pi/T\). Convolution with one shifted impulse creates one shifted copy of the original spectrum; convolution with the full train therefore creates repeated copies:

Everything that follows is a consequence of one question: do the copies overlap?

3.4  The Nyquist–Shannon theorem, and aliasing

Assume \(x_c(t)\) is band-limited, so \(X_c(j\Omega)=0\) for \(|\Omega|\ge\Omega_N\). The original spectrum has total width \(2\Omega_N\), and sampling places a copy at every multiple of \(\Omega_s\). The original and its nearest copy do not overlap when

The theorem says more than “the samples are a good approximation.” For a band-limited signal sampled above the Nyquist rate, the samples determine the continuous waveform exactly, including every value between sampling instants. Band-limiting is essential because it restricts how quickly the waveform can change.

3.4.1  Where the factor of two comes from

The factor of two follows directly from the spectral copies in Eq. (3.4). Each copy extends from \(-\Omega_M\) to \(+\Omega_M\) around its centre, so its total width is \(2\Omega_M\).

The replica centred at \(\Omega=0\) extends up to \(+\Omega_M\). Its neighbour, centred at \(\Omega_s\), extends down to \(\Omega_s-\Omega_M\). They stay clear of each other precisely when

The replicas must therefore be separated by at least \(2\Omega_M\). The detailed shape of the spectrum does not matter for this baseband result; only its highest occupied frequency matters.

When the condition holds, the baseband copy sits alone in \(|\Omega|<\Omega_s-\Omega_M\), so an ideal lowpass filter isolates it exactly, and the \(1/T\) in Eq. (3.4) is undone by a gain of \(T\):

Equation (3.7) shows two useful results. First, reconstruction is exact when the Nyquist condition is satisfied. Second, the reconstruction-filter cutoff may lie anywhere between the original spectrum and its nearest copy. This empty region is the guard band. Sampling exactly at the Nyquist rate leaves no guard band, while oversampling makes practical filter design easier by widening it.

3.4.2  Aliasing: where the energy actually goes

“The replicas overlap” is geometrically clear but does not yet tell you what you will hear. Take the simplest possible signal, a cosine at \(\Omega_0\), whose spectrum is one impulse pair at \(\pm\Omega_0\) (Section 2.9.1). Sample it and Eq. (3.4) puts a pair at \(\pm\Omega_0+k\Omega_s\) for every \(k\). Now reconstruct with an ideal low-pass cutting at \(\Omega_s/2\), and ask which line survives:

• \(\Omega_0<\Omega_s/2\). The \(k=0\) pair is inside the cutoff and every other pair is outside. Out comes \(\cos(\Omega_0 t)\). Correct.
• \(\Omega_s/2<\Omega_0<\Omega_s\). Now the \(k=0\) pair is outside the cutoff and gets discarded — but the \(k=-1\) pair, sitting at \(\Omega_0-\Omega_s\), has moved inside it. Out comes a cosine at a frequency that was never in the input:

The high tone has taken on the alias of a lower one. Geometrically the frequency has been folded about \(\Omega_s/2\), which is why \(f_s/2\) is also called the folding frequency: \(f_s/2+\varepsilon\) comes back as \(f_s/2-\varepsilon\). This is the frequency-domain twin of Lecture 2's Problem 2.3 — there, two cosines produced an identical sequence; here we see which of them a reconstructor will pick, and it always picks the one below \(f_s/2\).

MATLAB corner

This script gives three views of sampling. First, it sweeps a tone through the Nyquist frequency and measures the alias that appears in the sampled data. Next, it plays the same folding effect as a chirp. Finally, it uses Eq. (3.17) to reconstruct a band-limited signal and checks that the sinc interpolation passes through every sample.

Goal: plot measured vs. true frequency through and past \(f_s/2\) to see the fold happen numerically, then verify ideal sinc interpolation reconstructs the samples exactly.

%% aliasing_demo.m - Part 1: watch a tone fold back about fs/2  (Eq. 3.6)
fs = 1000;                       % sampling rate, so folding is at 500 Hz
f0 = 100:50:1900;                % sweep the true tone well past fs/2
n  = 0:2047;
fa = zeros(size(f0));
for k = 1:numel(f0)
    x = cos(2*pi*f0(k)/fs*n);    % the SAMPLES - all we ever get
    X = abs(fft(x .* hann(numel(n))'));
    [~,i] = max(X(1:end/2));
    fa(k) = (i-1)*fs/numel(n);   % frequency actually measured
end
figure('Color','w'); plot(f0, fa, 'LineWidth', 1.6); grid on; hold on;
plot(f0, f0, '--'); xline(fs/2, 'r:', 'f_s/2'); xline(fs, 'r:', 'f_s');
xlabel('true frequency (Hz)'); ylabel('measured frequency (Hz)');
title('The triangle wave of aliasing: everything folds into [0, f_s/2]');
% -> rises to 500, folds back to 0 at 1000, rises again... forever.
%% Part 2 - the same thing you can HEAR (needs a sound card)
fsa = 8000; t = 0:1/fsa:3;
chirp_true = chirp(t, 0, 3, 12000);   % sweeps up to 12 kHz > fsa/2 = 4 kHz
soundsc(chirp_true, fsa);
% The pitch rises, hits 4 kHz, then falls again - it cannot exceed fs/2.
%% Part 3 - sinc interpolation really does pass through every sample
T = 1; n = -20:20; x = sinc(n/4);   % a band-limited test sequence
t = linspace(-20, 20, 4000);
xr = zeros(size(t));
for k = 1:numel(n)
    xr = xr + x(k) * sinc((t - n(k)*T)/T);   % Eq. (3.17)
end
figure('Color','w'); plot(t, xr, 'LineWidth', 1.4); hold on; grid on;
stem(n, x, 'filled');
legend('sinc-interpolated x_r(t)', 'samples x[n]');
title('Ideal band-limited interpolation');

3.4.3  The general folding rule, and the strobe

Example 3.2 worked three cases of one wheel by folding each spoke frequency against \(f_s\) by hand. The rule behind all three is worth stating once, because it applies to any periodic motion under any sampler: the observer reconstructs the smallest motion consistent with the samples, so the apparent rate is the true rate folded into \(|f|\le f_s/2\):

where \(\lceil\cdot\rfloor\) rounds to the nearest integer. The consequences are the ones everyone has seen. At \(f_r=f_s\) the wheel advances exactly one full turn between frames and appears frozen. Just below \(f_s\) it falls slightly short each frame, so \(f_a\) is negative and the wheel appears to turn backwards. Just above, it creeps slowly forwards. None of this is an artefact of film or of the eye: it is Eq. (3.8) with a rotation in place of a sinusoid.

One consequence is worth naming, because it is the basis of an instrument and of a recurring mistake. A strobe flashing at \(f_s\) freezes a shaft turning at \(f_r\) whenever \(f_a=0\), which happens not only at \(f_r=f_s\) but at every \(f_r=kf_s\). A stationary image therefore fixes the shaft speed only up to an unknown integer \(k\), and a strobe tachometer read at the first frozen image can be wrong by an exact multiple. The standard procedure — halve the flash rate and check whether the image freezes again — is Eq. (3.9) being used to resolve \(k\).

3.5  Undersampling on purpose: when \(f_s<2f_{\max}\) is legitimate

Section 3.4 is routinely misquoted as "sample at twice the highest frequency". What Eq. (3.6) actually requires is that the replicas do not overlap, and for a signal whose energy occupies a band well away from zero, that is a much weaker requirement than twice the highest frequency. Exploiting the difference is called bandpass sampling, or undersampling, and it is standard practice in radio receivers.

Let the signal occupy \(f_L \le |f| \le f_H\), with bandwidth \(B = f_H - f_L\). Sampling still places replicas at every multiple of \(f_s\), but the band and its mirror image now have gaps between them for the replicas to land in. Non-overlap holds when \(f_s\) is chosen from one of a set of allowed windows:

Each \(n\) is a Nyquist zone: the band folds down to baseband from the \(n\)-th zone up. Taking the largest admissible \(n\) gives the lowest usable rate,

so the true floor is twice the bandwidth, not twice the top frequency, and it is reached only when \(f_H\) happens to be an integer multiple of \(B\). Note that \(n=1\) recovers ordinary baseband sampling, \(f_s \ge 2f_H\), so Eq. (3.10) contains the Nyquist rule rather than contradicting it.

3.6  The exchange rate: \(\omega=\Omega T\), and the debt from Lecture 2

Lecture 2's glossary promised that Lecture 3 would supply the conversion between continuous-time frequency and discrete-time frequency. It falls out of two lines. Take the CTFT of the impulse train in Eq. (3.2), using the sifting property to evaluate each integral:

Set that beside the DTFT of the sequence, Eq. (2.40), \(X(e^{j\omega})=\sum_n x[n]e^{-j\omega n}\). With \(x[n]=x_c(nT)\), the two sums are term by term identical the moment we identify \(\omega\) with \(\Omega T\):

The CTFT of the sampled signal and the DTFT of the sequence are the same function, drawn on two differently labelled axes. Substituting into Eq. (3.4) gives the relation in the form you will use most often:

Once the signal is written as \(x[n]\), the physical time interval \(T\) is no longer visible on the sample index. The frequency axis must therefore be normalized by the sampling rate. This gives the conversion

The sampling frequency \(f_s\) maps to \(2\pi\) rad/sample, so the spectral copies repeat every \(2\pi\) on the discrete-time axis. The Nyquist frequency \(f_s/2\) maps to \(\omega=\pi\). Thus, DTFT periodicity and the folding of analog frequencies around \(f_s/2\) are two descriptions of the same effect.

3.7  Reconstruction: the ideal filter and sinc interpolation

Figure 3.4.a(c) already contains the recipe: if the replicas do not overlap, keep the one at \(k=0\) and throw the rest away. That is a low-pass filter. Its gain must be \(T\) to undo the \(1/T\) in Eq. (3.4), and its cutoff \(\Omega_c\) may sit anywhere in the gap:

The usual choice is \(\Omega_c=\Omega_s/2=\pi/T\), which is legal for any rate satisfying Eq. (3.5). Now translate to the time domain. A rectangle in frequency inverts to a sinc in time — Lecture 2's Section 2.9.1, read right to left — so

Filtering the impulse train of Eq. (3.2) means convolving with \(h_r\), and convolving an impulse train just drops a copy of \(h_r\) at each impulse, weighted by its area. The result is the Whittaker–Shannon interpolation formula:

At a sampling instant \(t=mT\), the sinc for sample \(m\) equals 1 and every other sinc equals 0. The reconstructed curve therefore passes exactly through every sample. Between sample instants, many sinc terms contribute. The Nyquist condition determines whether those in-between values reproduce the original signal or an aliased alternative.

3.8  What real hardware does instead: AAF, ZOH, oversampling

3.8.1  The anti-aliasing filter, and why it must come first

Real signals are never perfectly band-limited. Noise and interference can extend well above the intended signal band, and these components will fold into the measurement when sampled.

An analog low-pass anti-aliasing filter is therefore placed before the sampler to remove unwanted high-frequency content. It cannot be replaced by a digital filter after the ADC, because aliasing has already mixed the frequencies by that stage and the original components can no longer be separated.

3.8.2  Oversampling: buying filter margin with sample rate

A real anti-aliasing filter needs a finite transition band between frequencies it passes and rejects. Sampling exactly at the Nyquist rate leaves no transition band and would require an unrealistically sharp analog filter.

Oversampling uses a sampling rate well above \(2f_N\). The spectral copies move farther apart, leaving more frequency space for a gradual analog roll-off. A simpler analog filter can then be followed by a sharp digital filter and decimation. This is why practical audio converters sample above 40 kHz for a 20 kHz audio band, and why sigma–delta converters often use much higher internal rates.

3.8.3  The zero-order hold, and its sinc droop

At the output end, a practical DAC does not emit impulses either. It emits the simplest thing a latch and a switch can produce: it holds each sample constant for a full period. This is the zero-order hold,

which is exactly Eq. (3.17) with the sinc kernel replaced by a rectangle — the crudest possible interpolator, a staircase. Its cost is immediate from the \(\mathrm{rect}\leftrightarrow\mathrm{sinc}\) pair: the ZOH is an LTI filter whose frequency response is a sinc,

The zero-order hold has two main effects. Its sinc-shaped magnitude response falls with frequency, reaching \(2/\pi=0.637\), or \(-3.9\,\text{dB}\), at \(f_s/2\); practical DAC systems often compensate for this droop. Its phase term \(e^{-j\Omega T/2}\) represents a delay of half a sample and does not distort the waveform shape. The sinc zeros at multiples of \(f_s\) also help suppress the spectral images.

3.9  Quantization: the loss that cannot be undone

Sampling handled the time axis. The amplitude axis is still continuous, and Section 3.1 warned that this half of the problem has no Nyquist theorem to rescue it. An \(R\)-bit quantizer maps each sample to the nearest of \(2^R\) levels spanning \([-X_m,X_m]\), so the step size is

Quantization is a non-linear, memoryless operation, which is bad news: none of the machinery of Lectures 1–2 applies to it. The standard escape is to stop treating it as an operation at all and model it as an added signal,

and to assume \(e[n]\) is white, uniformly distributed, and uncorrelated with \(x[n]\). Under those assumptions the noise power is a one-line integral:

3.9.1  SQNR and the 6 dB-per-bit rule

Take a signal that uses the converter's full range, uniformly distributed on \([-X_m,X_m]\). Its power is \(P_x=\frac{1}{2X_m}\int_{-X_m}^{X_m}x^2dx=X_m^2/3\), and dividing by Eq. (3.22) with \(\Delta\) from Eq. (3.20):

Each additional bit improves the ideal SQNR by about 6 dB. This gives approximately 48 dB for 8 bits, 72 dB for 12 bits, 96 dB for 16 bits and 144 dB for 24 bits. These values assume that the input uses the converter's full range. If the signal amplitude is halved, the SQNR falls by 6 dB, so correct input gain is important as well as bit depth.

MATLAB corner

Equation (3.24) predicts how SQNR changes with word length. The script first tests that prediction using a signal that spans many quantizer levels. It then uses a low-level sinusoid to show when the usual additive-noise model becomes unreliable, and how dither changes the error.

Goal: measure SQNR vs. word length against the \(6.02R\) dB line, then show a low-level sinusoid's quantization error collapse into audible tones without dither.

%% quantization_sqnr.m - measure the 6.02R dB rule instead of trusting it  (Eq. 3.24)
N = 2^16;  Xm = 1;
x = Xm * (2*rand(1,N) - 1);      % uniform on [-Xm, Xm], as Eq. (3.23) assumes
bits = 4:2:16;  sqnr = zeros(size(bits));
for k = 1:numel(bits)
    R = bits(k);
    D = Xm / 2^(R-1);                        % step size, Eq. (3.20)
    xq = D * round(x/D);                     % the quantizer
    e  = xq - x;
    sqnr(k) = 10*log10(mean(x.^2)/mean(e.^2));
    fprintf('R=%2d  measured %6.2f dB   var(e)=%.3e  D^2/12=%.3e\n', ...
            R, sqnr(k), var(e), D^2/12);
end
figure('Color','w'); plot(bits, sqnr, 'o-', 'LineWidth',1.5); hold on; grid on;
plot(bits, 6.02*bits, '--');
xlabel('word length R (bits)'); ylabel('SQNR (dB)');
legend('measured', '6.02R dB', 'Location','northwest');
% var(e) matches D^2/12 to ~1%: the model of Eq. (3.22) is sound HERE.
%% Now break the model: quantize a slow, low-level sinusoid
n = 0:N-1;  R = 6;  D = Xm/2^(R-1);
xs  = 0.02 * sin(2*pi*n/512);            % tiny signal, spans few levels
es  = D*round(xs/D) - xs;
xd  = xs + (rand(1,N)-0.5)*D;            % + 1 LSB of dither
ed  = D*round(xd/D) - xs;
figure('Color','w');
subplot(2,1,1); plot(20*log10(abs(fft(es.*hann(N)'))/N)); ylim([-160 -20]);
  title('No dither: error spectrum is a comb of harmonics (audible tones)');
subplot(2,1,2); plot(20*log10(abs(fft(ed.*hann(N)'))/N)); ylim([-160 -20]);
  title('With dither: a flat noise floor, slightly higher but benign');

3.10  Discrete-time processing of continuous-time signals

In a practical system, sampling and reconstruction appear on either side of a digital filter. We can therefore ask what equivalent continuous-time response the complete ADC–DSP–DAC chain presents to the original analog signal.

Take the chain one stage at a time. The converter samples, and Section 3.6's exchange rate \(\omega=\Omega T\) maps the continuous spectrum onto the discrete one:

The digital filter multiplies, as any LTI system does in its own frequency domain:

The reconstructor then keeps only \(|\Omega|<\pi/T\) and applies the gain \(T\) of Eq. (3.7). If \(x_c(t)\) was band-limited to \(|\Omega|<\pi/T\) in the first place — that is, if it was sampled above the Nyquist rate — only the \(k=0\) term of Eq. (3.25) survives inside the passband, the \(1/T\) and the \(T\) cancel, and the chain collapses to a single continuous-time system:

3.11  Sampling a sequence: the discrete-time sampling theorem

Sampling is not confined to the boundary between the analogue and digital worlds. A sequence can be sampled too — keep every \(N\)-th value, set the rest to zero — and the result obeys a theorem with the same shape as Section 3.4's, in the same way and for the same reason.

Define \(x_p[n] = x[n]\) when \(n\) is a multiple of \(N\), and \(x_p[n]=0\) otherwise. This is multiplication by a discrete impulse train, so by the windowing property of Lecture 2 it convolves the spectrum with that train's transform, producing \(N\) replicas within one \(2\pi\) period:

The replicas are spaced \(2\pi/N\) apart instead of \(2\pi\). If \(x[n]\) is band-limited in the discrete sense, meaning \(X(e^{j\omega})=0\) for \(\omega_M<|\omega|\le\pi\), they stay clear of each other exactly when

and \(x[n]\) can then be recovered from \(x_p[n]\) by a discrete-time lowpass filter of gain \(N\) and cutoff \(\pi/N\). The argument is Eq. (3.6) transposed: same replication, same non-overlap condition, same ideal filter, with \(2\pi\) in the role of \(\Omega_s\).

3.12  Changing the rate: decimation and interpolation

Digital systems change sampling rate by removing samples or inserting zeros. The index operations are simple, but each one changes the spectrum and therefore requires filtering.

3.12.1  Decimation: keeping one sample in M

The decimator forms \(x_d[n]=x[nM]\), discarding the rest. Compressing a sequence in time stretches its spectrum, so the replicas of Eq. (3.28) spread out by a factor \(M\):

After downsampling by \(M\), the spectral copies become \(M\) times wider. They overlap unless the input is first limited to \(|\omega|<\pi/M\). A complete decimator therefore consists of a low-pass anti-alias filter followed by keeping one sample in every \(M\). Omitting the filter can produce plausible-looking but incorrect output.

3.12.2  Interpolation: inserting L−1 zeros

The expander runs the other way: \(x_e[n]=x[n/L]\) when \(n\) is a multiple of \(L\), zero otherwise. Stretching in time compresses the spectrum, and does so without any loss:

Zero insertion does not lose information, but it creates \(L-1\) extra spectral copies, called images, inside \(|\omega|\le\pi\). A low-pass interpolation filter with cutoff \(\pi/L\) removes these images. Its passband gain is \(L\) to restore the signal level after zero insertion.

% rate_conversion_demo.m - decimation and interpolation, done wrong and right
% (Eq. 3.30-3.32)
fs = 8000; t = (0:1/fs:0.05-1/fs).';
x = sin(2*pi*300*t) + 0.6*sin(2*pi*3200*t);   % 3.2 kHz tone is the trap
M = 4;
% --- wrong: throw samples away with no anti-alias filter
x_bad = x(1:M:end);
% --- right: band-limit to pi/M first, then decimate
h = fir1(96, 1/M);
x_good = filter(h, 1, x); x_good = x_good(1:M:end);
% --- interpolation by L: insert zeros, then remove the images
L = 3; x_up = zeros(L*numel(x_good), 1); x_up(1:L:end) = x_good;
g = L * fir1(96, 1/L);        % gain L restores the amplitude
x_interp = filter(g, 1, x_up);
figure
pwelch(x_bad,  [], [], [], fs/M);  hold on
pwelch(x_good, [], [], [], fs/M);
legend('no anti-alias filter', 'filtered then decimated')