Sampling and Reconstruction

3.1  Introduction

Under certain conditions, a continuous-time signal can be represented completely by its values at points equally spaced in time, and can be recovered from them. This property follows from the sampling theorem. A motion picture is a familiar illustration. It consists of a sequence of frames, each an instantaneous view of a continuously changing scene, and when the frames are shown at a sufficiently high rate we perceive continuous motion. A printed picture provides a spatial analogue, since it consists of a fine grid of dots, each a sample of the continuous scene.

The sampling theorem is also the link between continuous-time and discrete-time signals. Because a suitably band-limited signal can be recovered exactly from its samples, it can be represented by a sequence, processed by a discrete-time system and converted back to continuous time. Discrete-time systems are flexible, inexpensive and easily reproduced, and this is the main reason why so much continuous-time signal processing is now carried out digitally.

Two distinct operations are involved when a physical signal is converted into a sequence of numbers. They have very different properties, and it is important to keep them separate.

In the design of a converter the two operations are therefore treated differently. The sampling rate is chosen to satisfy the conditions of the sampling theorem, while the word length is chosen to meet a specified noise level.

The lecture is organized as follows. Sections 3.2 to 3.5 introduce impulse-train sampling, derive the sampling theorem and examine the effect of undersampling, including the deliberate undersampling of bandpass signals. Section 3.6 establishes the relation \(\omega=\Omega T\) between continuous-time and discrete-time frequency. Sections 3.7 and 3.8 consider reconstruction by interpolation and the holds used in practice, and Section 3.9 treats quantization. Sections 3.10 and 3.11 extend the ideas of sampling to discrete-time signals, leading to decimation, interpolation and rational sample-rate conversion.

3.2  Representation of a signal by its samples: impulse-train sampling

Without further information, a set of equally spaced samples does not determine a signal uniquely. Figure 3.2.a shows three different continuous-time signals that have identical values at integer multiples of \(T\); that is, \(x_1(kT)=x_2(kT)=x_3(kT)\) for every integer \(k\). Clearly, infinitely many signals pass through the same set of samples.

The situation changes if the signal is band-limited, that is, if its Fourier transform is zero outside a finite band of frequencies, and if the samples are taken sufficiently close together in relation to the highest frequency present. The samples then specify the signal uniquely, and it can be reconstructed from them exactly. This is the content of the sampling theorem, which we now develop.

A convenient representation of sampling at regular intervals is obtained by multiplying the signal by a periodic impulse train,

This mechanism is known as impulse-train sampling and is shown in Figure 3.2.b. The impulse train \(p(t)\) is called the sampling function, its period \(T\) the sampling period, and its fundamental frequency \(\Omega_s=2\pi/T\) the sampling frequency. The sampled signal is the product

By the sifting property of the unit impulse (Lecture 1), \(x_c(t)\delta(t-t_0)=x_c(t_0)\delta(t-t_0)\). Applying this to each impulse of \(p(t)\) gives

That is, \(x_p(t)\) is an impulse train whose impulse areas are the values of \(x_c(t)\) at intervals of \(T\).

Two points about this representation should be noted. First, \(x_p(t)\) is still a continuous-time signal: it is a function of \(t\) that is zero except at multiples of \(T\). The continuous-time Fourier transform can therefore be applied to it, and in Section 3.3 we use it to determine the effect of sampling in the frequency domain. Second, \(x_p(t)\) is not the same object as the sequence of samples. The discrete-time sequence

is indexed by the integer \(n\), and consecutive values are separated by one unit of \(n\) rather than by \(T\) seconds. The conversion from \(x_p(t)\) to \(x[n]\) may therefore be regarded as a normalization of the time axis, as Figure 3.2.c illustrates for two sampling periods. We shall see in Section 3.6 that this normalization is responsible for the scaling of the frequency axis between continuous and discrete time.

Impulse-train sampling is a mathematical model. Narrow pulses of large amplitude, which approximate impulses, are difficult to generate and transmit, and practical samplers instead acquire the signal over a short aperture and hold the value on a capacitor. The model is nevertheless the appropriate starting point, because it predicts the frequency-domain consequences of sampling in their simplest form. The hold operation is considered in Section 3.7.

3.3  The effect of sampling in the frequency domain

Since \(x_p(t)\) is the product of \(x_c(t)\) and \(p(t)\), the multiplication property of the Fourier transform (Lecture 2, Eq. (2.31)) gives \(X_p(j\Omega)=\frac{1}{2\pi}X_c(j\Omega)*P(j\Omega)\). The Fourier transform of the periodic impulse train is itself an impulse train, with impulses at the multiples of the sampling frequency:

Convolution with an impulse simply shifts a function, \(X_c(j\Omega)*\delta(\Omega-\Omega_0)=X_c\big(j(\Omega-\Omega_0)\big)\), and it follows that

Thus \(X_p(j\Omega)\) is a periodic function of \(\Omega\), consisting of a superposition of shifted replicas of \(X_c(j\Omega)\) scaled by \(1/T\). This is illustrated in Figure 3.3 for a signal band-limited to \(|\Omega|<\Omega_M\). In part (c), \(\Omega_M<\Omega_s-\Omega_M\), and the replicas do not overlap; the spectrum of \(x_c(t)\) is then reproduced exactly, apart from the factor \(1/T\), around each integer multiple of the sampling frequency. In part (d), \(\Omega_s<2\Omega_M\), and adjacent replicas overlap.

Whether \(x_c(t)\) can be recovered from \(x_p(t)\) depends, therefore, on whether the replicas in Eq. (3.6) overlap. This question leads directly to the sampling theorem.

3.4  The sampling theorem and the effect of undersampling

Let \(x_c(t)\) be band-limited, so that \(X_c(j\Omega)=0\) for \(|\Omega|>\Omega_M\). The replica of Eq. (3.6) centred at the origin extends to \(\Omega_M\), and the replica centred at \(\Omega_s\) begins at \(\Omega_s-\Omega_M\). The two do not overlap provided \(\Omega_s-\Omega_M>\Omega_M\), that is,

When this condition holds, \(x_c(t)\) can be recovered from \(x_p(t)\) by a lowpass filter that retains the replica at \(\Omega=0\) and rejects all the others. The filter must have gain \(T\), to cancel the factor \(1/T\) in Eq. (3.6), and a cutoff frequency between \(\Omega_M\) and \(\Omega_s-\Omega_M\):

The complete system for sampling and reconstruction is shown in Figure 3.4.a, together with the spectra at each point. The result may be stated as follows.

Three remarks complete the statement. First, the theorem asserts more than an approximate representation: under its conditions the samples determine every value of \(x_c(t)\), including the values between the sampling instants. Second, the inequality is strict; Example 3.1 shows that sampling a sinusoid at exactly twice its frequency is not sufficient. Third, the band between \(\Omega_M\) and \(\Omega_s-\Omega_M\) contains no signal energy and may be used as the transition band of the reconstruction filter. This guard band has zero width when the signal is sampled at the Nyquist rate, and it widens as the sampling rate is increased.

Ideal filters cannot be realized, and in any practical application the filter of Eq. (3.8) is replaced by a nonideal filter that approximates the desired characteristic closely enough for the problem at hand. Any departure from the ideal characteristic produces some discrepancy between \(x_r(t)\) and \(x_c(t)\), and the acceptable level of this distortion then determines the choice of filter. To emphasize the basic principles, we continue to assume ideal filters in the development that follows.

3.4.1  Aliasing

If \(\Omega_s<2\Omega_M\), the replicas in Eq. (3.6) overlap. Where they overlap their values add, \(X_c(j\Omega)\) is no longer reproduced within \(|\Omega|<\Omega_s/2\), and \(x_c(t)\) can no longer be recovered by lowpass filtering. This effect is called aliasing. The reconstructed signal then differs from \(x_c(t)\), but the two remain equal at the sampling instants. With the cutoff frequency \(\Omega_c=\Omega_s/2\), each term of the reconstruction is zero at every sampling instant except its own (Section 3.7.1), so that for any choice of sampling frequency

Some insight into the relation between \(x_c(t)\) and \(x_r(t)\) when aliasing occurs is obtained from the simple case of a sinusoid,

whose spectrum consists of a pair of impulses at \(\pm\Omega_0\). Figure 3.4.b shows \(X_p(j\Omega)\) for four values of \(\Omega_0\) with \(\Omega_s\) fixed, together with the passband of the reconstruction filter for \(\Omega_c=\Omega_s/2\). To make the movement of the replicas visible, the impulse originating at \(+\Omega_0\) is drawn with a solid line and that originating at \(-\Omega_0\) with a dashed line. For the four cases the reconstructed output is as follows.

When aliasing occurs, the original frequency \(\Omega_0\) takes on the identity of the lower frequency \(\Omega_s-\Omega_0\), so that

As \(\Omega_0\) increases from \(\Omega_s/2\) towards \(\Omega_s\), the output frequency \(\Omega_s-\Omega_0\) decreases, and when \(\Omega_0=\Omega_s\) the reconstructed signal is a constant. This is to be expected, since sampling once per cycle gives equal samples, identical to those obtained from a constant signal. Frequencies above \(\Omega_s/2\) are thus reflected, or folded, about \(\Omega_s/2\), which accounts for the term folding frequency. Figure 3.4.c shows the four cases in the time domain. In each case the lowpass filter fits to the samples a sinusoid whose frequency is less than \(\Omega_s/2\).

3.4.2  Phase reversal, and sampling at exactly the Nyquist rate

Consider next the sinusoid \(x_c(t)=\cos(\Omega_0t+\phi)\). Its spectrum has the same form as before, except that the impulse at \(+\Omega_0\) now has area \(\pi e^{j\phi}\) and the impulse at \(-\Omega_0\) has area \(\pi e^{-j\phi}\). For the cases (b) and (c) of Figure 3.4.b the reconstruction is again \(x_c(t)\). In cases (d) and (e), however, the impulses that fall within the passband originate from the opposite sides of the original spectrum, so that their phases are interchanged. The result is

in which the sign of the phase has been reversed. This phase reversal is responsible for the apparent reversal of rotation in the stroboscopic effect described in Section 3.4.3.

The sampling theorem requires the sampling frequency to be greater than twice the highest frequency in the signal, not merely greater than or equal to it. The following example shows why.

MATLAB corner

The following script illustrates aliasing numerically. It sweeps the frequency of a sampled tone through and beyond the folding frequency and measures the frequency at which the tone appears in the sampled data. It then plays the same effect as an audible chirp, and finally reconstructs a band-limited signal with the interpolation formula of Eq. (3.22) to verify that the result passes through every sample.

Goal: plot the measured frequency against the true frequency through and beyond \(f_s/2\), and verify that ideal band-limited interpolation reproduces the samples exactly.

%% aliasing_demo.m - Part 1: watch a tone fold back about fs/2  (Eq. 3.7)
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.22)
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 stroboscopic effect

The reflection of higher frequencies into lower ones by undersampling is the principle on which the stroboscopic effect is based. Consider the arrangement in Figure 3.4.f, in which a disc with a single radial line rotates at a constant rate and is illuminated by a strobe. The strobe acts as a sampling system, since it illuminates the disc for very brief intervals at a periodic rate. When the strobe frequency is much higher than the rotational frequency, the rotation is perceived correctly. When the strobe frequency becomes less than twice the rotational frequency, the disc appears to rotate more slowly than it does, and, because of the phase reversal of Eq. (3.12), it may appear to rotate in the opposite direction. At one flash per revolution the line appears stationary, since the rotational frequency has been aliased to zero.

In general, the apparent rate of a periodic motion observed by a sampler is obtained by folding the true rate \(f_r\) into the interval \(|f|\le f_s/2\):

A negative value of \(f_a\) corresponds to apparent rotation in the reverse direction. Figure 3.4.g shows three successive strobe positions for three rotational rates. A similar effect is commonly observed in films, in which the wheels of a moving vehicle may appear to rotate slowly or backwards. A film is a sequence of frames taken at a rate of typically 24 frames per second, and this rate is the sampling frequency.

Equation (3.13) also explains a limitation of the stroboscope as an instrument for measuring rotational speed. The image is stationary whenever \(f_a=0\), which occurs not only for \(f_r=f_s\) but for every \(f_r=kf_s\) with \(k\) an integer. A stationary image therefore determines the speed only to within an integer multiple of the flash rate. The ambiguity is resolved by noting that the highest flash rate giving a single stationary image is equal to the rotational rate.

Aliasing is also used deliberately in the sampling oscilloscope, an instrument intended for observing very high-frequency periodic waveforms. The waveform is sampled at a rate well below its frequency, and the aliased waveform, at a much lower frequency, can be displayed with ordinary circuitry. The deliberate undersampling of bandpass signals, considered next, rests on the same principle.

3.5  Sampling of bandpass signals

The sampling theorem as stated in Section 3.4 applies to signals whose spectra extend from zero frequency to \(\Omega_M\). The essential requirement in its derivation, however, is that the replicas in Eq. (3.6) should not overlap. For a signal whose spectrum occupies a band well away from zero frequency, this requirement can be met with a sampling rate much lower than twice the highest frequency. Sampling in this way is called bandpass sampling, or undersampling, and is widely used in radio receivers.

Let the signal occupy \(f_L\le|f|\le f_H\), with bandwidth \(B=f_H-f_L\). Sampling again places replicas at every multiple of \(f_s\), but the positive-frequency band and its mirror image at negative frequencies leave gaps into which the replicas can fall without overlap. This occurs when \(f_s\) lies in one of the intervals

The integer \(n\) identifies the Nyquist zone, of width \(f_s/2\), from which the band is folded into the baseband interval \([0,f_s/2]\). At the end points of each interval the replicas touch, and a practical design keeps away from them. Taking the largest admissible value of \(n\) gives the lowest usable rate,

which equals \(2B\) only when \(f_H\) is an integer multiple of \(B\). For \(n=1\), Eq. (3.14) reduces to the lowpass condition \(f_s\ge2f_H\), so that the bandpass result contains the lowpass result as a special case.

3.6  Relation between continuous-time and discrete-time frequency

We now relate the spectrum of the sequence \(x[n]\) to that of the continuous-time signal from which it was obtained. Section 3.3 gave \(X_p(j\Omega)\) as a sum of replicas of \(X_c(j\Omega)\). Transforming the sampled signal \(x_p(t)=\sum_n x_c(nT)\,\delta(t-nT)\) of Eq. (3.3) term by term, using the fact that the transform of \(\delta(t-nT)\) is \(e^{-j\Omega nT}\), gives a second expression for the same spectrum, in terms of the sample values:

Eqs. (3.6) and (3.16) therefore describe the same \(X_p(j\Omega)\): the first in terms of the spectrum \(X_c(j\Omega)\), the second in terms of the samples \(x_c(nT)\). The DTFT of \(x[n]\) is \(X(e^{j\omega})=\sum_n x[n]e^{-j\omega n}\) (Lecture 2, Eq. (2.41)). With \(x[n]=x_c(nT)\), the two sums are identical term by term when \(\omega\) is identified with \(\Omega T\). Consequently,

Substituting \(\Omega=\omega/T\) into the replica form of Eq. (3.6), \(X_p(j\Omega)=\frac{1}{T}\sum_k X_c\big(j(\Omega-k\Omega_s)\big)\) with \(\Omega_s=2\pi/T\), expresses the DTFT of the sequence directly in terms of the spectrum of \(x_c(t)\):

The relations among \(X_c(j\Omega)\), \(X_p(j\Omega)\) and \(X(e^{j\omega})\) are illustrated in Figure 3.6 for two sampling rates. The spectrum of the sequence is obtained from that of \(x_c(t)\) by periodic replication at multiples of \(\Omega_s\), Eq. (3.6), followed by the linear scaling \(\omega=\Omega T\) of the frequency axis, Eq. (3.17). The replication is a consequence of impulse-train sampling. The frequency scaling may be regarded as a consequence of the normalization of time in the conversion from \(x_p(t)\) to \(x[n]\) shown in Figure 3.2.c: by the time-scaling property of the Fourier transform, a scaling of the time axis by \(1/T\) produces a scaling of the frequency axis by \(T\).

Under the correspondence \(\omega=\Omega T\), the following frequencies correspond:

Since the sampling frequency corresponds to \(\omega=2\pi\), the replicas \(X_c\big(j(\Omega-k\Omega_s)\big)\) of Eq. (3.6) repeat with period \(2\pi\) on the discrete-time frequency axis, and the folding frequency corresponds to \(\omega=\pi\). The periodicity of the DTFT and the folding of continuous-time frequencies about \(\Omega_s/2\) are therefore two descriptions of the same effect.

3.7  Reconstruction of a signal from its samples using interpolation

Interpolation, the fitting of a continuous signal to a set of sample values, is the general procedure by which a signal is reconstructed from its samples, either exactly or approximately. In the system of Figure 3.4.a, the reconstructed signal is the output of a filter with impulse response \(h(t)\) whose input is \(x_p(t)\). With \(x_p(t)\) given by Eq. (3.3),

Equation (3.19) describes how a continuous curve is fitted between the sample points, and different choices of \(h(t)\) give different interpolation procedures. We consider three of them: band-limited interpolation, which is exact for band-limited signals sampled above the Nyquist rate, and the zero-order and first-order holds, which are simpler approximations.

3.7.1  Band-limited interpolation

For the ideal lowpass filter of Eq. (3.8), \(H_r(j\Omega)=T\) for \(|\Omega|<\Omega_c\) and \(0\) otherwise, with \(\Omega_M<\Omega_c<\Omega_s-\Omega_M\), the impulse response is

so that the interpolation formula becomes

With the usual choice \(\Omega_c=\Omega_s/2=\pi/T\), the factor \(\Omega_cT/\pi\) is equal to one, and Eq. (3.21) takes the form known as the Whittaker–Shannon interpolation formula:

In this case \(h_r(0)=1\) and \(h_r(nT)=0\) for \(n\ne0\). At a sampling instant \(t=mT\), therefore, every term of the sum vanishes except the term \(n=m\), and \(x_r(mT)=x[m]\), which is the property stated in Eq. (3.9). Between the sampling instants all of the terms contribute. The reconstruction is illustrated in Figure 3.7.a. Interpolation with the impulse response of an ideal lowpass filter is called band-limited interpolation, since it reconstructs \(x_c(t)\) exactly whenever \(x_c(t)\) is band-limited and the sampling frequency satisfies the conditions of the sampling theorem.

3.7.2  The zero-order hold

Narrow pulses that approximate impulses are difficult to generate, and it is often more convenient to produce the sampled signal in the form given by a zero-order hold. Such a system samples \(x_c(t)\) at a given instant and holds that value until the next sampling instant, as shown in Figure 3.7.b. The output is a staircase approximation to \(x_c(t)\), and digital-to-analog converters commonly produce an output of this form.

The output \(x_0(t)\) of the zero-order hold can in principle be generated by impulse-train sampling followed by an LTI system with a rectangular impulse response, as depicted in Figure 3.7.c. In the notation of Eq. (3.19),

The frequency response of the zero-order hold follows from the transform of a rectangular pulse together with the time-shifting property:

Its magnitude falls from \(T\) at \(\Omega=0\) to \(2T/\pi\) at \(\Omega=\Omega_s/2\), a reduction by the factor \(2/\pi\), or 3.92 dB, and it has zeros at the nonzero multiples of \(\Omega_s\). The linear-phase factor corresponds to a delay of half a sampling period, which does not distort the waveform. Figure 3.7.d compares \(|H_0(j\Omega)|\) with the magnitude of the ideal interpolating filter. The zero-order hold is evidently a coarse approximation: it attenuates the upper part of the signal band, and it passes a considerable part of each spectral image.

In some applications the output of the zero-order hold is adequate without further filtering, particularly when lowpass filtering is applied naturally at a later stage. Figure 3.7.e illustrates this for an image, sampled and reconstructed in two dimensions. Part (a) is the original image, and part (b) shows its samples taken with spatially narrow pulses at a spacing of 12 pixels in each direction. Part (c) is the result of applying a two-dimensional zero-order hold to (b), which produces a mosaic. The human visual system introduces lowpass filtering with a cutoff frequency that decreases with viewing distance, so that when the image is viewed from a distance the discontinuities of the mosaic are smoothed. Part (d) is obtained with one-fourth of the sample spacing in each direction; with normal viewing the mosaic is then much less evident, although the edges of the lettering and the fence remain visibly stepped.

3.7.3  The first-order hold

The zero-order hold produces a discontinuous output. A smoother reconstruction is obtained by linear interpolation, in which adjacent sample points are joined by straight lines, as in Figure 3.7.f. The result is continuous, although its derivative is discontinuous at the sample points because of the changes in slope.

Linear interpolation, also called the first-order hold, again has the form of Eq. (3.19), with the triangular impulse response

illustrated in Figure 3.7.g. Since \(h_1(t)\) is a rectangular pulse of duration \(T\) convolved with itself and scaled by \(1/T\) (Problem 3.5), its frequency response is the square of a sinc function:

At \(\Omega=\Omega_s/2\) its magnitude is \(4T/\pi^2\), a reduction of 7.84 dB relative to its value \(T\) at zero frequency, and it has double zeros at the nonzero multiples of \(\Omega_s\). Compared with the zero-order hold, the first-order hold attenuates the images more strongly, but it also attenuates the upper part of the signal band more. As defined here \(h_1(t)\) is noncausal; a causal implementation delays it by \(T\), which introduces the linear-phase factor \(e^{-j\Omega T}\) without changing the magnitude.

Figure 3.7.h compares the zero-order and first-order holds applied to the same samples of the test image. The first-order hold, which in two dimensions interpolates linearly between the four nearest samples, removes the blocks of the mosaic, although fine detail such as the lettering and the fence, which the sampling has aliased, is not restored by either hold. In an analogous manner, second- and higher-order holds can be defined that produce reconstructions of greater smoothness. The output of a second-order hold, for example, is continuous and has a continuous first derivative, although its second derivative is discontinuous.

The properties of the three interpolating functions are compared below.

3.8  Practical sampling and reconstruction

3.8.1  The anti-aliasing filter

Physical signals are never exactly band-limited. Noise and interference extend above the band of interest, and any component above the folding frequency is aliased into the band when the signal is sampled. A continuous-time lowpass filter, the anti-aliasing filter, is therefore placed before the sampler to attenuate these components. It cannot be replaced by a discrete-time filter after the converter, because at that point the aliased components already occupy the same frequencies as the wanted signal (Section 3.4.1).

3.8.2  Oversampling

A realizable anti-aliasing filter has a transition band of nonzero width between the frequencies that it passes and those that it rejects. If the signal is sampled at the Nyquist rate, the guard band between the signal band and the first replica has zero width, and the filter would have to be ideal. In oversampling, the sampling rate is chosen well above \(2f_M\), so that the replicas are further apart and the analog filter can have a gradual transition, as shown in Figure 3.8.a. A sharp discrete-time filter, followed by decimation (Section 3.11), then completes the band-limiting. For this reason audio converters sample above 40 kHz for a 20 kHz band, and sigma–delta converters operate internally at much higher rates.

3.8.3  Reconstruction from the output of a zero-order hold

When a digital-to-analog converter produces the output of a zero-order hold, the reconstruction filter must compensate for the frequency response of the hold. Figure 3.8.b shows the representation of the zero-order hold from Figure 3.7.c followed by a reconstruction filter \(H_{\mathrm{eq}}(j\Omega)\). The overall output \(r(t)\) is equal to \(x_c(t)\) if the cascade of \(h_0(t)\) and the reconstruction filter is the ideal lowpass filter of Eq. (3.8), that is, if \(H_0(j\Omega)H_{\mathrm{eq}}(j\Omega)=H_r(j\Omega)\). From Eq. (3.24), this requires

For a cutoff frequency \(\Omega_c=\Omega_s/2\), the magnitude and phase of \(H_{\mathrm{eq}}(j\Omega)\) are shown in Figure 3.8.c, together with those of the hold that it compensates. Within the passband the magnitude rises from 1 at \(\Omega=0\) to \(\pi/2\) at the band edge, compensating the droop of the hold, and the phase increases linearly, cancelling the half-sample delay. The phase advance cannot be realized causally, and it is not needed, since a delay of \(T/2\) does not distort the waveform. The magnitude correction is approximated in practice, frequently by a discrete-time filter with a response proportional to \((\Omega T/2)/\sin(\Omega T/2)\) placed before the converter. In many applications the output of the zero-order hold is used without any correction.

Oversampling reduces the need for compensation. If the hold operates at a rate several times the Nyquist rate, the signal occupies only the lower part of the interval \(|\Omega|<\Omega_s/2\), where the droop is small, as Figure 3.8.d illustrates.

3.9  Quantization

Sampling discretizes the time axis, but the amplitude of each sample is still a real number. A quantizer with a word length of \(R\) bits maps each sample to the nearest of \(2^R\) levels spanning the range \([-X_m,X_m]\), as in Figure 3.9.a, so that the step size is

Quantization is a nonlinear, memoryless operation, and the methods of LTI analysis do not apply to it directly. It is usual, therefore, to represent the quantized value as the sum of the sample and an error signal,

and to model \(e[n]\) as a white sequence, uniformly distributed over one quantization step, as in the bottom panel of Figure 3.9.a, and uncorrelated with \(x[n]\). Under these assumptions the noise power follows in three steps.

Collecting the steps,

3.9.1  Signal-to-quantization-noise ratio

Consider a signal that uses the full range of the converter and is uniformly distributed over \([-X_m,X_m]\). Its power is \(P_x=\frac{1}{2X_m}\int_{-X_m}^{X_m}x^2\,dx=X_m^2/3\), and dividing by Eq. (3.30), with \(\Delta\) from Eq. (3.28), gives

Each additional bit therefore improves the SQNR by approximately 6 dB, giving about 48 dB for 8 bits, 72 dB for 12 bits, 96 dB for 16 bits and 144 dB for 24 bits. For a full-scale sinusoid the corresponding figure is \(6.02R+1.76\) dB. These values assume that the signal fills the range of the converter. If the amplitude of the signal is halved, its power falls by a factor of four and the SQNR by 6 dB, so that the gain preceding the converter is as important as its word length.

MATLAB corner

The script below compares the measured SQNR of a uniform quantizer with Eq. (3.32), using a signal that spans many quantization levels. It then quantizes a low-level sinusoid with and without dither to show the conditions under which the additive-noise model fails.

Goal: measure the SQNR as a function of word length and compare it with the \(6.02R\) dB line, and show that without dither the quantization error of a low-level sinusoid consists of discrete harmonics.

%% quantization_sqnr.m - measure the 6.02R dB rule instead of trusting it  (Eq. 3.32)
N = 2^16;  Xm = 1;
x = Xm * (2*rand(1,N) - 1);      % uniform on [-Xm, Xm], as Eq. (3.31) 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.28)
    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.30) 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  Sampling of discrete-time signals

The analysis developed for the sampling of continuous-time signals has a close counterpart for discrete-time signals, with a number of important applications. In analogy with Figure 3.2.b, a sequence is sampled by multiplying it by a periodic impulse train, as shown in Figure 3.10.a. The sampled sequence \(x_p[n]\) is equal to \(x[n]\) at integer multiples of the sampling period \(N\) and is zero at the intermediate samples:

By the multiplication property of the DTFT, the spectrum of \(x_p[n]\) is the periodic convolution of \(X(e^{j\omega})\) with the transform of \(p[n]\). The latter is an impulse train with impulses at the multiples of the discrete-time sampling frequency \(2\pi/N\),

and it follows that

Equation (3.35) is the discrete-time counterpart of Eq. (3.6) and is illustrated in Figure 3.10.b. The replicas are spaced by \(2\pi/N\). If \(X(e^{j\omega})=0\) for \(\omega_M<|\omega|\le\pi\), they do not overlap provided \(2\pi/N-\omega_M>\omega_M\), that is,

In the absence of aliasing, \(X(e^{j\omega})\) is reproduced around \(\omega=0\) and around the integer multiples of \(2\pi\), and \(x[n]\) can be recovered from \(x_p[n]\) by a lowpass filter with gain \(N\) and a cutoff frequency greater than \(\omega_M\) and less than \(2\pi/N-\omega_M\), as in Figure 3.10.c. If aliasing occurs, the reconstructed sequence \(x_r[n]\) differs from \(x[n]\), but, as in continuous time, the two remain equal at the multiples of the sampling period: \(x_r[kN]=x[kN]\) for every integer \(k\).

With \(h[n]\) the impulse response of the lowpass filter, the reconstruction is \(x_r[n]=x_p[n]*h[n]\). For an ideal filter with cutoff \(\omega_c\), \(h[n]=(N\omega_c/\pi)\sin(\omega_cn)/(\omega_cn)\), and

This is the discrete-time counterpart of the band-limited interpolation formula, Eq. (3.21). In practice the ideal filter is approximated, and discrete-time versions of the zero-order and first-order holds of Section 3.7 are frequently used.

3.11  Changing the sampling rate: decimation, interpolation and rational conversion

In many applications it is inefficient to store or transmit the sampled sequence \(x_p[n]\) directly, since it is known to be zero between the sampling instants. The sequence is instead replaced by one that contains only every \(D\)th value, an operation called decimation. The reverse operation, which increases the sampling rate, is called interpolation. Both are simple operations on the index \(n\), and their effect is understood through the corresponding spectra. In this section \(U\) denotes the upsampling factor and \(D\) the downsampling factor, the sampling period \(N\) of Section 3.10 playing the role of \(D\). A cascade of the two operations changes the sampling rate by the factor \(U/D\). Some texts, and the MATLAB example below, use \(L\) and \(M\) for the two factors.

3.11.1  Decimation

The decimated sequence retains every \(D\)th sample of \(x_p[n]\):

the second form following because \(x_p[n]\) and \(x[n]\) are equal at the multiples of \(D\). The relation among \(x[n]\), \(x_p[n]\) and \(x_b[n]\) is shown in Figure 3.11.a.

The DTFT of the decimated sequence is \(X_b(e^{j\omega})=\sum_kx_b[k]e^{-j\omega k}=\sum_kx_p[kD]e^{-j\omega k}\). Letting \(n=kD\), and noting that \(x_p[n]=0\) when \(n\) is not a multiple of \(D\), the sum may be taken over all \(n\):

Decimation therefore expands the spectrum of \(x_p[n]\) by the factor \(D\), so that the replicas spaced by \(2\pi/D\) in \(X_p(e^{j\omega})\) become replicas spaced by \(2\pi\). Figure 3.11.b shows this for \(D=2\). Combining Eq. (3.39) with Eq. (3.35) expresses the result in terms of the original spectrum:

The decimated sequence is free of aliasing if \(X(e^{j\omega})=0\) for \(|\omega|\ge\pi/D\), which is the condition of Eq. (3.36). Figure 3.11.c shows what happens when it is violated. A decimator therefore consists of a lowpass anti-aliasing filter with cutoff \(\pi/D\) followed by the downsampler. If the input is known to contain empty bands, the filter may use them as transition bands; what matters is that no unwanted occupied band overlaps the wanted band after the frequency axis has been expanded by \(D\).

3.11.2  Interpolation

The upsampler inserts \(U-1\) zeros between successive samples: \(x_u[n]=x[n/U]\) when \(n\) is an integer multiple of \(U\), and \(x_u[n]=0\) otherwise. The original sample values are unchanged; the sequence is expanded in time, and its spectrum is compressed in frequency:

The argument \(U\omega\) compresses each period of the spectrum by the factor \(U\), so that \(U\) copies, called images, appear in each interval of length \(2\pi\). The peak value of the DTFT is not divided by \(U\): a spectrum with unit peak still has images with unit peak, each of them narrower. A lowpass interpolation filter removes the images, and for a full-band input its ideal cutoff is \(\pi/U\).

The interpolation filter normally has passband gain \(U\). A constant input shows why. After zero insertion only one sample in \(U\) is nonzero, so the average value of the expanded sequence is \(1/U\), and a filter with a gain of \(U\) at zero frequency restores the constant output. The three stages are shown in Figure 3.11.d for \(U=2\).

3.11.3  Rational conversion by \(U/D\)

A rational sample-rate converter is the cascade shown in Figure 3.11.e. Its lowpass filter performs two functions: it removes the images created by the upsampler, and it removes any remaining components that would be aliased by the downsampler.

The upsampler raises the sampling rate from \(F_{s,\mathrm{in}}\) to \(UF_{s,\mathrm{in}}\), the filter operates at this intermediate rate without changing it, and the downsampler lowers it by the factor \(D\). The output rate is therefore \(F_{s,\mathrm{out}}=(U/D)\,F_{s,\mathrm{in}}\), and the ratio

is the sampling-rate conversion factor. A factor greater than 1 raises the rate, and a factor less than 1 lowers it. To convert between two given rates, the ratio \(F_{s,\mathrm{out}}/F_{s,\mathrm{in}}\) is reduced to lowest terms and \(U\) and \(D\) are read from its numerator and denominator. A common factor would leave the conversion factor unchanged but raise the intermediate rate, and with it the cost of the filter. For example, conversion from 48 kHz to 36 kHz has \(36/48=3/4\), so that \(U=3\) and \(D=4\).

Suppose that the input spectrum is zero outside \(|\omega|\le\omega_{\max}\). The filter specification is then obtained from the spectra as follows. Figure 3.11.f illustrates each step for the values of Example 3.7.

1. Upsample by \(U\). The wanted copy is compressed to \(|\omega|\le\omega_{\max}/U\), which gives the passband edge
   $$\boxed{\omega_p=\frac{\omega_{\max}}{U}}.$$
2. Locate the nearest image. The images are centred at \(2\pi k/U\), and the nearest positive-frequency image begins at
   $$\boxed{\omega_{\mathrm{st}}=\frac{2\pi-\omega_{\max}}{U}=\frac{2\pi}{U}-\omega_p}.$$
   If this image is the first unwanted occupied band, it defines the stopband edge \(\omega_{\mathrm{st}}\). The interval between \(\omega_p\) and \(\omega_{\mathrm{st}}\) contains no input energy and may be used as the transition band.
3. Complete the lowpass specification. With \(\omega_p\) and \(\omega_{\mathrm{st}}\) known,
   $$\boxed{\Delta\omega=\omega_{\mathrm{st}}-\omega_p,\qquad\omega_c=\frac{\omega_p+\omega_{\mathrm{st}}}{2}}.$$
   The midpoint \(\omega_c\) is the natural cutoff frequency for a window-method FIR design (Lecture 7).
4. Set the gain. A passband gain of \(U\) preserves the amplitude of the signal in the time domain. This gain should be distinguished from the peak value of a plotted DTFT.
5. Check the downsampler. Downsampling expands the wanted band by \(D\), so that
   $$\boxed{\omega_{\max,\mathrm{out}}=D\omega_p=\frac{D}{U}\,\omega_{\max}}.$$
   The output is free of aliasing only if \(D\omega_p<\pi\), and any unwanted occupied band that would be aliased onto the wanted band must already have been attenuated by \(H(z)\).

MATLAB corner

The script below demonstrates decimation with and without an anti-aliasing filter, and interpolation by zero insertion followed by the removal of images. The input, sampled at 8 kHz, contains a wanted 300 Hz tone and an unwanted 3.2 kHz tone.

Goal: show that decimating by \(M=4\) without filtering aliases the 3.2 kHz tone to 800 Hz, and that lowpass filtering to \(\pi/M\) before decimation removes it.

%% rate_conversion_demo.m - decimation and interpolation, with and without filtering
% (Eqs. 3.48-3.50)
fs = 8000; t = (0:1/fs:0.05-1/fs).';
x = sin(2*pi*300*t) + 0.6*sin(2*pi*3200*t);   % the 3.2 kHz tone will alias
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')

?Problems with solutions

Summary and key formulas

The principal results of the lecture are collected below.

• Sampling and quantization. Sampling is reversible under the conditions of the sampling theorem; quantization is not, and its effect is described by the noise power \(\sigma_e^2=\Delta^2/12\).
• Spectral replication. Impulse-train sampling replicates the spectrum of \(x_c(t)\) at intervals of \(\Omega_s\) and scales it by \(1/T\), Eq. (3.6).
• The sampling theorem. A band-limited signal is uniquely determined by its samples if \(\Omega_s>2\Omega_M\), and it is recovered by an ideal lowpass filter of gain \(T\). The inequality is strict, as Example 3.1 shows.
• Aliasing. If \(\Omega_s<2\Omega_M\) the replicas overlap. A sinusoid with \(\Omega_s/2<\Omega_0<\Omega_s\) is reconstructed at \(\Omega_s-\Omega_0\) with its phase reversed, and the stroboscopic effect is a visible example. Aliasing must be prevented by an anti-aliasing filter before sampling.
• Bandpass sampling. A bandpass signal can be sampled at a rate determined by its bandwidth rather than by its highest frequency, provided \(f_s\) lies in one of the intervals of Eq. (3.14).
• Frequency relation. The DTFT of \(x[n]\) is the spectrum of \(x_p(t)\) with the frequency axis scaled by \(T\), so that \(\omega=\Omega T\). The sampling frequency corresponds to \(\omega=2\pi\) and the folding frequency to \(\omega=\pi\).
• Interpolation. Band-limited interpolation is exact but noncausal. The zero-order and first-order holds are simpler approximations, with attenuations of 3.92 dB and 7.84 dB at \(\Omega_s/2\). The droop of the zero-order hold can be compensated, or reduced by oversampling.
• Quantization. Each additional bit improves the SQNR by about 6 dB at full scale. Dither removes the correlation between the error and the signal.
• Sampling of sequences and rate conversion. Sampling a sequence with period \(N\) replicates its spectrum at intervals of \(2\pi/N\), and decimation then expands the spectrum by \(N\). Upsampling compresses the spectrum and creates images, which the interpolation filter removes. In a rational converter the single filter is placed between the upsampler and the downsampler.
• Next, in Lecture 4: the \(z\)-transform. Sampling a complex exponential gives \(z=e^{sT}\), which maps the left half of the \(s\)-plane into the unit circle and the \(j\Omega\)-axis onto it.