import numpy as np
import matplotlib.pyplot as pltLecture 3 - Sampling & Reconstruction
Summary
In this lecture, we will continue to refresh understanding on the background material from signals and systems. Specifically, we focus on sampling and reconstruction, multirate signal processing, and related practical aspects for analog-to-digital convertors (ADCs) and digital-to-analog converters (DACs).
Motivation in Course Context
Based upon the modern radio architecuture introduced in the first lecture, we are now focused on the DACs and ADCs that provide us with the interface between the CT and DT worlds.
Outline
- Analog-to-Digital Converter (ADC) Architecture
- Digital-to-Analog Converter (DAC) Architecture
- Sampling and Reconstruction of CT Signals
- Frequency-Domain View of Sampling
- Ideal Bandlimited Reconstruction and the Sampling Theorem
- Aliasing, Nyquist Zones, and Anti-Alias Filtering
- Decimation and Interpolation of DT Sequences
- Downsampling
- Interpolation
- DT Sampling Theorem
- Practical Interpolation Filters Require Oversampling
Analog-to-Digital Conversion
At a high level, the ADC architecture in Figure 2 works well provided the following:
First, the desired frequency components of the signal \(r(t)\) are sufficiently bandlimited, e.g., maximum frequency component less than \(f_s/(2M)\).
Second, the filters \(h(t)\) and \(g[n]\) are practical approximations to low pass filters (LPFs) that pass the desired frequency content of \(r(t)\) and \(s[n]\), respectively, and significantly attenuate other undesired frequency components. The wider the transition band of these filters, generally the more bandlimited we require \(r(t)\) to be.
Third, the quantizer does not introduce too much distortion, e.g., \(|y-q(y)|<D\) for some small \(D\) for “most” values of the input \(y\).
Note that quantization issues will be addressed in our next lecture.
Digital-to-Analog Conversion
At a high level, the DAC architecture in Figure 3 works well provided the following:
First, the quantizer does not introduce too much distortion, e.g., \(|y-q(y)|<D\) for some small \(D\) for “most” values of the input \(y\).
Second, the input signal \(y[n]\) is sufficiently bandlimited that its own periodic images, or the additional ones created through upsampling to produce \(x[n]\), do not enter into the transition band or passband of the interpolation filters.
Third, the filters \(g[n]\) and \(h(t)\) are practical approximations to low pass filters (LPFs) that pass the desired frequency content of \(y[n]\) and \(s(t)\), respectively, and significantly attenuate other undesired frequency components. The wider the transition band of these filters, generally the more bandlimited we require \(y[n]\) to be.
Note that quantization issues will be addressed in our next lecture.
Sampling and Reconstruction of CT Signals
The time domain relationship for sampling is rather simple: given a CT signal \(s(t)\) and sampling frequency \(f_s > 0\), we define the sampling period as \(T_s = 1 / f_s\) and set
\[ s[n]=\left. s(t) \right|_{t=nT_s}=s(nT_s) \tag{1}\]
The frequency domain relationship between the two samples that is more complicated.
Sampling in the Frequency Domain
Let \(s(t)\leftrightarrow S(f)\) be a Fourer transform pair, and let \(s[n]=s(n T_s)\), \(n\in\mathbb{Z}\) for some sampling period \(T_s > 0\).
The frequency domain relationship between the CT signal and the DT sequence is as follows:
\[ S(e^{j2\pi u}) = \frac{1}{T_s} \sum_{k=-\infty}^{\infty} S(\left (u - k ) / T_s \right) = f_s \sum_{k=-\infty}^{\infty} S(\left (u - k)f_s \right) \tag{2}\]
This relationship demonstrates that sampling introduces three effects:
- Amplitude scaling by \(1/T_s\)
- Periodic replication of the spectrum at integer multiples of the sampling frequency, i.e., replicas of \(S(f)\) shifted by \(kf_s \mid k\in\mathbb{Z}\) and added up
- Scaling of the frequecny axis \(u = f/f_s\)
After we illustrate these effects via an example, we provide two alternative ways to develop the mathematical relationship.
Illustration
In the example below, \(s(t)\) is a signal bandlimited to maximum frequency \(f_m > 0\) and with Fourier transform \(S(f)\).
We see from Figure 4 that there is no overlap between the periodic replicas \(S(f-kf_s)\), \(k\in\mathbb{Z}\) in \(S_p(f)\) if \[ f_s > 2f_m \tag{3}\]
Informal Proof
We define an intermediate CT signal \(s_p(t) = s(t)\cdot p(t)\), where \(p(t)\) is an impulse train
\[ p(t)=\sum_{n=-\infty}^{\infty} \delta(t-nT_s) \]
so that
\[ s_p(t) = \sum_{n=-\infty}^{\infty} s(nTs) p(t-nT_s) \]
In the frequency domain, since \(p(t)\) is periodic with generalized Fourier transform
\[ P(f) = \frac{1}{T_s} \sum_{k=-\infty}^{\infty} \delta(f-k/T_s) \]
convolution in the frequency domain yields
\[ S_p(f) = \frac{1}{T_s} \sum_{k=-\infty}^{\infty} S(f-k/T_s) \]
Finally, scaling the time access between CT and DT by \(1/T_s\) via \(s[n]=s(nT_s)\) corresponds to scaling the frequency axis by \(T_s\) via
\[ S(e^{j2\pi u})=S_p(u/T_s)=\frac{1}{T_s} \sum_{k=-\infty}^{\infty} S((u-k)/T_s) \]
Formal Proof
To establish the frequency-domain relationship, start with the CT inverse Fourier transform relationship \[ s(t)=\int_{-\infty}^{\infty} S(f) e^{-2\pi f t} df \]
so that
\[ s[n] = s(nT_s) = \int_{-\infty}^{\infty} S(f) e^{-2\pi (f/f_s) n} df \]
The change of variables \(u=f/f_s\) yields
\[ s[n] = f_s \int_{-\infty}^{\infty} S(u f_s ) e^{-2\pi u n} du \]
which shows \(s[n]\) as a superposition of an infinite continuum of DT complex exponentials. However, these complex exponentials only live on the finite normalized frequency axis \(u\in[-1/2,1/2]\) because \(e^{-2\pi (u+k) n}=e^{-2\pi u n}\) for all \(k\in\mathbb{Z}\).
We therefore break up the infinite integral into a sum of integrals, one over each unit interval \([k-1/2,k+1/2]\), \(k\in\mathbb{Z}\). Specifically,
\[ s[n] = f_s \sum_{k=-\infty}^{\infty} \int_{k-1/2}^{k+1/2} S(u f_s ) e^{-2\pi u n} du \]
Finally, we define the changes of variables \(v=u-k\) and \(l=-k\) so that
\[ s[n] = f_s \sum_{k=-\infty}^{\infty} \int_{-1/2}^{+1/2} S((v+k) f_s ) e^{-2\pi (v+k) n} dv = \int_{-1/2}^{+1/2} \left[ f_s \sum_{l=-\infty}^{\infty} S((v-l) f_s ) \right] e^{-2\pi v n} dv \]
The last relationship is clearly the inverse DT Fourier transform relationship for \(s[n]\), so that by inspection we have the desired result
\[ S(e^{2\pi v n})=f_s \sum_{l=-\infty}^{\infty} S((v-l) f_s ) \tag{4}\]
Bandlimited Interpolation and the Sampling Theorem
We say that a signal \(s(t)\) is bandlimited to maximum frequency \(f_m\) if \(S(f)=0\) for \(|f|>f_m\).
The example and informal proof above suggests that, despite periodic replication through the sampling process, we can recover a bandlimited CT signal from its samples provided there is no overlap of the replicas in the frequency domain, i.e., \(f_s > 2 f_m\). This lower limit on the sampling frequency is called the Nyquist rate for the bandlimited signal.
Bandlimited Interpolation
Exact recovery simply consists of filtering the intermediate signal \(s_p(t)\) with a lowpass filter of the appropriate cutoff frequency and passband gain, and illustrated in the figure below.
The lowpass filter passes and scales the replica centered at \(f=0\) and eliminates the replicas centered at \(f=kf_s\), \(k\in\mathbb{Z}-\{0\}\).
Strictly speaking, the cutoff frequency for such a filter can be any frequency between \(f_m\) and \(f_s-f_m\); however, a convention is to set the cutoff frequency to be \(f_s / 2\), as illustrated in Figure 5. With this convention, the ideal lowpass filter has frequency response
\[ H(f)=T_s \mathrm{rect}(fT_s) \tag{5}\]
and impulse response
\[ h(t)=\mathrm{sinc}(t/T_s) \tag{6}\]
The resulting reconstruction equation is called bandlimited interpolation, and takes the form
\[ \tilde{s}(t)=s_p(t) * h(t) = \left(\sum_{n=-\infty}^{\infty} s(nT_s) \delta(t-nT_s)\right) * \mathrm{sinc}(t/T_s) = \sum_{n=-\infty}^{\infty} s(nT_s) \mathrm{sinc}\left( \frac{t-nT_s}{T_s} \right) \tag{7}\]
Sampling Theorem
Putting these observations together leads to the following result.
Theorem 1 (Sampling Theorem) If a continuous signal \(s(t)\) is bandlimited to maximum frequency \(f_m > 0\) and sampled with frequency \(f_s > 2 f_m\), corresponding to the period \(T_s < 1/(2f_m)\), then it can be reconstructed from its \(T_s\)-spaced samples as:
\[ \lim_{N \rightarrow \infty} \sum_{n=-N}^{N} s(nT_s)\ \mathrm{sinc}\left( \frac{t-nT_s}{T_s} \right)=s(t) \]
Aliasing, Nyquist Zones, and Anti-Alias Filtering
If a CT signal \(s(t)\) is not bandlimited, or even if the signal is bandlimited but the sampling frequency \(f_s\) does not exceed its Nyquist rate, the periodic replicas \(S(f-kf_s)\) will overlap in the frequency domain. This effect is called aliasing.
Aliasing can be better understood if we define the Nyquist Zones, which consist of successive intervals of length \(f_s/2\) along the CT frequency axis. The Nyquist Zones are illustrated in the figure below.
Specifically, for positive \(k\in\mathbb{Z}\), the \(k\)-th Nyquist Zone is the interval \([(k-1)f_s/2,kf_s/2]\). For negative \(k\mathbb{Z}\), the \(k\)-th Nyquist Zone is the interval \([kf_s/2,(k+1)f_s/2)\). All Nyquist Zones for \(k\) odd and positive alias into the \(k=1\) Nyquist Zone, and all Nyquist Zones for \(k\) even and positive alias into the \(k=-1\) Nyquist Zone. Conversely, all Nyquist Zones for \(k\) odd and negative alias into the \(k=-1\) Nyquist Zone, and all Nyquist Zones for \(k\) even and negative alias into the \(k=1\) Nyquist Zone.
To avoid aliasing, we often put a lowpass filter in front of the sampler, to limit the frequency content. Ideally, this would be a lowpass filter with cutoff frequency \(f_s/2\). Such a filter is called an anti-alias filter.
Decimation and Interpolation of DT Sequences
In this section, we develop similar ideas for changing the sampling rate of DT sequences as for sampling of CT signals. The concepts are largely the same, although the notation is slightly different, and we have to keep in mind that the DT normalized frequency axis is periodic with period \(1\).
Decimation
In some applications, we have fixed analog hardware that operates at a high sampling rate \(f_s\), and then we reduce the effective sampling rate in reprogrammable digital hardware or software. Decimation or Downsampling by a factor \(M\in\mathbb{Z}^+\) has input-output relationship
\[ y[n]=x[nM] \tag{8}\]
which looks very similar to the relationship for sampling a CT signal to produce a DT sequence. The downsampler keeps every \(M\)-th sample and discards those inbetween.
Figure 7 shows a block diagram for the downsampling operation and illustrates a good way to think about its operation relative to CT sampling.
In particular, if we imagine a CT signal \(x(t)\) for which the sequence \(x[n]\) represents the samples \(x(nT_s)\), downsampling yields the samples \(y[n]=x(nMT_s)\), i.e., samples spaced with larger period \(MT_s\) or smaller sampling frequency \(f_s/M\).
Using the previous frequency-domain relationships for sampling from CT to DT, we can, after some algebraic manipulations, obtain the following DT relationship for downsampling
\[ Y(e^{j2\pi u})=\frac{1}{M} \sum_{m=0}^{M-1} X(e^{j2\pi (u-m)/M}) \tag{9}\]
Again, we have replication of the spectrum, except there are only \(M\) replicas within a period of length \(1\) in the DT frequency axis.
As in the CT case, aliasing can occur, we can define DT Nyquist zones, and we can apply anti-alias or decimation filters before downsampling to reduce the effects of aliasing.
Interpolation
We often want to recover a DT sequence from its downsampled version, or to increase the effective sampling rate in reprogrammable hardware or software, partiularly in order to produce samples for fixed analog hardware that operates at a high sampling rate. This interpolation is accomplised in two steps, upsampling followed by filtering.
Upsampling
Upsampling by a factor \(M\in\mathbb{Z}^+\) has input-output relationship
\[ w[n]= \begin{cases} x[n/M] & n=0,\pm M, \pm 2M, \ldots \\ 0 & \mathrm{otherwise} \end{cases} \tag{10}\]
The upsampler inserts \(M-1\) zero samples between input samples. Another way to write down this relationship mathematically is
\[ w[n] = \sum_{k=-\infty}^{\infty} x[kM]\ \delta[n-kM] \]
which is analogous to CT pulse train conversation.
Figure 8 shows a block diagram for the upsampling operation and illustrates a good way to think about its operation relative to CT pulse train conversion.
In particular, we imagine a DT sequence signal \(x[n]\) that we plan to convert back into continuous time with sample spacing \(MT_s\). Intead of getting there direclty with an impulse train with spacing \(MT_s\), we can achieve the same result by first upsampling and then converting to an impulse train with spacing \(T_s\). Equating the two resulting spectra \(W_p(f)\), we can derive the input-output relationship for upsampling to be
\[ W(e^{j2\pi u})=X(e^{j2\pi u M}) \tag{11}\]
Since upsampling scales the time axis by the factor \(M\), it scales the frequency axis by the factor \(1/M\). This is yet another instance in which we have to be conscious about the periodic nature of the DT frequency axis and the DT Fourier transform. This is illustrated in Figure 9 for \(M=2\).
Imagining that \(x[n]\) represents samples of a CT signal \(x(t)\), this example illustrates that upsampling \(x[n]\) by \(2\) is not the same as sampling \(x(t)\) at twice the frequency. There are extra frequency components in \(X(e^{j2\pi u 2})\), for example, those that are centered at \(\pm 1/2\), that are called images.
Interpolation Filtering
To eliminate the images introduced by upsampling, we again apply a lowpass filter, however, this time in discrete time. For interpolation by a factor of \(M\), the convention is to apply an ideal DT lowpass filter with cutoff frequency \(1/(2 M)\). The frequency response of this filter is \[ H(e^{j2\pi u})= M \mathrm{rect}(u M) \tag{12}\]
for \(u\in[-1/2,1/2]\) and periodic with period \(1\), and its impulse response is
\[ h[n] = \mathrm{sinc}(n /M). \tag{13}\]
DT Sampling Theorem
A DT sequence \(x[n]\) is bandlimited to normalized frequency \(u_m > 0\) if \(X(e^{j2\pi u})=0\) for \(u_m < |u| < 1/2\). Note that this definition is slightly different that in CT, again because the DT frequency axis is periodic with period \(1\).
Theorem 2 (Reconstruction) If a DT sequence \(x[n]\) is bandlimited to maximum frequency \(0 < u_m 1/2\) and downsampled by the integer factor \(M\) such that \(M < 1/(2u_m)\), then it can be reconstructed through bandlimited interpolation of its \(M\)-spaced samples as
\[ \lim_{K\rightarrow \infty} \sum_{k=-N}^{N} x[kM]\ \mathrm{sinc}\left( \frac{n-kM}{M} \right)=x[n] \]
General Sample Rate Conversion
In many cases, we want to convert to a fractional sampling rate. If the sample rate conversion that we want to achieve is a factor \(M / N\), the usual approach is to first interpolate by the factor \(M\) and then decimate by the factor \(N\).
Practical Interpolation Filters Require Oversampling
The final, and practically most important, topic for this lecture addresses practical implementation of interpolation filters. The ideal lowpass filters discussed so far are non-causal and infinite duration, making them impossible to implement in practice.
In this section we consider more general filters for the CT interpolatio \[ \tilde{s}(t) = \mathrm{l.i.m.}_{N\rightarrow \infty} \sum_{n=-N}^{N} s(nT_s)\ h(t-nT_s). \tag{14}\]
Two classic examples are the zero-order hold, for which \(h(t)=\mathrm{rect}(t/T_s)\), and linear interpolation, which is developed in the problem set.
Another possibility is to time-limit the \(\mathrm{sinc}\) pulse by setting
\[ h(t)=\mathrm{rect}(t/(2 L T_s))\ \mathrm{sinc}(t/T_s), \]
which maintains \(2L\) sampling periods about \(t=0\). This multiplication in the time domain corresponds to convolution in the frequency domain. The smaller we seek to make the integer parameter \(L\), and thereby reduce the filtering delay of the interpolation, the wider the dispersion of the ideal lowpass response in the frequency domain.
For any one of these non-ideal filters, let \(f_1\) denote the largest frequency in the passband, and let \(f_2\) denote the smallest frequency in the stopband. The band of frequencies \([f_1,f_2]\) is often call the transition band of the lowpass filter.
Returning to our first illustration for a signal \(s(t)\) bandlimited to maximum frequency \(f_m\), to ensure that the spectrum replica \(S(f)\) is passed by the filter, we require \(f_1 > f_s\). However, to ensure that the replica \(S(f-f_s)\) is not passed by the filter, we require \(f_2 < f_s - f_m\). If \(f_m\) is fixed, and we want to widen the transition band by increasing \(f_2\), then the only way to ensure good reconstruction, i.e., \(\tilde{s}(t) \approx s(t)\), is to increase the sampling rate \(f_s\) as we increase \(f_2\). This often translates into \(f_s\) being considerably larger than the Nyquist rate \(2 f_m\), and is in general called oversampling.
Quantizer Parameters and Input-Output Relationships
A quantizer is a function that takes a continuous-valued input and products a discrete-valued output. It is generally a non-invertible function.
Figure 10 illustrates a fairly general quantizer with \(N\in\mathbb{Z}^+\) quantization levels. Often \(N\) is chosen to be a power of \(2\).
General Parameterization
Mathematically, we define a set of \(N+1\) values \(b_0 = -\infty < b_1 < b_2 < \ldots < b_{N-1} < b_N = \infty\) along the input axis, and we refer to them as the quantization region boundaries. These values partition the input space into the intervals
\[ \mathcal{R}_n = (b_{n-1},b_n] \]
for \(n=1,2,\ldots,N\), which are called the quantization regions.
We also define the set of \(N\) values \(a_1 < a_2 < \ldots a_N\) along the output axis, and refer to them as the quantizer representation points.
The functional definition of the quantizer is then \(q(y) = a_n\) if \(y\in \mathcal{R}_n\), \(n=1,2,\ldots,N\).
The minimum representation point \(a_1\) and the maximum representation point \(a_N\) define the output range \(a_N-a_1\) of the quantizer.
Minimum Mean-Square Error Quantization
If we are given a probability distribution on the input \(f_Y(y)\), then we can consider choosing the quantization region boundaries and representation points to minimze \(\mathbb{E}[|q(Y)-Y|^2]\). This is called minimum mean-square quantization.
Although a solution to this problem does not exist in general, there is an iterative approach that will converge to at least a local minimum that is called the Lloyd-Max Algorithm. The iterative algorithm alternates between improving the quantization region boundaries and the representation points. Specifically:
- Given the quantization representation points \(\{a_1,a_2,\ldots,a_N\}\), we select the quantization region boundaries according to correspond to the midpoints between representation points, i.e., \[ b_n = \frac{a_1+a_n}{2} \]
for \(n=1,2,\ldots,N-1\).
- Given the quantization region boundaries \(\{b_0,b_1,\ldots,b_N \}\), we select the quantization representation points to be the conditional mean of the random variable over the quantization region, i.e.,
\[ a_n = \mathbb{E}[Y | Y \in \mathcal{R}_n] \]
for \(n=1,2,\ldots,N\).
Uniform Quantization
For many practical implementations, the quantization regions \(\mathcal{R}_2, \mathcal{R}_3,\ldots,\mathcal{R}_{N-1}\) are set to be the same length \(\Delta > 0\), so that the quantization region boundaries satisfy
\[ b_n = b_{k-1} + \Delta \tag{15}\]
for \(k=2,3,\ldots,N_1\).
Furthermore, the representation points of the internal quantization regions are the midpoints of these intervals. Specifically,
\[ a_n = \begin{cases} \frac{b_n + b_{n-1}}{2} = b_{n-1} + \frac{\Delta}{2} & n=2,\ldots,N-2 \\ b_n - \frac{\Delta}{2} & n= 1 \\ b_n + \frac{\Delta}{2} & n=N-1 \end{cases} \tag{16}\]
The output range of such a quantizer is then \(a_N - a_1 = (N-1)\Delta\).
We can define a standard, \(K\)-bit uniform quantizer with output range \(1\) as the function
\[ q_K(y) = \begin{cases} \lfloor r \left( y + \frac{1}{2} + 1/(2r) \right) \rfloor / r - \frac{1}{2} & -\frac{1}{2} < y < +\frac{1}{2} \\ +\frac{1}{2} & y \ge +\frac{1}{2}\\ -\frac{1}{2} & y \le -\frac{1}{2} \end{cases} \tag{17}\]
where \(r=2^K-1\). For convenience, the quantizer representation points can be easily represented the integers \(0,1,\ldots,2^K-1\) or their binary representations.
This quantizer limits the input to \(\pm 1/2\). If we want to limit the input and output to \(\pm A\), then we can scale the quantizer by applying \(2A q_K(y / (2A))\).
The input-output relationships for this standardized and scaled standardized uniform quantizer as illustrated in the Python code below.
def q(y,K):
z = np.zeros(np.size(y))
for i in range(np.size(z)):
if y[i] >= 0.5:
z[i] = 0.5
elif y[i] <= -0.5:
z[i] = -0.5
else:
r = np.power(2.0,K)-1.0
z[i] = np.floor( r*(y[i]+0.5 + 1.0/r/2.0 ) ) / r - 0.5
return zK = 4 # 4 bits, 16 levels
M = 100*np.power(2.0,K) # 100 samples per quantization region
A = 1 # Expand range of quantizer to \pm Ay_values = np.arange(-2*A,+2*A,A/M)
qKy = 2.0*A*q(y_values/(2.0*A),K)plt.plot(y_values,qKy,'-')
plt.plot(y_values,y_values,'-.')
plt.xlabel("$y$")
plt.ylabel("$2 A q_K(y/(2A))$")
plt.show()
Effect of Uniform Quantization on Sinusoidal Signals
u0 = 1/4
n = np.arange(0,128.0/u0,1)
y = np.cos(2*np.pi*u0*n)K=8
A=2.5
yq=2.0*A*q(y/(2.0*A),K)nfft=512
u=np.fft.fftshift(np.fft.fftfreq(nfft))
Y=np.fft.fftshift(np.fft.fft(y,nfft))
Yq=np.fft.fftshift(np.fft.fft(yq,nfft))
plt.plot(u,10.0*np.log10(np.abs(Y)**2),'-')
plt.plot(u,10.0*np.log10(np.abs(Yq)**2),'--')
plt.xlabel("$u$")
plt.ylabel("$Y(e^{j2\pi u})$ (dB)")
plt.axis([-1/2.0, 1/2.0, -80.0, 60.0])
plt.show()
Take Aways
Quantization converts a continuous-amplitude input into a discrete-valued output, and it inherently introduces some distortion.
It is important to consider the output range of the quantizer relative to our input signal amplitude to ensure that we obtain a low error representation
- If the input signal amplitude is too small, quantization noise will dominate the output.
- If the input signal amplitude is too large, clipping distortion will dominate the output.