Now consider the following multirate processing scheme in which $L(z)$ is an ideal \emph{lowpass} filter with cutoff frequency $\pi/2$ and $H(z)$ is an ideal \emph{highpass} filter with cutoff frequency $\pi/2$:
+ Consider the input-output characteristic of the following multirate systems. Remember that, technically, one cannot talk of transfer functions in the case of multirate systems since sampling rate changes are not time-invariant. It may happen, though, that by carefully designing the processing chain, the input-output characteristic does indeed implement a time-invariant transfer function.
+ \begin{enumerate}
+ \item Find the overall transformation operated by the following system:
+ \item Assume $H(z) = A(z^2)+z^{-1}B(z^2)$ for arbitrary $A(z)$ and $B(z)$. Show that the transfer function of the following system is equal to $A(z)$.
- \item Assume $H(z) = A(z^2)+z^{-1}B(z^2)$ for arbitrary $A(z)$ and $B(z)$. Show that the transfer function of the following system is equal to $A(z)$:
The standard model for the error introduced by quantization is that of additive white noise, where the noise is i.i.d. and uncorrelated with the signal:
\[
\mathcal{Q}\{x[n]\} = x[n] + e[n].
\]
If the quantizer is uniform and the input signal is also uniformly distributed, the probability distribution of the noise samples is also uniform over the interval $[-\Delta/2, \Delta/2]$ where $\Delta$ is the quantization step size.
The quantized signal $\hat{x}[n]$ is processed by a digital filter with impulse response
\[
h[n]=\frac{1}{2}[a^n+(-a)^n]u[n].
\]
Determine the variance of the noise at the output of the filter and determine the signal-to-noise ratio at the output assuming that the original signal is also white and i.i.d.
\begin{exercise}{Digital processing of continuous-time signals}
-
- In your grandmother's attic you just found a treasure: a collection of super-rare 78rpm vinyl jazz records. The first thing you want to do is to transfer the recordings to compact discs, so you can listen to them without wearing out the originals. Your idea is obviously to play the record on a turntable and use an A/D converter to convert the line-out signal into a discrete-time sequence, which you can then burn onto a CD. The problem is, you only have a ``modern'' turntable, which plays records at 33rpm. Since you're a DSP wizard, you know you can just go ahead, play the 78rpm record at 33rpm and sample the output of the turntable at 44.1~KHz. You can then manipulate the signal in the discrete-time domain so that, when the signal is recorded on a CD and played back, it will sound right.
-
- Design a system which performs the above conversion. If you need to get on the right track, consider the following:
+ In your grandmother's attic you just found a treasure: a collection of super-rare $78$ rpm vinyl\index{vinyl} jazz records. The first thing you want to do is to transfer the recordings to digital format, so you can listen to them without wearing out the originals. Your idea is obviously to play the record on a turntable and use an A/D converter to convert the line-out signal into a discrete-time sequence, which you can then store or burn onto a CD. The problem is, you only have a ``modern'' turntable, which plays records at $33$ rpm. Since you're a DSP wizard, you know you can just go ahead, play the $78$ rpm record at $33$~rpm and sample the output of the turntable at $44.1$~KHz. You can then manipulate the signal in the discrete-time domain so that, when the signal is recorded on a CD and played back, it will sound right.
+
+ In order to design a system that performs the above conversion consider the following:
\begin{itemize}
- \item Call $s(t)$ the continuous-time signal encoded on the 78rpm vinyl (the jazz music) \item Call $x(t)$ the continuous-time signal you obtain when you play the record at 33rpm on the modern turntable
- \item Let $x[n] = x(nT_s)$, with $T_s = 1/44100$.
+ \item Call $s(t)$ the continuous-time signal encoded on the $78$ rpm vinyl (the jazz music).
+ \item Call $x(t)$ the continuous-time signal you obtain when you play the record at $33$ rpm on the modern turntable.
+ \item Let $x[n] = x(nT_s)$, with $T_s = 1/44,100$.
\end{itemize}
- and answer the following questions:
+ Answer the following questions:
\begin{enumerate}
- \item Express $x(t)$ in terms of $s(t)$. \item Sketch the Fourier transform $Xf)$ when $S(f)$ is as in the following figure. The highest nonzero frequency of $S(f)$ is $F_{\max} = 16000$~Hz (old records have a smaller bandwidth than modern ones).
+ \item Express $x(t)$ in terms of $s(t)$.
+ \item Sketch the Fourier transform $X(f)$ when $S(f)$ is like so:
- \item Design a system to convert $x[n]$ into a sequence $y[n]$ so that, when you interpolate $y[n]$ to a continuous-time signal $y(t)$ with interpolation period $T_s$, you obtain $Y(f) = S(f)$.
- \item What if you had a turntable which plays records at 45rpm? Would your system be different? Would it be better?
+ The highest nonzero frequency of $S(f)$ is $F_{\max} = 16,000$~Hz (old records have a smaller bandwidth than modern ones).
+ \item Design a system to convert $x[n]$ into a sequence $y[n]$ so that, when you interpolate $y[n]$ to a continuous-time signal $y(t)$ with interpolation period $T_s$, you obtain $Y(f) = S(f)$.
+ \item What if you had a turntable which plays records at $45$ rpm? Would your system be different? Would it be better?
A single sample of $x[n]$ may have been corrupted and we would like to approximately or exactly recover it. We denote $n_0$ the time index of the corrupted sample and $\hat{x}[n]$ the corresponding corrupted sequence.
\begin{enumerate}
\item Specify a strategy for exactly recovering $x[n]$ from $\hat{x}[n]$ if $n_0$ is known.
\item What would you do if the value of $n_0$ is not known?
\item Now suppose we have $k$ corrupted samples: what is the condition that $X(e^{j \omega })$ must satisfy to be able to exactly recover $x[n]$?
+ We will operate in the $z$-domain; recall that, given the $z$-transform of the original sequence $X(z)$, the $z$-transform of a 2-upsampled sequence is $X(z^2)$ whereas the $z$-transform of a 2-downsampled sequence is $[X(z^{1/2}) + X(-z^{1/2})]/2$. With this:
\begin{enumerate}
- \item In the $z$-domain, let us denote the downsampling by 2 and upsampling by 2 operations by $D_2\{\cdot\}$ and $U_2\{\cdot\}$ respectively, so that
- \begin{align*}
- U_2\{X(z)\} &= X(z^2) \\
- D_2\{X(z)\} &= [X(z^{1/2}) + X(-z^{1/2})]/2
- \end{align*}
- \begin{itemize}
- \item Using $z$-transforms, downsampling by 2 followed by filtering by $H(z)$ can be written as
+ The two operations are thus equivalent. (Please note in $H(z^2)$ how it's the \textit{argument} $z$ of the $z$-transform that is replaced by $\pm z^{1/2}$, so that the minus sign goes away. This makes sense: if a sequence has been upsampled, downsampling simply returns it to its original form.)
+ \item Filtering by $H(z)$ followed by upsampling by 2 can be written as
+ \[
+ Y_1(z)= H(z^2)X(z^2).
+ \]
Upsampling by 2 followed by filtering by $H(z^2)$ can be written as
- \begin{align*}
- Y(z) &= H(z^2)U_2\{X(z)\} \\
- &= H(z^2)X(z^2).
- \end{align*}
- The two operations are thus equivalent.%
- \end{itemize}
+ \[
+ Y_2(z) = H(z^2)X(z^2).
+ \]
+ The two operations are thus equivalent.
+ \end{enumerate}
+\end{solution}
- \item Using the identities proven in (a), the system can be redrawn as
- The lower branch contains an upsampler followed by a delay and a downsampler. The output of such a system is easily seen to be $0$. Thus only the upper branch remains and the final transfer function of the system is given by
+ The lower branch contains an upsampler by two followed by a delay and a downsampler by two. The output of such a processing chain is easily seen to be equal to $0$ for all $n$. Thus only the upper branch remains and the final transfer function of the system is given by
\begin{eqnarray*}
- \frac{Y(z)}{X(z)}=A(z).
+ H(z) = A(z).
\end{eqnarray*}
- \item System 1 is described by the following equation %
+ \item In both systems, call $w[n]$ the signal just before the final downsampler by two; the $z$-transform of the output will therefore be
+ \[
+ Y(z) = [W(z^{1/2}) + W(-z^{1/2})]/2.
+ \]
+ For the first system we have that $W(z) = H(z)G(z)X(z^2)$ so
- \item We need to change the sampling rate so that, when $y[n]$ is interpolated at 44.1~KHz its spectrum is equal to $S(f)$. The rational sampling rate change factor is clearly $33/78$ which is simply $11/26$ after factoring. The processing scheme is as follows:
- where $L(z)$ is a lowpass filter with cutoff frequency $\pi/26$ and gain $L_0 = 11$; both the sampler and interpolator work at $F_s = 1/T_s = 44100$~Hz. We have:
+ \item We need to change the sampling rate so that, when $y[n]$ is interpolated at 44.1~KHz its spectrum is equal to $S(f)$. The rational sampling rate change factor is clearly $33/78$ which is simply $11/26$ after factoring. The processing scheme is as follows:
+ where $L(z)$ is a lowpass filter with cutoff frequency $\pi/26$ and gain $L_0 = 11$; both the sampler and interpolator work at $F_s = 1/T_s = 44100$~Hz. We have:
- \item The sampling rate change scheme stays the same except that now $45/78 = 15/26$. Therefore, the final upsampler has to compute more samples than in the previous scheme. The computational load of the sampling rate change is entirely dependent on the filter $L(z)$. If we upsample more before the output, we need to compute more filtered samples and therefore at 45rpm the scheme is less efficient.
-\end{enumerate}
-
+ \item The sampling rate change scheme stays the same except that now $45/78 = 15/26$. Therefore, the final upsampler has to compute more samples than in the previous scheme. The computational load of the sampling rate change is entirely dependent on the filter $L(z)$. If we upsample more before the output, we need to compute more filtered samples and therefore at 45rpm the scheme is less efficient.
Given that $X(e^{j \omega })=0$ for $\frac{\pi}{3} \leq | \omega | \leq \pi$, $x[n]$ can be thought of as a signal that has been sampled at $3$ times the Nyquist frequency. Therefore, we can downsample the signal without losing information at least by a factor of 3.
\begin{enumerate}
\item Assume $n_0$ is odd; we can then downsample $x[n]$ by 2, without loss of information and the corrupted sample will be discarded in the downsampling operation. We can then upsample by 2 and recover the original signal eliminating the error. If $n_0$ is even, simply shift the signal by 1 and perform the same operation.
\item If the value of $n_0$ is not known, we need to determine whether $n_0$ is odd or even. We can write
since, by hypothesis, $X(e^{j \frac{\pi}{2} })=0$. Therefore, If $\hat{X}(e^{j \frac{\pi}{2} })$ is real, $n_0$ is even and if it is imaginary, $n_0$ is odd.
\item If there are $k$ corrupted samples, the worst case is when the corrupted samples are consecutive. In that case we need to downsample $\hat{x}[n]$ by a factor of $k$ and then upsample it back. To do so without loss of information it must be:
i.e. a signal obtained from the discrete-time sequence using a zero-centered zero-order hold with interpolation period $T_s = 1$s. Let $X_0(f)$ be the Fourier transform of $x_0(t)$.
\begin{enumerate}
\item Express $X_0(f)$ in terms of $X(e^{j\omega})$.
\item Compare $X_0(f)$ to $X(f)$, where $X(f)$ is the spectrum of the continuous-time signal obtained using an ideal sinc interpolator with $T_s=1$:
Comment on the result: you should point out two major problems.
\item The signal $x(t)$ can be obtained back from the zero-order hold interpolation via a continuous-time filtering operation:
\[
x(t)=x_0(t)\ast g(t).
\]
Sketch the frequency response of the filter $g(t)$.
\item Propose two solutions (one in the continuous-time domain, and another in the discrete-time domain) to eliminate or attenuate the distortion due to the zero-order hold. Discuss the advantages and disadvantages of each.
\end{enumerate}
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
\begin{enumerate}
\item
\begin{align*}
X_0(f) &= \int_{-\infty}^{\infty} x_0(t)e^{-j2\pi f t} dt\\
&= \int_{-\infty}^{\infty} \sum_{n=-\infty}^{\infty}x[n]\mbox{rect}\,(t-n)e^{-j2 \pi f t} dt \\
&= \sum_{n=-\infty}^{\infty} x[n] \int_{-\infty}^{\infty} \mbox{rect}\,(t-n)e^{-j2\pi f t} dt \\
&= \sum_{n=-\infty}^{\infty}x[n] e^{-j2\pi f n} \, \int_{-1/2}^{1/2} e^{-j2\pi f \tau} d\tau \\
\item To understand the effects of the zero-order hold consider for instance a discrete-time signal with a triangular
spectrum, as shown in the left panel below. We know that sinc interpolation will exactly preserve the shape of the spectrum and return a continuous-time signal that is strictly bandlimited to the $[-F_s/2, F_s/2]$ interval (with $F_s = 1/T_s = 1$), that is:
Conversely, the spectrum of the continuous-time signal interpolated by the zero-order hold is the product of $X(e^{j2\pi f})$, which is periodic with period $F_s = 1$ Hz, and the term $\sinc(f)$, whose first spectral null is for $f=1$ Hz. Here are the two terms, and their product, in magnitude:
\dspFunc{x \dspPeriodize \dspTri{0}{1} x \dspSinc{0}{2} mul abs}
\dspCustomTicks[axis=x]{0 0 2 1 4 2}
\end{dspPlot}
\end{center}
There are two main problems in the zero-order hold interpolation as compared to the sinc interpolation:
\begin{itemize}
\item The zero-order hold interpolation is NOT bandlimited: the $2\pi$-periodic replicas of the digital spectrum leak through in
the continuous-time signal as high frequency components. This is due to the sidelobes of the interpolation function in the frequency domain (rect in time $\leftrightarrow$ sinc in frequency) and it represents an undesirable high-frequency content which is typical of all local interpolation schemes.
\item There is a distortion in the main portion of the spectrum (that between$-F_s/2$ and $F_s/2 = 0.5$ Hz) due to the non-flat frequency response of the interpolation function. It can be seen if we zoom in the baseband portion:
\item A first solution is to compensate for the distortion introduced by $G(f)$ in the discrete-time domain. This is
equivalent to pre-filtering $x[n]$ with a discrete-time filter of magnitude $1/G(e^{j2\pi f})$ . The advantages of this method is that digital filters such as this one are relatively easy to design and that the filtering can be done in the discrete-time domain. The disadvantage is that this approach does not eliminate or attenuate the high frequency leakage outside the baseband.
In continuous time, one could cascade the interpolator with an analog lowpass filter to eliminate the leakage. The disadvantage is that it is hard to implement an analog lowpass which can also compensate for the in-band distortion introduced by $G(f)$; such a filter will also introduce unavoidable phase distortion (no analog filter has linear phase).
Consider a local interpolation scheme as in the previous exercise but now the characteristic of the interpolator is:
\[
i(t) = \begin{cases}
1 -2|t| & |t| \leq 1/2 \\
0 & \mbox{otherwise}
\end{cases}
\]
This is a triangular characteristic but the same unit support as the zero-order hold. If we pick an interpolation interval $T_s$ and interpolate a given discrete-time signal $x[n]$ with $I(t)$, we obtain the continuous-time signal
Assume that the spectrum of $x[n]$ between $-\pi$ and $\pi$ is
\[
X(e^{j\omega}) = \begin{cases}
1 & |\omega| \leq 2 \pi/3 \\
0 & \mbox{otherwise}
\end{cases}
\]
(with the obvious $2\pi$-periodicity over the entire frequency axis).
\begin{enumerate}
\item Compute and sketch the Fourier transform $I(f)$ of the interpolating function $i(t)$. Recall that the triangular function can be expresses as the convolution of a suitably scaled rect with itself.
\item Sketch the Fourier transform $X(f)$ of the interpolated signal $x(t)$; in particular, clearly mark the Nyquist frequency $F_s/2$
\item The use of $i(t)$ instead of a sinc interpolator introduces two types of errors: briefly describe them.
\item To eliminate the error in the baseband $[-F_s/2,F_s/2]$ we can pre-filter the signal $x[n]$ before interpolating with $i(t)$. Write the frequency response of the required discrete-time pre-filter $H(e^{j\omega})$.
\end{enumerate}
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
\begin{enumerate}
\item We know that the convolution of $\rect(t)$ with itself produces a rectangular function between $-1$ and $1$ and with height 1. To obtain a rectangular function of support 1 and height 1, we need to shrink the rects by a factor of two and compensate for the value in $t=0$ so that
\item like before, here are two types of error, in-band and out-of-band:
\begin{itemize}
\item $\mathbf{In-band}$: the spectrum between $[-F_s/2, F_s/2]$ (the baseband) is distorted by the non-flat response of the interpolating function over the baseband.
\item $\mathbf{Out-of-band}$: the periodic copies of $X(e^{j2\pi f T_s})$ outside of $[-F_s/2, F_s/2]$ are not eliminated by the interpolation filter, since it is not an ideal lowpass.
\end{itemize}
\item As before, e need to undo the linear distortion introduced by the nonflat response of the interpolation filter in the baseband. The idea is to have a modified spectrum $H(e^{j\omega})X(e^{j\omega})$ so that, over $[-F_s/2, F_s/2]$, it is
\[
X(f) = X(e^{j2\pi f T_s}).
\]
If we use $H(e^{j\omega})X(e^{j\omega})$ in the interpolation formula, we have
\[
X(f) = \frac{T_s}{2} \, H(e^{j2\pi f T_s})X(e^{j2\pi f T_s})\, \sinc^2\left(\frac{f T_s}{2}\right)
\]
so that
\[
H(e^{j2 \pi f T_s}) = [\frac{T_s}{2}\sinc^2\left(\frac{f T_s}{2}\right)]^{-1}.
\]
Therefore, the frequency response of the digital filter will be
An alternative way of describing the sampling operation relies on the concept of \textit{modulation by a pulse train}. Given a sampling interval $T_s$, a continuous-time pulse train $p(t)$ is an infinite collection of equally spaced Dirac deltas:
\[
p(t) = \sum_{k=-\infty}^{\infty}\delta(t-kT_s).
\]
The pulse train is the used to modulate a continuous-time signal:
\[
x_s(t) = p(t)\,x(t).
\]
Intuitively, $x_s(t)$ represents a ``hybrid'' signal where the nonzero values are only those of the discrete time samples that would be obtained by raw-sampling $x(t)$ with period $T_s$; however, instead of representing the samples a countable sequence (i.e. with a different mathematical object) we are still using a continuous-time signal that is nonzero only over infinitesimally short instants centered on the sampling times. Using Dirac deltas allows us to embed the instantaneous sampling values in the signal.
Note that the Fourier Transform of the pulse train is
(where, as per usual, $F_s = 1/T_s$). This result is a bit tricky to show, but the intuition is that a periodic set of pulses in time produces a periodic set of pulses in frequency and that the spacing between pulses frequency is inversely proportional to the spacing between pulses in time.
Derive the Fourier transform of $x_s(t)$ and show that if $x(t)$ is bandlimited to $F_s/2$, where $F_s = 1/T_s$, then we can reconstruct $x(t)$ from $x_s(t)$.
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
By using the modulation theorem, the product in time becomes a convolution in frequency:
In other words, the spectrum of the delta-modulated signal is the periodization (with period $F_s=1/T_s$) of the original spectrum. If the latter is bandlimited to $F_s/2$ there will be no overlap between copies in the periodization and therefore $x(t)$ can be obtained simply by lowpass filtering $x_s(t)$ in the continuous-time domain.
\item What is the bandwidth of the signal? What is the minimum sampling frequency that satisfies the sampling theorem?
\item If we sample the signal with a sampling frequency $F_a = 2f_0$, clearly there will be aliasing. Plot the DTFT of the resulting discrete-time signal $x_a[n] = x_c(n/F_a)$.
\item Suggest a way to perfectly reconstruct $x_c(t)$ from $x_a[n]$.
\item From the previous example it would appear that we can exploit ``gaps'' in the original spectrum to reduce the sampling frequency necessary to losslessly sample a bandpass signal. In general, what is the minimum sampling frequency that we can use to sample with no loss a real-valued signal whose frequency support on the positive axis is $[f_0, f_1]$ (with the usual symmetry around zero, of course)?
\end{enumerate}
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
\begin{enumerate}
\item The highest nonzero frequency is $2f_0$ and therefore $x_c(t)$ is $4f_0$-bandlimited. The minimum sampling frequency that satisfies the sampling theorem is $F_s = 4f_0$. Note however that the support over which the (positive) spectrum is nonzero is the interval $[f_0, 2f_0]$ so that one could say that the total \emph{effective} bandwidth of the signal is only $2f_0$.
\item The digital spectrum will be the periodized continuous-time spectrum, rescaled to $[-\pi, \pi]$; the periodization after sampling at a frequency $F_a = 2f_0$, yields
for instance, as $\omega$ goes from zero to $\pi$, the nonzero contribution to the DTFT will be the term $X_c(\frac{\omega}{\pi}f_0 - 2f_0)$ where the argument goes from $-2f_0$ to $-f_0$. The spectrum is represented here with periodicity explicit:
\item filter in continuous time $x_p(t)$ with an ideal bandpass filter with (positive) passband equal to $[f_0, 2f_0]$ to obtain $x_c(t)$.
\end{itemize}
\item The effective \emph{positive} bandwidth of a signal whose spectrum is nonzero only over $[-f_1, -f_0] \cup [f_0, f_1]$ is $W = f_1 - f_0$. Obviously the sampling frequency must be at least equal to the \textit{total} effective bandwidth, so a first condition on the sampling frequency is $F_s \geq 2W.$ We can now distinguish two cases.
\begin{itemize}
\item[1)] assume $f_1$ is a multiple of the positive bandwidth, i.e.\ $f_1 = MW$ for some integer $M > 1$ (for $x[n]$ above, it was $M = 2$). Then the argument we made before can be easily generalized: if we pick $F_s = 2W$ and sample we have that
Since $f_0 = f_1 - W = (M-1)W$, this translates to
\begin{align*}
(2k+M-1)W & \leq f \leq (2k+M)W \\
(2k-M)W & \leq f \leq (2k-M+1)W
\end{align*}
which, once again, are non-overlapping intervals.
\item[2)] if $f_1$ is \emph{not} a multiple of $W$ the easiest thing to do is to decrease the lower frequency $f_0$ to a new {\em smaller} frequency $f_0'$ so that the new positive bandwidth $W' = f_1 - f_0'$ divides $f_1$ exactly. In other words we set a new lower frequency $f_0'$ so that it will be $f_1 = M(f_1-f_0')$ for some integer $M>1$; it is easy to see that
\[
M = \biggl\lfloor \frac{f_1}{f_1 - f_0} \biggr\rfloor.
\]
since this is the maximum number of copies of the $W$-wide spectrum which fit \emph{with no overlap} in the $[0, f_0]$ interval. If $W > f_0$ obviously we cannot hope to reduce the sampling frequency and we have to use normal sampling. This artificial change of frequency will leave a small empty ``gap'' in the new bandwidth $[f_0', f_1]$, but that's no problem. Now we use the previous result and sample at $F_s = 2(f_1 - f_0')$ with no overlaps. Since $f_1 - f_0' = f_1/M$, we have that, in conclusion, the minimum sampling frequency is
- \item Assume $H(z) = A(z^2)+z^{-1}B(z^2)$ for arbitrary $A(z)$ and $B(z)$. Show that the transfer function of the following system is equal to $A(z)$:
+ \item Downsampling by $2$ followed by filtering by $H(z)$ is equivalent to filtering by $H(z^2)$ followed by downsampling by $2$.
+ \item Filtering by $H(z)$ followed by upsampling by $2$ is equivalent to upsampling by $2$ followed by filtering by $H(z^2)$.
+ \end{enumerate}
\end{exercise}
}\fi
\ifanswers{%
-
\begin{solution}{}
+ We will operate in the $z$-domain; recall that, given the $z$-transform of the original sequence $X(z)$, the $z$-transform of a 2-upsampled sequence is $X(z^2)$ whereas the $z$-transform of a 2-downsampled sequence is $[X(z^{1/2}) + X(-z^{1/2})]/2$. With this:
\begin{enumerate}
- \item In the $z$-domain, let us denote the downsampling by 2 and upsampling by 2 operations by $D_2\{\cdot\}$ and $U_2\{\cdot\}$ respectively, so that
- \begin{align*}
- U_2\{X(z)\} &= X(z^2) \\
- D_2\{X(z)\} &= [X(z^{1/2}) + X(-z^{1/2})]/2
- \end{align*}
- \begin{itemize}
- \item Using $z$-transforms, downsampling by 2 followed by filtering by $H(z)$ can be written as
- The two operations are thus equivalent. (Please note in $H(z^2)$ how it's the \textit{argument} $z$ of the $z$-transform that is replaced by $\pm z^{1/2}$, so the minus sign goes away. This makes sense: if a sequence has been upsampled, downsampling simply returns it to its original form.)
- \item Filtering by $H(z)$ followed by upsampling by 2 can be written as
+ \item Downsampling by 2 followed by filtering by $H(z)$ yields
+ The two operations are thus equivalent. (Please note in $H(z^2)$ how it's the \textit{argument} $z$ of the $z$-transform that is replaced by $\pm z^{1/2}$, so that the minus sign goes away. This makes sense: if a sequence has been upsampled, downsampling simply returns it to its original form.)
+ \item Filtering by $H(z)$ followed by upsampling by 2 can be written as
+ \[
+ Y_1(z)= H(z^2)X(z^2).
+ \]
Upsampling by 2 followed by filtering by $H(z^2)$ can be written as
+ Consider the input-output characteristic of the following multirate systems. Remember that, technically, one cannot talk of transfer functions in the case of multirate systems since sampling rate changes are not time-invariant. It may happen, though, that by carefully designing the processing chain, the input-output characteristic does indeed implement a time-invariant transfer function.
+ \begin{enumerate}
+ \item Find the overall transformation operated by the following system:
+ \item Assume $H(z) = A(z^2)+z^{-1}B(z^2)$ for arbitrary $A(z)$ and $B(z)$. Show that the transfer function of the following system is equal to $A(z)$.
- The lower branch contains an upsampler followed by a delay and a downsampler. The output of such a system is easily seen to be $0$. Thus only the upper branch remains and the final transfer function of the system is given by
+ The lower branch contains an upsampler by two followed by a delay and a downsampler by two. The output of such a processing chain is easily seen to be equal to $0$ for all $n$. Thus only the upper branch remains and the final transfer function of the system is given by
\begin{eqnarray*}
- \frac{Y(z)}{X(z)}=A(z).
+ H(z) = A(z).
\end{eqnarray*}
- \item System 1 is described by the following equation %
+ \item In both systems, call $w[n]$ the signal just before the final downsampler by two; the $z$-transform of the output will therefore be
+ \[
+ Y(z) = [W(z^{1/2}) + W(-z^{1/2})]/2.
+ \]
+ For the first system we have that $W(z) = H(z)G(z)X(z^2)$ so
Compute the values of $y[n]$ for $0 \le n \le 6$, showing your calculation method.
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
Intuitively, the upsampling by four followed by $I_1(z)$ creates a linear interpolation over three ``extra'' samples between the original values the (``connect-the-dots'' strategy). The delay by two followed by the downsampler selects the midpoint of each interpolation interval. As a whole, the chain implements a fractional delay of half a sample using a linear interpolator so that
Since the interpolator $I(z)$ has finite support of length 7, we can concentrate on the interval $[-4, 4]$ and extend the result to the other points. Linear interpolation fills in the gaps while the delay shifts the interpolated signal by two towards the right:
Consider a stationary i.i.d. random process $x[n]$ whose samples are uniformly distributed over the $[-1, 1]$ interval. Consider a quantizer $\mathcal{Q}\{\cdot\}$ with the following characteristic:
where the integral is computed over the range of the input and $f_x(x)$ is the probability density function of the samples. We can split the integral over the non-overlapping quantization intervals $i_k$ as
Consider a stationary i.i.d. random process $x[n]$ whose samples are uniformly distributed over the $[-1, 2]$ interval. The process is uniformly quantized with a 1-bit quantizer with the following characteristic:
\[
\mathcal{Q}\{x\} = \begin{cases}
-1 & \text{if } x < 0 \\
1 & \text{if } x \geq 0
\end{cases}
\]
Compute the signal to noise ratio at the output of the quantizer
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
The input is uniformly distributed so its power will be
A {\em deadzone} quantizer is a quantizer that has a quantization interval centered around zero. To see the effects of the deadzone quantizer on SNR consider an i.i.d. discrete-time process $x[n]$ whose values are in the $[-1, 1]$ interval. Consider the following uniform 2-bit quantizers for the interval:
\multido{\n=0+1}{3}{\pscircle*[linecolor=gray](!0 \n\space 1 sub 2 mul 3 div){2pt}}
\psplot[plotpoints=800,linewidth=2pt,linecolor=gray]{-0.99}{0.99}{x 3 mul 2 div round 2 mul 3 div}
\end{pspicture}
\end{tabular}
\end{center}
Both quantizers operate at two bits per sample but the deadzone quantizer "wastes" a fraction of a bit since it has only 3 quantization intervals instead of 4; for a uniformly distributed input, therefore, the SNR of the deadzone quantizer is smaller than the SNR of the standard quantizer. Assume now that the probability distribution for each input sample is the following:
\[
P[x[n] = \alpha] = \begin{cases}
0 & \mbox{if $|\alpha| > 1$} \\
p & \mbox{if $|\alpha| = 0$} \\
(1-p)/2 & \mbox{otherwise}
\end{cases}
\]
In other words, each sample is either zero with probability $p$ or drawn from a uniform distribution over the $[-1, 1]$ interval; we can express this distribution as a pdf like so:
\[
f(x) = \frac{1-p}{2} + p\delta(x)
\]
Determine the minimum value of $p$ for which it is better to use the deadzone quantizer, i.e. the value of $p$ for which the SNR of the deadzone quantizer is larger that the SNR of the uniform quantizer.
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
First of all, since the input signal is the same, in order to compare SNRs we just need to compare the mean square errors of the two quantizers and find out the value of $p$ for which the deadzone quantizer's MSE is smaller than the normal quantizer's MSE. The formula for the MSE of a scalar quantizer over the $[-1, 1]$ interval (under the usual hypotheses of iid samples with pdf $f(x)$) is
The number of quantization levels in the two quantizers are $M=M_n=4$ for the normal 2-bit quantizer and $M=M_d=3$ for the deadzone quantizer. Let's compute the MSE for the normal quantizer using the composite pdf for the input
As we have seen, a continuous-time signal can be converted to a discrete-time sequence via sampling. By changing the value of the sampling rate we can obtain an arbitrary number of discrete-time signals from the same original continuous-time source; the number of samples per second will increase or decrease linearly with the sampling rate and, according to whether the conditions of the sampling theorem are satisfied or not, the resulting discrete-time sequences will be an exact representation of the original signal or will be affected by aliasing. Mutirate theory studies the relationship between such sequences; or, in other words, it addresses the question of whether we can transform a sampled representation into another with a different underlying sampling frequency purely from within discrete time.
The primary application of multirate theory is digital sampling rate conversion, an operation that becomes necessary when the original sampling frequency of a stored signal is incompatible with the working rate of the available D/A device. But we will see that multirate also plays an important role in a variety of applied signal processing algorithms; in particular, the technique of \textit{oversampling} is often used in situations where an increase of the data rate is used to improve the performance of the analog elements in the processing chain. Since speeding up digital computations is in general much cheaper than using high-performance analog circuits, oversampling is commonly used in consumer-level devices. Other applications of multirate range from efficient filter design, to spectral shaping for communication systems, to data compression standards. And, finally, from a more theoretical point of view, multirate is a fundamental ingredient of advanced analysis techniques which go under the name of time-frequency decompositions.
-Upsampling a sequence by an integer factor $N$ produces a higher-rate sequence by creating $N-1$ new samples for every sample in the original signal. In its basic form, an upsampler simply inserts $N-1$ zeros after every input sample, as shown in Figure~\ref{fig:mr:up}. If we denote by $\mathcal{U}_N$ the upsampling operator, we have
+Upsampling a signal by $N$ produces a higher-rate sequence by creating $N-1$ new samples for every original data point. In its basic form, an upsampler simply inserts $N-1$ zeros after every input sample, as shown in Figure~\ref{fig:mr:up}. If we denote by $\mathcal{U}_N$ the upsampling operator, we have
No information is lost via upsampling and the original signal can be easily recovered by discarding the extra samples (as we will see in the next section, upsampling by $N$ followed by downsampling by $N$ returns the original signal). The spectrum of an upsampled signal can be easily obtained by first considering its \ztrans\:
In the frequency domain, therefore, upsampling is simply a contraction of the frequency axis by a factor of $N$. The inherent $2\pi$-periodicity of the spectrum must be taken into account so that, in this contraction, the periodic repetitions of the base spectrum are ``pulled in'' the $[-\pi, \pi]$ interval. The effects of upsampling on a signal's spectrum are shown graphically for a simple signal in Figures~\ref{fig:mr:upA} and~\ref{fig:mr:upB}; in all figures the top panel shows the original spectrum $X(e^{j\omega})$ over $[-\pi, \pi]$; the middle panel shows the same spectrum over a wider range to make the $2\pi$-periodicity explicitly; the last panel shows the upsampled spectrum $X_{N\uparrow}(e^{j\omega})$, highlighting the rescaling of the $[-N\pi, N\pi]$ interval.\index{upsampling|)} As a rule of thumb, upsampling ``brings in'' exactly $N$ copies of the original spectrum over the $[-\pi, \pi]$ interval even if, in the case of an even upsampling factor, one copy is split between the negative and positive frequencies.
\caption{Upsampling by $4$ in the time domain: original signal (top panel); upsampled signal, where 3 zeros have been appended to each original sample (bottom panel). Note the difference in time indexes between top and bottom panels. }\label{fig:mr:up}
The upsampled signal in~(\ref{eq:mr:up}), with its $N-1$ zeros between original samples, exhibits two problems. In the time domain, the upsampled signal looks ``jumpy'' because of the periodic zero runs. This is the discrete-time equivalent to a lack of smoothness, since the signal keeps dropping back to zero, and it is apparent in the bottom panel of Figure~\ref{fig:mr:up}. In the frequency domain, simple upsampling has ``brought in'' copies of the original spectrum in the $[-\pi, \pi]$ interval, creating spurious high frequency content. These two issues are actually one and the same and they can be solved, as we will see, by using an appropriate filter.
The problem of filling the gaps between nonzero samples in an upsampled sequence is, in many ways, similar to the discrete- to continuous-time interpolation problem of Section~\ref{sec:is:interp}, except that now we are operating entirely in discrete time. If we adapt the interpolation schemes that we have already studied, we have the following cases\index{interpolation!in multirate}:
In this discrete-time interpolation scheme, also known as \emph{piecewise-constant interpolation}, after upsampling by $N$, we use a filter with impulse response
\begin{equation}\label{eq:mr:zoh}
h_0 [n] =
\begin{cases}
1 & \ n = 0,1, \ldots, N-1 \\
0 & \mbox{ otherwise}
\end{cases}
\end{equation}
which is shown in Figure~\ref{fig:mr:zoh}-(a). This interpolation filter simply repeats the last original input samples $N$ times, giving a staircase approximation as shown in the top panel of Figure~\ref{fig:mr:upinterp}.
We know that, in continuous time, the smoothest interpolation is obtained by using a sinc function. This holds in discrete-time as well, and the resulting interpolation filter is the discrete-time sinc:
\begin{equation}
h[n] = \sinc\left(\frac{n}{N}\right)
\end{equation}
Note that the sinc above is equal to one for $n = 0$ and is equal to zero at all integer multiples of $N$, $n = kN$; this fulfills the interpolation condition, that is, after interpolation, the output equals the input at multiples of $N$: $(h \ast x_{N\downarrow})[kN] = x_{N\downarrow}[kN] = x[k]$.
-The three impulse responses above are all lowpass filters; in particular, the sinc interpolator is an ideal lowpass with cutoff frequency $\omega_c = \pi/N$ while the others are approximations of the same. As a consequence, the effect of the interpolator in the frequency domain is the removal of the $N-1$ spectral copies ``drawn in'' the $[-\pi, \pi]$ interval by the upsampler. An example is shown in
-Figure~\ref{fig:mr:upfilt} where the signal in Figure~\ref{upsamplingFigC} is filtered by an ideal lowpass filter with cutoff $\pi/4$.
+The three impulse responses above are all lowpass filters; in particular, the sinc interpolator is an ideal lowpass with cutoff frequency $\omega_c = \pi/N$ while the others are approximations of the same. As a consequence, the effect of the interpolator in the frequency domain is the removal of the $N-1$ spectral copies ``drawn in'' the $[-\pi, \pi]$ interval by the upsampler. An example is shown in Figure~\ref{fig:mr:upfilt} where the upsampled signal is filtered by an ideal lowpass filter with cutoff $\pi/4$.
It turns out that the smoothest possible interpolation in the time domain corresponds to the perfect removal of the spectral repetitions in the frequency domain. Interpolating with a zero-order or a first-order kernel, by contrast, only attenuates the replicas instead of performing a full removal, as we can readily see by considering their frequency responses. Since we are in discrete-time, however, there are no difficulties associated to the design of a digital lowpass filter which closely approximates an ideal filter, so that alternate kernel designs (such as optimal FIRs) can be employed.
This is in contrast to the design of discrete---to continuous---time interpolators, which are analog designs. That is why sampling rate changes are much more attractive in the discrete-time domain.
Downsampling by~$N$\index{downsampling|mie}\index{decimation} (also called subsampling or decimation)
produces a lower-rate sequence by keeping only one out of $N$ samples in the original signal. If we denote by $\mathcal{S}_N$ the downsampling operator\footnote{
We use the letter $\mathcal{S}$ rather than $\mathcal{D}$ since the latter is used for the delay operator.},
Downsampling, as shown in Figure~\ref{fig:mr:down} effectively {\em discards} $N-1$ out of $N$ samples and therefore a loss of information may be incurred. Indeed, in terms of the underlying sampling frequency, decimation produces the signal that would have been obtained by sampling $N$ times more slowly. The potential problems with this data reduction will take the form of aliasing, as we will see shortly.
\subsection{Properties of the Downsampling Operator}
The downsampling operator is linear but it is not time invariant. This can be easily verified with an example; if we subsample by~2 the signal $\mathbf{x}$ we obtain the sequence
One of the consequences of the lack of time-invariance is that complex sinusoids are not eigensequences for the downsampling operator; for instance, if $x[n] = e^{j \pi n} = (-1)^n$, we have
\begin{equation}
\label{eq:mr:12}
x_{2\downarrow}[n] = x[2n] = e^{j 2\pi n} = 1;
\end{equation}
in other words, the highest digital frequency has been mapped to the lowest frequency. This looks very much like aliasing, as we will now explore in detail.
\caption{Downsampling by $4$ in the time domain: original signal (top panel); samples ``killed'' by the downsampler (middle panel); downsampled signal (bottom panel). Note the difference in time indexes between top and bottom panels.}\label{fig:mr:down}
where $\boldsymbol{\xi}_N$ is a ``kill sequence'' defined as
\[
\xi_N[n] =
\begin{cases}
1 & \mbox{ for $n$ multiple of $N$} \\
- 0 & \mbox{ otherwise}
+ 0 & \mbox{ otherwise;}
\end{cases}
\]
-and shown in Figure~\label{fig:mr:xi}; indeed, multiplication by $\boldsymbol{\xi}_N$ will ``kill off'' all the terms in the sum for which the index is not a multiple of $N$. The sequence $\boldsymbol{\xi}_N$ is $N$-periodic and one way to represent it is as the inverse DFS of size $N$ of a vector of all ones, as in~(\ref{eq:fa:unitDFT1}):
+indeed, multiplication by $\boldsymbol{\xi}_N$ will ``kill off'' all the terms in the sum for which the index is not a multiple of $N$. The sequence $\boldsymbol{\xi}_N$ is $N$-periodic and one way to represent it is as the inverse DFS of size $N$ of a vector of all ones, as in~(\ref{eq:fa:unitDFT1}):
we have the scaled sum of $N$ superimposed copies of the original spectrum $X(e^{j\omega})$ where each copy is shifted in frequency by a multiple of $2\pi/N$. We are in a situation similar to that of equation~(\ref{eq:is:periodizedFT}) where sampling created a periodization of the underlying spectrum; here the spectra are already inherently $2\pi$-periodic, and downsampling creates $N-1$ additional interleaved copies. The final spectrum $X_{N\downarrow} (e^{j\omega})$ is simply a stretched version of $A(e^{j\omega})$, so that the interval $[-\pi/N, \pi/N]$ becomes $[-\pi, \pi]$.
Because of the superposition, aliasing\index{aliasing!in multirate} can occur; this is a consequence of the potential loss of information that occurs when samples are discarded. For baseband signals, it is easy to verify that in order for the spectral copies in~(\ref{eq:mr:dss}) not to overlap, the maximum (positive) frequency $\omega_M$ of the original spectrum\footnote{
Here, for simplicity, we are imagining a lowpass real signal whose spectral magnitude is symmetric. More complex cases exist and some examples will be described next.}
must be less than $\pi/N$; this is the \emph{non-aliasing condition} for the downsampling operator. Conceptually, fulfillment of the non-aliasing condition indicates that the discrete-time representation of the original signal is intrinsically redundant; $(N-1)/N$ of the information can be safely discarded and this is mirrored by the fact that only $1/N$ of the spectral frequency support is nonzero. We will see shortly that, in this case, the original signal can be perfectly reconstructed with an upsampling and filtering operation.\index{downsampling|)}
In Figures~\ref{fig:mr:exA} to~\ref{fig:mr:exE}) the top panel shows the original spectrum $X(e^{j\omega})$; the second panel shows the same spectrum but plotted over a wider interval so as to make its periodic nature explicit; the third panel shows in different colors the individual terms in the sum in~(\ref{eq:mr:nonscaled}); the fourth panel shows $A(e^{j\omega})$ \emph{before} scaling and stretching by $N$; finally, the last panel shows $X_{N\downarrow}(e^{j\omega})$ over the usual $[-\pi, \pi]$ interval.
\itempar{Downsampling by 2.} If the downsampling factor is $2$, the \ztrans\ and the Fourier transform of the output are simply
Figure~\ref{fig:mr:exA} shows an example for a lowpass signal whose maximum frequency is $\omega_M = \pi/2$ (i.e.\ a half-band signal). The non-aliasing condition is fulfilled and, in the superposition, the two shifted versions of the spectrum do not overlap. As the frequency axis stretches by a factor of $2$, the original half-band signal becomes full band.
Figure~\ref{fig:mr:exB} shows an example in which the non-aliasing condition is violated. In this case, $\omega_M = 2\pi/3 > \pi/2$ and the spectral copies do overlap. We can see that, as a consequence, the downsampled signal loses its lowpass characteristics. Information is irretrievably lost and the original signal cannot be reconstructed.
\itempar{Downsampling by 3.} If the downsampling factor is $ 3 $ we have
\begin{equation*}
X_{3\downarrow}(e^{j\omega}) = \frac{1}{3}\, \left[ X \bigl(e^{j\frac{\omega}{3}} \bigr) + X \bigl(e^{j(\frac{\omega - 2\pi}{3})} \bigr) + X\bigl(e^{j(\frac{\omega - 4\pi}{3})} \bigr) \right]
\end{equation*}
Figure~\ref{fig:mr:exB} shows an example in which the non-aliasing condition is violated ($\omega_M = 2\pi/3 > \pi/3$).
\itempar{Downsampling a Highpass Signal.} Figure~\ref{fig:mr:exD} shows an example of downsampling by $2$ applied to a half-band {\em highpass} signal. Since the signal occupies only the upper half of the $[0, \pi]$ frequency band (and, symmetrically, only the lower half of the $[-\pi, 0]$ interval), the interleaved copies do not overlap and, technically, there is no aliasing. The shape of the signal, however, is changed by the downsampling operation and what started out as a highpass signal is transformed into a lowpass signal. To make the details of the transformation clearer in this example we have used a {\em complex-valued} highpass signal for which the positive and negative parts of the spectrum have different shapes; it is apparent how the original left and right spectral halves are end up in reverse positions in the final result. The original signal can be exactly reconstructed (since there is no destructive overlap between spectral copies) but the required procedure is a bit more creative and will be left as an exercise.
\subsection{Antialiasing Filters in Downsampling}
We have seen in Section~\ref{sec:is:antialias} that, in order to control the error when sampling a non-bandlimited signal, our best strategy is to bandlimit the signal using a lowpass filter. The same holds when downsampling by $N$ a signal whose spectral support extends beyond $\pi/N$: before downsampling we should apply a lowpass with cutoff $\omega_c = \pi/N$ as shown in Figure~\ref{fig:mr:downfilt}. While a loss of information is still unavoidable, the filtering operation allows us to control said loss and prevent the disrupting effects of aliasing.
An example of the processing chain is shown in Figure~\ref{fig:mr:exE} for a downsampling factor of~$2$; a half-band lowpass \index{half-band filter} filter is used to truncate the signal's spectrum outside of the $[-\pi/2, \pi/2]$ interval and then downsampling proceeds as usual with non-overlapping spectral copies.
\caption{Downsampling by $3$; the highest frequency is larger than $\pi/3$ (here, $\omega_M = 2\pi/3$) and aliasing occurs. Notice how three spectral replicas contribute to the final spectrum.}\label{fig:mr:exC}
\caption{Downsampling by~$2$ of a {\em complex} highpass signal; the asymmetric spectrum helps to understand how non-destructive aliasing works.}\label{fig:mr:exD}
The conversion from one sampling rate to another can always take the ``obvious'' route of interpolating to continuous time and resampling the resulting signal at the desired rate; but this would require dedicated analog equipment and would introduce some quality loss because of the inevitable analog noise. We have just seen, however, that we can increase or decrease the implicit rate of a sequence by an integer factor entirely in the discrete-time domain, by using an upsampler or a downsampler and a digital lowpass filter. For fractional sampling rate changes we simply need to cascade the two operators.
The order in which upsampling and downsampling are performed in a rate changer is crucial since, in general, the operators are not commutative. It is easy to appreciate this fact by means of
Intuitively it's clear that, since no information is lost when using an upsampler, in a fractional sampling rate change the upsampling operation will come first. Typically, a rate change by $N/M$ is obtained by cascading an upsampler by~$N$, a lowpass filter, and a downsampler by~$M$. Since normally the upsampler is followed by a lowpass with cutoff $\pi/N$ while the downsampler is preceded by a lowpass with cutoff $\pi/M$, we can use a single lowpass whose cutoff frequency is the minimum of the two. A block diagram of this system is shown in Figure~\ref{fig:mr:frac}.\index{fractional sampling rate change} \index{rational sampling rate change}
As an example, suppose we want to increase the rate of a sequence originally sampled at 24~KHz up to 32~KHz. For this we need a fractional change of $32/24$ which, after simplification, corresponds to an upsampler by 4 followed by a downsampler by 3, as shown in the top panel of Figure~\ref{fig:mr:fracex}; the lowpass filter's cutoff frequency is $\pi/4$ and, in this case, the lowpass filter acts solely as an interpolator since the overall rate is increased. Conversely, if we want to convert a 32~KHz signal to a 24~KHz, that is, apply a sampling rate change by $3/4$, we can use the cascade shown in the bottom panel of Figure~\ref{fig:mr:fracex}; the cutoff frequency of the filter does not change but the filter, in this case, acts as an anti-alias.
If the ratio between the two sampling frequencies cannot be decomposed into a ratio of small coprime factors, the intermediate rates in a fractional rate changer can be prohibitively high. That was the rationale, for instance, behind an infamous engineering decision taken by the audio industry in the early 90's when the first digital audio recorders (call DAT) were introduced on the market; in order to make it impossible for users to create perfect copies of CDs on digital tapes, the sampling frequency of the recorders was set to 48~KHz. Since CDs are encoded at 44.1~KHz, this required fractional change rate of 160/147. At the time, a 160-fold upsampling was simply not practical to implement so users would have to necessarily go through the analog route to copy CDs. Incidentally, although digital audio tapes have quickly faded into obsolescence, the problem of converting audio from 44.1~KHz to 48~KHz remains relevant today since the sampling frequency used in DVDs is also 48~KHz. The good news is that fractional resamplers can be implemented using local interpolation techniques rather than via a formal upsampling/downsampling chain. To understand the procedure, let's first analyze the practical version of a subsample interpolator and then apply the idea to a resampler.
\itempar{Subsample Interpolation.}
\index{subsample interpolation}
Consider an $F_s$-bandlimited continuous-time signal $\mathbf{x}_c$ and its sampled version $\mathbf{x}$ defined by $x[n] = x_c(nT_s)$, with $T_s \leq 1/F_s$. Given a fractional delay $\tau T_s$, with $|\tau|< 1/2$, we want to determine the sequence
\[
x_\tau[n] = x_c(nT_s + \tau T_s)
\]
using only discrete-time processing; for simplicity, let's just assume $T_s=1$. We know from Section~\ref{sec:is:duality} that the theoretical way to obtain $\mathbf{x}_\tau$ from $\mathbf{x}$ is to use an ideal fractional delay filter:
\[
\mathbf{x}_\tau = \mathbf{d}_\tau \ast \mathbf{x}
\]
where
\begin{align*}
D_\tau(e^{j\omega}) &= e^{j\omega\tau} \\
d_\tau[n] &= \sinc(n - \tau).
\end{align*}
As we have seen in Section~\ref{sec:is:sincinterp}, the sinc interpolator originates as the limit of polynomial interpolation when the number of interpolation points goes to infinity. In this case we can work backwards, and replace the sinc with a low-degree, \textit{local} Lagrange interpolation as in Section~\ref{sec:is:lagrange}. Given $L$ samples to the left and the right of $x[n]$, we can build the continuous-time signal\index{Lagrange interpolation}
where $L_k^{(N)}(t)$ is the $k$-th Lagrange polynomial of order $2N$ defined in~(\ref{eq:is:lagPoly}). With this, we can use the approximation
\begin{equation} \label{eq:mr:subApprox}
x_\tau[n] = l_x(n; \tau) \approx x(n + \tau).
\end{equation}
Figure~\ref{fig:mr:lagPoly} shows for instance the quadratic Lagrange polynomials that would be used for a three-point local interpolation ($L=1$); an example of the interpolation and approximation procedures are shown graphically in Figure~\ref{fig:mr:approx}.
which is the convolution, computed in $n$, between the input signal and the $(2N+1)$-tap FIR
\begin{equation}\label{eq:mr:FIRcoeffs}
\hat{d}_\tau[n] = \begin{cases}
L_k^{(N)}(\tau) & \mbox{for $|n| \leq L$} \\
0 & \mbox{otherwise}
\end{cases}
\end{equation}
For instance, for a quadratic interpolation as in Figure~\ref{fig:mr:approx}, the nonzero coefficients are
\begin{align*}
\hat{d}_\tau [-1] & = \tau \frac{ \tau -1 }{2} \\
\hat{d}_\tau [0] & = - (\tau +1)(\tau -1) \\
\hat{d}_\tau [ 1] & = \tau \frac{\tau +1}{2}
\end{align*}
The FIR interpolator is expressed in noncausal form purely out of convenience; in practical implementations an additional delay would be used to make the whole processing chain causal.
\itempar{Local resampling.} The principle of the subsample interpolator we just described can be used to perform fractional resampling. Again, consider an $F_s$-bandlimited continuous-time signal $\mathbf{x}_c$ and its sampled version $x[n] = x(n T_1)$, with $T_1 \leq 1/F_s$. Given a second sampling period $T_2 \leq 1/F_s$, we want to estimate $y[n] = x(nT_2)$ using only the discrete-time samples $x[n]$. Call
\begin{equation}
\beta = \frac{T_2}{T_1} = \frac{F_1}{F_2}
\end{equation}
the ratio between sampling frequencies; if $\beta < 1$ we are interpolating to a higher underlying rate (i.e. creating more samples overall) whereas if $\beta > 1$ we are downsampling (discarding samples overall). For any output index $n$ we can always write
\begin{equation}%\label{eq:mr:FIRcoeffs}
nT_2 = [m(n) + \tau(n)]\,T_1
\end{equation}
with\footnote{``nint'' denotes the nearest integer function, so that $\mbox{nint}(0.2) = \mbox{nint}(-0.4) = 0$.}
that is, we can associate each output sample $x_2[n]$ to a reference index $m(n)$ in the input sequence plus a fractional delay $\tau(n)$ which is kept less than one half in magnitude. With this, we can use subsample interpolation to approximate each sample
where the coefficients $\hat{d}_{\tau(n)}[k]$, defined in~(\ref{eq:mr:FIRcoeffs}), are now dependent on the output index $n$. In theory, a new set of FIR coefficients is needed for each output sample but, if $F_1$ and $F_2$ are commensurable, we only need to precompute a finite set of interpolation filters. Indeed, if we can write
\[
\beta = \frac{A}{B}, \qquad A, B \in \mathbb{N}
\]
then it is easy to verify from~(\ref{eq:mr:taun}) that
\[
\tau(n + kB) = \tau(n) \quad \mbox{for all $k\in\mathbb{Z}$}.
\]
In other words, there are only $B$ distinct values of $\tau$ that are needed for the subsample interpolation. In the case of our CD to DVD conversion, for instance, we need 147~three-tap filters.
In practical algorithms, which need to work causally, resampling takes place incrementally; this is particularly important in digital communication system design where resampling is vital in order to compensate for slight timing differences between the transmitter's D/A and the receiver's A/D, as we will see in Chapter~\ref{ch:cs}. In general, when $\beta < 1$ we will need to sometimes reuse the same anchor sample with a different $\tau$ in order to produce more samples; conversely, for $\beta > 1$, sometimes we will need to skip one or more input samples. The first case is illustrated in Figure~\ref{fig:mr:resampleup} for $\beta = 0.78$, i.e. for a 28\% sampling rate increase, and the computation of the first few resampled values proceeds as follows:
\begin{itemize}
\item $n=0$: with initial synchronism, $m(0) = 0$, $\tau(0) = 0$ so $y[0] = x[0]$.
The second case is illustrated in Figure~\ref{fig:mr:resampledown} for $\beta = 1.22$ (i.e. an 18\% rate reduction) and the computation of the first few resampled values proceeds as follows:
\begin{itemize}
\item $n=0$: with initial synchronism, $m(0) = 0$, $\tau(0) = 0$ so $y[0] = x[0]$.
The term ``oversampling'' describes a situation in which a signal's sampling rate is deliberately increased beyond the minimum value required by the sampling theorem in order to improve the performance of A/D and D/A converters at a lower cost than would be required by the use of more sophisticated analog circuitry.
The sampling theorem guarantees that we can losslessly convert an $F_s$-band\-limited signal $\mathbf{x}_c$ into a discrete-time sequence, provided that the sampling frequency is larger than $F_s$. In a full A/D conversion, therefore, the only remaining source of distortion is introduced by quantization. Assuming a uniformly distributed input and a matching uniform quantizer, in Section~\ref{sec:da:quantization}, we have modeled this distortion as an additive noise source,
\begin{equation}\label{eq:mr:overAD}
\hat{\mathbf{x}} = \mathbf{x + e},
\end{equation}
where $\mathbf{e}$ is a white process, \textit{uncorrelated with $\mathbf{x}$}, and whose uniform power spectral density
depends only on $\Delta$, the width of quantization interval or, equivalently, on the number of bits per sample. This is represented pictorially in the top panel of Figure~\ref{fig:mr:overAD} which shows the PSDs of a critically sampled signal and that of the associated quantization noise; the total SNR is the ratio between the areas under the two curves. One way to improve the SNR, as we know, is to increase $R$, the number of bits per sample used by the quantizer; unfortunately, the number of electronic components in an A/D converter grows proportionally to $R^2$, and so this is an expensive option.
A clever digital workaround is provided by the observation that the sampling frequency used to convert $\mathbf{x}_c$ into a digital signal does not appear explicitly in~(\ref{eq:mr:overAD}) and therefore if we oversample the analog signal by, say, a factor of two, we will ``shrink'' the support of the signal's PSD without affecting the noise; this is shown in the middle panel of Figure~\ref{fig:mr:overAD}. Increasing the sampling rate does not modify the power of the signal, so the overall SNR does not change: in the figure, note how the shrinking support is matched by a proportional increase in amplitude for the signal's PSD\footnote{
Although we haven't explicitly introduced a sampling theorem for random processes, the formula for the spectrum of a sampled signal in~(\ref{eq:is:noAliasingSpecEq}) formally holds for power spectral densities as well, and the change in amplitude is due to the sampling frequency appearing as a scaling factor.}.
At this point, however, we are in the digital domain and therefore it is simple (and cheap) to filter the out-of-band noise with a sharp lowpass filter with a magnitude response close to the dashed line in Figure~\ref{fig:mr:overAD}. In the example, this will halve the total power of the noise, increasing the SNR by a factor of 2 (that is, by 3~dB). We can now digitally downsample the result to obtain a critically sampled input signal with an improved SNR as illustrated in the bottom panel of Figure~\ref{fig:mr:overAD}. To foster the intuition, we can look at the process from the time domain: at a high sampling rate, neighboring samples can be seen as repeated measurements of the same value (the signal varies slowly compared to the speed of sampling) affected by quantization noise that is supposed to be independent from sample to sample; the lowpass filter acts as a local averager which reduces the variance of the noise for the subsampled sequence.
\caption{Oversampling for A/D conversion: signal's PSD and quantization error's PSD (gray) for critically sampled signal (top panel); oversampled signal and lowpass filter (middle panel); filtered and downsampled signal (bottom panel).}\label{fig:mr:overAD}
In general, the block diagram for an oversampled A/D is as in Figure~\ref{fig:mr:overADblock} and, theoretically, the SNR of the quantized signal improves by 3~dB for every doubling of the sampling rate with respect to the baseline SNR provided by the quantizer. However, as we just illustrated with a time-domain analysis, the technique is predicated on two fundamental assumptions, the statistical independence between signal and noise and the absence of correlation between successive noise samples, both of which become invalid as the sampling rate increases. With high oversampling factors the correlation between successive noise samples increases and therefore most of the quantization noise will leak in the band of the signal. More efficient oversampling technique use feedback and nonlinear processing to push more and more of the quantization noise out of band; these converters, known under the name of Sigma-Delta quantizers, are however very difficult to analyze theoretically.
All practical D/A converters use a kernel-based interpolator; recall from Section~\ref{sec:is:locInterp} that the interpolated continuous-time signal in this case is
with, as usual, $F_s = 1/T_s$. The above expression is the product of two terms; the first is the periodic digital spectrum, rescaled so that $\pi \rightarrow F_s/2$, and the second is the frequency response of the interpolation kernel.
An ideal D/A converter would require a sinc interpolator, which we know not to be realizable in practice, but which would completely remove the out-of-band components from the spectrum interpolated signal as shown in Figure~\ref{fig:mr:sincInterp}, where the top panel shows a digital spectrum, the middle panel shows the two terms of Equation~(\ref{eq:mr:interpSpec}), and the bottom panel shows the final analog spectrum.
At the other end of the interpolation gamut lies the zero-order hold which, as we have seen, is easy to implement but has terrible spectral properties; the problems are shown in Figure~\ref{fig:mr:sincInterp}, in which the contents of the three panels are as in Figure~\ref{fig:mr:sincInterp}. The ZOH-interpolated spectrum is affected in two ways:
\begin{enumerate}
\item there is significant out-of-band energy due to the spectral copies that are only attenuated and not eliminated by the interpolator;
\item there is in-band distortion due to the fact that the interpolator has a non-flat passband.
\end{enumerate}
Yet, the zero-order-hold is so easy and cheap to implement that most D/A circuits use it exclusively.
In theory, to compensate for the resulting problems, we would need an \textit{analog} filter whose frequency response is sketched in Figure~\ref{fig:mr:zohComp}; such a filter, however, even assuming we could design it exactly, would be a costly device since an analog filter with a sophisticated response requires high-precision electronic components.
Again, rather than using complicated and expensive analog filters, the performance of the D/A converter can be dramatically improved if we are willing to perform the conversion at a higher rate than the nominal sampling frequency, as shown in Figure~\ref{fig:mr:overDAblock}. The digital signal is first oversampled by a factor of $N$ and then filtered with a discrete-time filter $F(e^{j\omega})$ that combines the out of band frequency rejection of a sharp lowpass with cutoff $\pi/N$ with an in-band frequency shaping that mimics the characteristic of the filter in Figure~\ref{fig:mr:zohComp}. Since the filter is digital, there is no difficulty in designing, say, an unconditionally stable FIR with the desired characteristic.
The oversampled D/A procedure using a zero-order hold and an oversampling factor $N=2$ is illustrated in Figure~\ref{fig:mr:overDA}. The top panels shows the DTFT of the original signal; the spectrum that enters the interpolator is
We use an ideal filter for convenience but of course, in practical implementations, this would be a realizable sharp lowpass.}
matches the upsampler's rate and where $C(e^{j\omega})$ compensates for the nonflat in-band characteristics of the interpolator; the resulting spectrum is shown in the third panel. Note how oversampling has created some ``room'' to the left and the right of the spectrum; this will be important for the analog part of the D/A. If we now interpolate $x_o[n]$ at a rate of $F_o = NF_s$~Hz, we have
\begin{align}
X_o(f) &= \frac{1}{F_o}\,X_o(e^{j2\pi f / F_o})\, I_0\left(\frac{f}{F_o}\right) \nonumber \\
The fourth panel in Figure~\ref{fig:mr:overDAblock} shows the digital spectrum mapped on the real frequency axis and the magnitude response of the zero-order hold; note that now the first zero crossing for the latter occurs at $NF_s$ (compare that to Figure~\ref{fig:mr:zohInterp}). Since $C(e^{j\omega})$ is designed to compensate for $I_0(f)$, we have that
\[
X_o(f) = X(f) \quad \mbox{for $f\in [-F_s, F_s]$}
\]
that is, the interpolation is equivalent to a sinc interpolation over the baseband. The only remaining problem is the spurious high frequency content at multiples of $F_o$, as shown in the last panel of Figure~\ref{fig:mr:overDAblock}. This needs to be eliminated with an analog filter but, because of the room between spectral replicas created by oversampling, the required filter can have a wide transition band as shown in Figure~\ref{fig:mr:overDAlast} and therefore can be implemented very cheaply using for instance a low-order Butterworth lowpass. Finally, note that the present analysis remains valid also with higher order kernels, such as the first- and third-order interpolators detailed in Section~\ref{sec:is:locInterp}, since their frequency response is similar to that of the zero-order hold.
+ \item Downsampling by $2$ followed by filtering by $H(z)$ is equivalent to filtering by $H(z^2)$ followed by downsampling by $2$.
+ \item Filtering by $H(z)$ followed by upsampling by $2$ is equivalent to upsampling by $2$ followed by filtering by $H(z^2)$.
+ \end{enumerate}
\end{exercise}
+}\fi
+\ifanswers{%
+\begin{solution}{}
+ We will operate in the $z$-domain; recall that, given the \ztrans\ of the original sequence $X(z)$, the \ztrans\ of a 2-upsampled sequence is $X(z^2)$ whereas the \ztrans\ of a 2-downsampled sequence is $[X(z^{1/2}) + X(-z^{1/2})]/2$. With this:
+ \begin{enumerate}
+ \item Downsampling by 2 followed by filtering by $H(z)$ yields
+ The two operations are thus equivalent. (Please note in $H(z^2)$ how it's the \textit{argument} $z$ of the $z$-transform that is replaced by $\pm z^{1/2}$, so that the minus sign goes away. This makes sense: if a sequence has been upsampled, downsampling simply returns it to its original form.)
+ \item Filtering by $H(z)$ followed by upsampling by 2 can be written as
+ \[
+ Y_1(z)= H(z^2)X(z^2).
+ \]
+ Upsampling by 2 followed by filtering by $H(z^2)$ can be written as
-Consider the input-output characteristic of the following multirate systems. Remember that, technically, one cannot talk of transfer functions in the case of multirate systems since sampling rate changes are not time-invariant. It may happen, though, that by carefully designing the processing chain, the input-output characteristic does indeed implement a time-invariant transfer function.
-\begin{enumerate}
- \item Find the overall transformation operated by the following system:
+ Consider the input-output characteristic of the following multirate systems. Remember that, technically, one cannot talk of transfer functions in the case of multirate systems since sampling rate changes are not time-invariant. It may happen, though, that by carefully designing the processing chain, the input-output characteristic does indeed implement a time-invariant transfer function.
+ \begin{enumerate}
+ \item Find the overall transformation operated by the following system:
+ \item Assume $H(z) = A(z^2)+z^{-1}B(z^2)$ for arbitrary $A(z)$ and $B(z)$. Show that the transfer function of the following system is equal to $A(z)$.
+ The lower branch contains an upsampler by two followed by a delay and a downsampler by two. The output of such a processing chain is easily seen to be equal to $0$ for all $n$. Thus only the upper branch remains and the final transfer function of the system is given by
-\begin{exercise}{Digital processing of continuous-time signals}
-In your grandmother's attic you just found a treasure: a collection of super-rare $78$ rpm vinyl\index{vinyl} jazz records. The first thing you want to do is to transfer the recordings to compact discs, so you can listen to them without wearing out the originals. Your idea is obviously to play the record on a turntable and use an A/D converter to convert the line-out signal into a discrete-time sequence, which you can then burn onto a CD. The problem is, you only have a ``modern'' turntable, which plays records at $33$ rpm. Since you're a DSP wizard, you know you can just go ahead, play the $78$ rpm record at $33$~rpm and sample the output of the turntable at $44.1$~KHz. You can then manipulate the signal in the discrete-time domain so that, when the signal is recorded on a CD
-and played back, it will sound right.
-
-In order to design a system that performs the above conversion consider the following:
-\begin{itemize}
- \item Call $s(t)$ the continuous-time signal encoded on the $78$ rpm vinyl (the jazz music).
- \item Call $x(t)$ the continuous-time signal you obtain when you play the record at $33$ rpm on the modern turntable.
- \item Let $x[n] = x(nT_s)$, with $T_s = 1/44,100$.
-\end{itemize}
-Answer the following questions:
-\begin{enumerate}
- \item Express $x(t)$ in terms of $s(t)$.
- \item Sketch the Fourier transform $X(f)$ when $S(f)$ is like so:
+\begin{solution}{}
+ Intuitively, the upsampling by four followed by $I_1(z)$ creates a linear interpolation over three ``extra'' samples between the original values the (``connect-the-dots'' strategy). The delay by two followed by the downsampler selects the midpoint of each interpolation interval. As a whole, the chain implements a fractional delay of half a sample using a linear interpolator so that
+ Since the interpolator $I(z)$ has finite support of length 7, we can concentrate on the interval $[-4, 4]$ and extend the result to the other points. Linear interpolation fills in the gaps while the delay shifts the interpolated signal by two towards the right:
- The highest nonzero frequency of $S(f)$ is $F_{\max} = 16,000$~Hz (old records have a smaller bandwidth than modern ones).
- \item Design a system to convert $x[n]$ into a sequence $y[n]$ so that, when you interpolate $y[n]$ to a continuous-time signal $y(t)$ with interpolation period $T_s$, you obtain $Y(f) = S(f)$.
- \item What if you had a turntable which plays records at $45$ rpm? Would your system be different? Would it be better?
-\end{enumerate}
+
+ The downsampler selects the points in $c[n]$ where $n$ is a multiple of four, which are the midpoints between original data values:
+\begin{exercise}{Digital processing of continuous-time signals}
+ In your grandmother's attic you just found a treasure: a collection of super-rare $78$ rpm vinyl\index{vinyl} jazz records. The first thing you want to do is to transfer the recordings to digital format, so you can listen to them without wearing out the originals. Your idea is obviously to play the record on a turntable and use an A/D converter to convert the line-out signal into a discrete-time sequence, which you can then store or burn onto a CD. The problem is, you only have a ``modern'' turntable, which plays records at $33$ rpm. Since you're a DSP wizard, you know you can just go ahead, play the $78$ rpm record at $33$~rpm and sample the output of the turntable at $44.1$~KHz. You can then manipulate the signal in the discrete-time domain so that, when the signal is recorded on a CD and played back, it will sound right.
+
+ In order to design a system that performs the above conversion consider the following:
+ \begin{itemize}
+ \item Call $s(t)$ the continuous-time signal encoded on the $78$ rpm vinyl (the jazz music).
+ \item Call $x(t)$ the continuous-time signal you obtain when you play the record at $33$ rpm on the modern turntable.
+ \item Let $x[n] = x(nT_s)$, with $T_s = 1/44,100$.
+ \end{itemize}
+ Answer the following questions:
+ \begin{enumerate}
+ \item Express $x(t)$ in terms of $s(t)$.
+ \item Sketch the Fourier transform $X(f)$ when $S(f)$ is like so:
+ The highest nonzero frequency of $S(f)$ is $F_{\max} = 16,000$~Hz (old records have a smaller bandwidth than modern ones).
+ \item Design a system to convert $x[n]$ into a sequence $y[n]$ so that, when you interpolate $y[n]$ to a continuous-time signal $y(t)$ with interpolation period $T_s$, you obtain $Y(f) = S(f)$.
+ \item What if you had a turntable which plays records at $45$ rpm? Would your system be different? Would it be better?
+ \end{enumerate}
\end{exercise}
+}\fi
+
+\ifanswers{%
+\begin{solution}{}
+ \begin{enumerate}
+ \item Playing the record at lower rpm slows the signal down by a factor $33/78$. Therefore
+ \item We need to change the sampling rate so that, when $y[n]$ is interpolated at 44.1~KHz its spectrum is equal to $S(f)$. The rational sampling rate change factor is clearly $33/78$ which is simply $11/26$ after factoring. The processing scheme is as follows:
+ where $L(z)$ is a lowpass filter with cutoff frequency $\pi/26$ and gain $L_0 = 11$; both the sampler and interpolator work at $F_s = 1/T_s = 44100$~Hz. We have:
+ \item The sampling rate change scheme stays the same except that now $45/78 = 15/26$. Therefore, the final upsampler has to compute more samples than in the previous scheme. The computational load of the sampling rate change is entirely dependent on the filter $L(z)$. If we upsample more before the output, we need to compute more filtered samples and therefore at 45rpm the scheme is less efficient.
-and show that this system implements a fractional delay (i.e.\ show that the transfer function of the system is that of a pure delay, where the delay is not necessarily an integer).
-
-To see a practical use of this structure, consider a data transmission system over an analog channel. The transmitter builds a discrete-time signal $s[n]$; this is converted to an analog signal $s_c(t)$ via an interpolator with period $T_s$, and finally $s_c(t)$ is transmitted over the channel. The signal takes a finite amount of time to travel all the way to the receiver; say that the
-transmission time over the channel is $t_0$~seconds: the received signal $\hat{s}_c(t)$ is therefore just a delayed version of the transmitted signal,
-\[
- \hat{s}_c(t) = s_c (t-t_0)
-\]
-At the receiver, $\hat{s}_c(t)$ is sampled with a sampler with period $T_s$ so that no aliasing occurs to obtain $\hat{s}[n]$.
-\begin{enumerate}
- \item Write out the Fourier Transform of $\hat{s}_c(t)$ as a function of $S_c(f)$.
- \item Write out the DTFT of the received signal sampled with rate $T_s$, $\hat{s}[n]$.
- \item Now we want to use the above multirate structure to compensate for the transmission delay. Assume $t_0 = 4.6\,T_s$; determine the values for $M$ and $L$ in the above block diagram so that $\hat{s}[n] = s[n - D]$, where $D \in \mathbb{N}$ has the smallest possible value (assume an ideal lowpass filter in the multirate structure).
+ and show that this system implements a fractional delay (i.e.\ show that the transfer function of the system is that of a pure delay, where the delay is not necessarily an integer).
+
+ To see a practical use of this structure, consider a data transmission system over an analog channel. The transmitter builds a discrete-time signal $s[n]$; this is converted to an analog signal $s_c(t)$ via an interpolator with period $T_s$, and finally $s_c(t)$ is transmitted over the channel. The signal takes a finite amount of time to travel all the way to the receiver; say that the
+ transmission time over the channel is $t_0$~seconds: the received signal $\hat{s}_c(t)$ is therefore just a delayed version of the transmitted signal,
+ \[
+ \hat{s}_c(t) = s_c (t-t_0)
+ \]
+ At the receiver, $\hat{s}_c(t)$ is sampled with a sampler with period $T_s$ so that no aliasing occurs to obtain $\hat{s}[n]$.
+ \begin{enumerate}
+ \item Write out the Fourier Transform of $\hat{s}_c(t)$ as a function of $S_c(f)$.
+ \item Write out the DTFT of the received signal sampled with rate $T_s$, $\hat{s}[n]$.
+ \item Now we want to use the above multirate structure to compensate for the transmission delay. Assume $t_0 = 4.6\,T_s$; determine the values for $M$ and $L$ in the above block diagram so that $\hat{s}[n] = s[n - D]$, where $D \in \mathbb{N}$ has the smallest possible value (assume an ideal lowpass filter in the multirate structure).
-One sample of $x[n]$ may have been corrupted and we would like to approximately or exactly recover it. We denote $n_0$ the time
-index of the corrupted sample and $\hat{x}[n]$ the corresponding corrupted sequence.
-\begin{enumerate}
- \item Specify a practical algorithm for exactly or approximately recovering $x[n]$ from $\hat{x}[n]$ if $n_0$ is known.
- \item What would you do if the value of $n_0$ is not known?
- \item Now suppose we have $k$ corrupted samples at either known or unknown locations. What is the condition that $X(e^{j \omega })$ must satisfy to be able to exactly recover $x[n]$? Specify the algorithm.
+ One sample of $x[n]$ may have been corrupted and we would like to approximately or exactly recover it. We denote $n_0$ the time index of the corrupted sample and $\hat{x}[n]$ the corresponding corrupted sequence.
+ \begin{enumerate}
+ \item Specify a practical algorithm for exactly or approximately recovering $x[n]$ from $\hat{x}[n]$ if $n_0$ is known.
+ \item What would you do if the value of $n_0$ is not known?
+ \item Now suppose we have $k$ corrupted samples at either known or unknown locations. What is the condition that $X(e^{j \omega })$ must satisfy to be able to exactly recover $x[n]$? Specify the algorithm.
Usually, discrete-time signals are introduced as the discretized version of real-world signals and this is because, traditionally, the mathematical models for natural phenomena are provided by functions of a continuous-time variable. Nowadays, however, the vast majority of signals are de facto ``born'' in discrete time and the actual act of sampling an analog measurement is made transparent by the ubiquity of digital acquisition devices. Furthermore, especially in the context of communication systems, we are just as concerned with the \emph{synthesis} of discrete-time signals as with their acquisition. With digital processing at the heart of virtually every professional and consumer device today, the discrete-time paradigm has largely supplanted its analog ancestors and discrete-time signals can therefore be introduced on their own merit from an abstract and self-contained point of view.
\section{On DSP Notation}
Discrete-time signal processing, as a discipline, sits at a crossroad between mathematics (think harmonic analysis), electrical engineering (think circuit design) and computer science (DSP initially appeared as a way to simulate analog electronics on a digital computer). Along the years, scholars and researchers have naturally borrowed a wide range of notational conventions (and associated quirks) from their domain of origin with the result that, sometimes, coherence and clarity are not exactly paramount in the signal processing literature.
The first problem concerns how to notate discrete-time signals, that is, \textit{sequences} of values indexed by an integer variable. A common convention is to use square brackets to highlight the fact that the ``time'' index is in fact an integer; as an example, the expression
\begin{equation}
x[n] = \cos(2\pi n / 3)
\end{equation}
clearly indicates the discrete-time nature of the sequence with square brackets, and uses the standard convention of parentheses for the argument of the trigonometric function (a real-valued quantity). The enduring advantage of this notation is that it underscores the algorithmic nature of DSP, since square brackets are used by most programming languages to select an element in an array. The problems begin to arise when, as is often the case, the expression $x[n]$ is used to denote the \textit{entire} signal and not just its value in $n$. This abuse of notation is hardly a prerogative of signal processing: in calculus it's very common to denote an entire function by $f(t)$, an expression that should only be used to indicate the value of $f$ in $t$. But in signal processing the problem is a bit more delicate: finite-length signals, for instance, are defined only for a precise subset of indexes (e.g., for $0 \leq n < N$); in these cases, we need to use extra care to indicate the legal domain for the expression $x[n]$ and the normal way to express a delay, $x[n-M]$, is in fact largely meaningless without additional provisos. Another source of problems is represented by the recurrent use of signal \textit{operators}, such as the Discrete Fourier Transform; a commonly used notation is
\[
X[k] = \mbox{DFT}\{x[n]\}
\]
but of course it's not clear what the indexes on the left and right-hand side of the expression stand for: are they just placeholders or are they related? Finally, we owe a special mention to the biggest source of confusion, the convolution operator, which is often notated as
\begin{equation}
y[n] = h[n] \ast x[n]. \label{eq:in:conv1}
\end{equation}
If we write the value of the \textit{sample} $y[n]$ explicitly,
it's clear that $n$ in $y[n]$ is the ``anchor point'' for the computation of the summation; but what is the meaning of $n$ in the right-hand side of~(\ref{eq:in:conv1})? Which of the following expressions is correct?
\begin{align*}
y[n-2] &= h[n - 2] \ast x[n-2] \\
y[n-2] &= h[n - 1] \ast x[n-1]
\end{align*}
In this book we are not going to try and fix all of these problems by introducing a new set of notational conventions; after all, the main purpose of a textbook is to prepare its readers so that they can ultimately tackle the ``raw'' research literature available on its subject, and part of this training includes the exposure to established practices. We will be rather pragmatic, and follow what we feel is a good compromise between ``tradition'' and clarity:
\begin{itemize}
\item for us, all discrete-time signals can be considered as vectors in an appropriate vector space. Therefore, we will use boldface characters when indicating the entire signal\footnote{
The real practical problem that prevented a wider adoption of vector notation is that it's really inconvenient to draw boldface characters by hand; while alternate and more handwriting-friendly vector symbols exist (e.g. $\vec{x}$), they are still ``busy-looking'' and can be easily mistook for other common symbols (e.g. $\bar{x}$ for a mean or $\hat{x}$ for an estimate). On the other hand, once the use of $\mathbf{x}$ is established to indicate the whole signal, one needs not be too strict with the boldface convention and the meaning of the lone symbol $x$ will be entirely clear.}
and normal type with square brackets to specify the value of a signal at a given discrete time:
\item signal operators will be indicated with calligraphic letters: the notation
\[
\mathbf{y} = \mathcal{H}\mathbf{x}
\]
indicates that $\mathbf{y}$ is the result of applying the operator $\mathcal{H}$ (e.g. a delay, a filter, or a transform) to the input $\mathbf{x}$. For finite-length signals and linear operators, this notation exactly mirrors the fact that the operator can be implemented as a matrix-vector multiplication. Sometimes, to enhance the algorithmic nature of some operators, we may choose to use explicit names, as in the expression
\item for the Discrete-Time Fourier Transform we will stick to the common notation that uses a complex exponential as the argument of the function, $X(e^{j\omega})$; the Fourier transform of a continuous-time function $x(t)$, on the other hand, will be denoted as $X(f)$, with $f$ expressed in Hertz.
A discrete-time signal\index{discrete-time signal|mie} is a complex-valued \emph{sequence}\index{sequence|see{discrete-time signal}}, that is, a complex-valued function of an \textit{integer} index $n \in \mathbb{Z}$; in its most general form, it is a two-sided, infinite and countable collection of values. As explained in Section~\ref{sec:in:notation}, we will use boldface characters to denote the entire signal and square brackets to indicate a specific element. For instance, we can define a sequence $\mathbf{x}$ analytically as:
which is a complex exponential of period~$40$ samples, plotted in Figure~\ref{fig:dt:compEx}. Another example, this time of a random sequence whose samples are uniformly distributed between $-1$ and $1$, is
\begin{equation}\label{eq:dt:rand}
x[n] = \mbox{the $n$-th output of a random source } \mathcal{U}(-1,1)
\end{equation}
a realization of which is plotted in Figure~\ref{fig:dt:rand}. Finally, an example of a sequence drawn from the real world and for which, therefore, no analytical expression exists, is
\begin{equation}\label{eq:dt:dow}
x[n] = \mbox{The average Dow-Jones index in year }n
\end{equation}
plotted in Figure~\ref{fig:dt:dow} from year~1900 to~2007.
\item The sequence index $n$ is best thought of as a measure of \emph{dimensionless time}; while it has no physical unit of measure, it imposes a chronological order on the values of the sequences.
\item We consider complex-valued discrete-time signals; while physical signals can be expressed by real quantities, the generality offered by the complex domain is particularly useful in designing systems which \emph{synthesize} signal, such as data communication systems.
\item In graphical representations, when we need to emphasize the discrete-time nature of the signal, we resort to stem (or {}``lollipop''{}\index{lollipop plot}) plots as in Figure~\ref{fig:dt:triangle}; this works well for plots that show a small number of data points. For lager data sets, when the discrete-time domain is implicitly understood, we will often use a function-like representation as in Figure~\ref{fig:dt:dow}. In the latter case, each ordinate of the sequence is graphically connected to its neighboring data points, giving the illusion of a smooth line; while this makes the plot easier on the eye, it must be remembered that the signal is defined only over a \emph{discrete} set of abscissas.
The examples described by~(\ref{eq:dt:triangle}) or~(\ref{eq:dt:compEx}) represent two-sided, infinite sequences. Of course, in the practice of signal processing, it is impossible to deal with an infinite
amount of data: for a processing algorithm to execute in a finite amount of time and to use a finite amount of storage, the input must be of finite length; even for algorithms that operate on the fly, i.e.\ algorithms that produce an output sample for each new input sample, an implicit limitation on the input data size is imposed by the necessarily limited life span of the processing device.\footnote{Or, in the extreme limit, of the supervising engineer $\ldots$} This limitation was clearly implicit in our attempts to plot these sequences, as we did in Figure~\ref{fig:dt:triangle} and~\ref{fig:dt:compEx}: what the diagrams show, in fact, is only a representative, meaningful portion of the signals, and we rely on the original formulas for their full description.
Signals sampled from the real world, for which no analytical description exists, always possess a finite time support because of the finite time spent recording the underlying phenomena; for the Dow Jones index, for instance, the data does not exist for years before 1884, when the index was introduced, and future values are certainly not known. More importantly (and more often), the finiteness of a discrete-time signal is explicitly imposed by design since we are interested in concentrating our processing efforts on a small portion of an otherwise longer signal; in a speech recognition system, for instance, the practice is to cut a speech signal into small segments and try to identify the phonemes associated to each one of them.\footnote{
Note that, in the end, phonemes are pasted together into words and words into sentences; therefore, for a complete speech recognition system, long-range dependencies become important again.}
These finite sets of data are not, strictly speaking, sequences in the mathematical sense, although they can be formally extended into one, as we will see.
Another case is that of \textit{periodic} signals, as in~(\ref{eq:dt:compEx}); even though these are indeed infinite sequences, it is clear that all the relevant information is contained in just a single period. By
describing one period (graphically or otherwise), we are, in fact, providing a full description of the entire sequence. This variety of signal types demands a little taxonomic effort that we will now undertake.
A finite-length discrete-time signal of length~$N$ is simply a collection (or a \textit{tuple}) of $N$ complex values. It is both natural and convenient to consider finite-length signals as Euclidean vectors in the $N$-dimensional space $\mathbb{C}^N$, and we will explore this viewpoint in detail in the next chapter. The key point for now is that the object that represents an $N$-point, finite-length signal is the vector
note the transpose operator that shows the convention of using \emph{column} vectors. The individual element of the signal is denoted by the expression
\[
x[n], \qquad n = 0,\, \ldots,\, N-1;
\]
the range of the index $n$ is finite and, for values outside of the $[0, N-1]$ interval, the notation $x[n]$ is \emph{not defined} .
As we said, finite-length signals are the only actual entities that we can manipulate in practical signal processing applications; it would be extremely awkward, however, to develop the whole theory of signal processing only in terms of finite-length signals. In fact, it is more convenient, when possible, to develop theoretical proofs using infinite sequences since such general results will hold for all finite-size signals as well.
Infinite-length signals are proper sequences in the mathematical sense; although these kind of signals in general lie beyond our processing and storage capabilities, they are useful from a theoretical point of view since they allow us to prove general theorems that apply to all other signal classes. We can define three subclasses as follows.
\itempar{Aperiodic Signals.} The most general type of discrete-time signal is represented by an aperiodic, two-sided complex sequence. In the theory of signal processing, we tend to like like sequences that are in some way \textit{summable}, since summability is associated with desirable mathematical properties. In particular, \emph{absolute summability} (a strong condition) is associated with system stability, while \emph{square summability} (a weaker condition) is associated to finite energy. In general, theoretical proofs will be easier when the signals involved are assumed to be absolutely summable\footnote{
Mostly because of Fubini's theorem, thanks to which we can freely reorder summations and integrals.}
whereas extending the proofs to square summable sequences is often quite a bit of work.
An $N$-periodic sequence $\tilde{\mathbf{x}}$ is an infinite, two-sided sequence for which
\begin{equation}
\tilde{x}[n] = \tilde{x}[n + pN], \quad\qquad p \in \mathbb{Z}.
\end{equation}
The tilde notation will be used whenever we want to explicitly stress a periodic behavior. Clearly an $N$-periodic sequence is completely defined by its $N$ values over a period; that is, a periodic sequence {}``carries no more information''{} than a finite-length signal of length $N$, in spite of its infinite support.
\index{periodic extension}
In this sense, periodic signals represent a first possible bridge between finite-length signals and infinite sequences. Indeed, we can always {}``convert''{} a finite-length signal $\mathbf{x} \in \mathbb{C}^N$ into a sequence via its periodic extension $\tilde{\mathbf{x}}$, built as
\begin{equation}
\tilde{x}[n] = x[n \mod N], \qquad \quad n \in \mathbb{Z}.
\end{equation}
We will soon see that this type of embedding for a finite-length signal is the {}``natural''{} one in most signal processing contexts.
A discrete-time sequence $\bar{\mathbf{x}}$ is said to have \emph{finite support} if its values are zero for all
indexes outside of a given range; that is, there exist two values~$N \in \mathbb{N}^+$ and $M \in \mathbb{Z}$ such that
\[
\bar{x}[n] = 0 \qquad \mbox{for } n < M \mbox{ or } n > M + N -1;
\]
Note that, although $\bar{\mathbf{x}}$ is an infinite sequence, knowledge of its $N$~nonzero values (and of the start time $M$) completely determines the entire signal. This suggests another way to embed a finite-length signal into a sequence $\bar{\mathbf{x}}$, by simply extending the original data with zeros to the left and to the right:
\begin{equation}\label{FiniteSupportExt}
\bar{x}[n] = \begin{cases}
x[n] & \mbox{ if } 0 \leq n < N - 1 \\
0 & \mbox{ otherwise}
\end{cases}
\qquad \quad n \in \mathbb{Z}.
\end{equation}
In general, we will use the bar notation $\bar{\mathbf{x}}$ for sequences defined as the finite support extension of a finite-length signal.
The following discrete-time signals are basic building blocks that will appear repeatedly throughout the rest of the book. The signals are defined for all $n \in \mathbb{Z}$ and are therefore two-sided infinite sequences; the expressions can however be used to define finite-length signals as well, simply by restricting the range of the index to the interval $[0, N-1]$.
\itempar{Impulse.}\index{impulse}
The discrete-time impulse $\boldsymbol{\delta}$ (also known, from its symbol, as the discrete-time \emph{delta}\footnote{Not to be confused with the Dirac delta functional, that we will encounter later.}) is the simplest discrete-time signal and it is defined as:
\begin{equation}
\delta[n] = \begin{cases}
1 & n = 0 \\
0 & n \neq 0;
\end{cases}
\end{equation}
the signal is shown in Figure~\ref{fig:dt:delta}. The delta signal has a single nonzero sample and is therefore the most ``concentrated'' signal in time; a shifted version of the delta, in which the nonzero sample occurs for $n=k$, is indicated by $\boldsymbol{\delta}_k$. The discrete-time delta models a physical phenomenon with the shortest possible duration.
an example is shown in Figure~\ref{fig:dt:expDecay}. The exponential decay is an important signal since it models the response of a discrete-time first order recursive filter, as we will study in detail in Chapter~\ref{ch:lti}. Exponential sequences are ``well-behaved'' only for values of $a$ less than one in magnitude; sequences in which $| a | > 1$ grow unbounded and represent an unstable behavior.
A complex exponential signal is defined by the expression:
\begin{equation}
x[n] = e^{j(\omega n + \phi)}
\end{equation}
where $\omega$ is the \textit{frequency} and $\phi$ is the initial \textit{phase}; both phase and frequency are dimensionless quantities expressed in radians. The complex exponential describes the prototypical oscillatory behavior in signal processing, as we will study in detail in Chapter~\ref{ch:fa}. The standard real-valued oscillations can always be obtained via Euler's formula, that is:
\begin{equation}
e^{j(\omega n + \phi)} = \cos(\omega_0 n + \phi) + j \sin(\omega_0 n + \phi)
\end{equation}
as shown by the example in Figure~\ref{fig:dt:compEx}.
The complex exponential (and, equivalently, its real-valued sinusoidal components) expresses an oscillatory behavior parametrized by a frequency and, possibly, an initial phase. In the continuous-time world, where time is measured in seconds ($s$), frequency has a physical dimension of $1/s$ and is measured in Hertz (Hz) so that a sinusoidal oscillation with frequency $f_0$~Hz is described by the function
\[
f(t) = \cos(2\pi f_0 t),
\]
where the $2\pi$ factor converts the argument of the trigonometric function to radians. In discrete time, where the index $n$ represents dimensionless time, the {}``digital''{} frequency is expressed in radians, itself a dimensionless quantity.\footnote{
A radian is dimensionless since it the ratio of two measures of length, that of the arc on a circle
subtended by the measured angle and the radius of the circle itself.}
The best way to appreciate how a discrete-time oscillation works is to consider a generic algorithm to generate successive samples of a discrete-time sinusoid at a given digital frequency $\omega_0$; from the formula
\[
x[n] = \cos(\omega_0 n + \phi)
\]
we can see that the argument of the trigonometric function (that is, the instantaneous phase) is incremented by $\omega_0$ at each step, yielding:
With this in mind, it is easy to see that \emph{the highest frequency of a discrete-time oscillator is $\omega_{\max} = 2\pi$}; for any value larger than this, the inner $2\pi$-periodicity of the trigonometric functions {}``maps back''{} the output values to a frequency between $0$ and $2\pi$:
for all values of $k \in \mathbb{Z}$. This modulo-$2\pi$ equivalence of digital frequencies, also known as \textit{aliasing}, is a pervasive concept in digital signal processing and it has many important consequences which we will study in detail in the next chapters.
Finally, please note that, contrary to what happens in continuous time, in discrete time \textit{not all sinusoidal sequences are periodic.} In fact, only a very small subset of digital frequencies give rise to periodic patterns, namely only those of the form
\begin{equation}
\omega = \frac{M}{N}\pi, \quad M, N \in \mathbb{Z}.
\end{equation}
For all frequencies that are not a rational multiple of $\pi$, sinusoidal signals are two-sided, non-periodic infinite sequences. We will analyze the properties of the complex exponential in more detail again in Section~\ref{sec:fa:cexp}.
Elementary signal manipulations include, among others, the classic operations used in linear combinations; note that the following definitions apply without changes to all classes of signals.
\itempar{Scaling.}\index{scaling}
We can scale a signal $\mathbf{x}$ by a scalar factor $\alpha \in \mathbb{C}$ simply by multiplying each sample by $\alpha$:
\begin{equation}
\left(\alpha\mathbf{x} \right)[n] = \alpha x[n]
\end{equation}
If $\alpha$ is real, then the scaling represents a simple amplification (for $\alpha > 1$) or an attenuation (for $\alpha < 1$) of the signal. If $\alpha$ is complex, amplification and attenuation are compounded with a phase shift, whose meaning will be clearer in the next chapter.
\itempar{Sum.}\index{sum!of signals}
The sum of two sequences is their term-by-term sum:
\begin{equation}
\left( \mathbf{x + y} \right)[n] = x[n]+y[n]
\end{equation}
Sum and scaling are linear operators; informally, this means that they behave {}``intuitively''{}:
Operators act on the entire signal and transform it into another signal. We will study several important signal operators in the rest of this book and we introduce here some elementary examples.
\itempar{Delay.}\index{shift}
The delay operator delays a sequence by one ``time'' step:
When dealing with finite-length signals, we need to adjust the definition of the delay operator since the notation $x[n-k]$ is not defined for all values of $n$ and $k$. There are two ways to do so, and they correspond to applying the shift operator to either a periodic extension or a finite-support extension of the original length-$N$ signal $\mathbf{x}$:
\begin{itemize}
\item if we use a finite-support extension, the delay operator $\mathcal{D}^k$ performs a \textit{logical shift} on the elements of the signal vector, that is, the elements are shifted $k$ times to the left or to the right and zeros are inserted as a replacement, as illustrated here:
\item if we use a periodic extension, the delay operator $\mathcal{D}^k$ performs a \textit{circular shift} on the elements of the signal vector, that is, the elements are shifted $k$ times to the left or to the right and zeros are inserted as a replacement, as illustrated here:
As we said, the periodic extension is the natural extension of a finite-length signal so that, unless otherwise stated, in the rest of this book we will assume that all shifts applied to finite-length signals are circular shifts. Mathematically, we can formulate the circular shift via a modulo operation:
Note that, given our representation of signals as vectors, the delay operator can be encoded as a matrix-vector multiplication; for example, for a signal of length four, the circular shift can be represented as:
In general, the $N\times N$ delay matrix is a matrix with ones under the main diagonal, a one in the top right element and zeros everywhere else; see exercises~\ref{dt:ex:matrixop1} and~\ref{dt:ex:matrixop2} at the end of the chapter for more examples and properties. The matrix notation can be extended to infinite-length signals as well by using infinite matrices, although we will not pursue this line of analysis further in this book.
\itempar{Time reversal.}\index{shift}
The time reversal operator flips a sequence in time around the origin $n=0$:
\begin{equation}
\left( \mathcal{R} \mathbf{x} \right)[n] = x[-n]
\end{equation}
For finite-length signals, the effect of the time reversal operator is consistent with the periodic extension paradigm we just described:
In other words, when reversing a finite-length signal, we assume that the signal is in fact periodic and, when this is flipped around zero, the first sample remains in place and the others are in reverse order.
for finite-length signals, the lower limit in the sum is replaced by zero. It is easy to verify that the unit step can be obtained by integrating the impulse:
The discrete-time differentiation operator is the first-order difference:\footnote{
We will see in Chapter~\ref{is} that a more precise approximation to actual differentiation. For many applications, however, the first-order difference is enough.}
The reproducing formula is a simple application of the basic operations we just described; given a signal $\mathbf{x}$, and recalling that we denote a shifted delta as $\boldsymbol{\delta}_n = \mathcal{D}^n\boldsymbol{\delta}$, we can always write
where the sum is taken over all valid values for $n$. The formula states that any signal can be expressed as a \textit{linear combination of shifted impulses}. The weights are simply the signal values themselves and, while self-evident, this formula is important because it introduces the idea that a signal can be expressed as a linear combination of elementary building blocks. We will see in the next chapter that this is an instance of basis decomposition for a signal; while shifted impulses represent an intuitive canonical set of building blocks, by changing the basis set we can decompose a signal in many other ways, in order to highlight specific properties.
\subsection{Energy and Power}\label{ErgPowSec}
We define the \emph{energy}\index{energy!of a signal} of a discrete-time signal as
the squared-norm notation will be clearer after the next chapter but here it means that the sum is taken over all legal values for the index $n$, that is, over $\mathbb{Z}$ for infinite sequences and over the range $[0, N-1]$ for finite-length signals. This definition is consistent with the idea that, if the values of the sequence represent a time-varying voltage, the above sum would be proportional to the total energy (in joules) dissipated over a $1\,\Omega$~resistor. Obviously, the energy is finite only if the above sum converges; for finite-length signals this is always the case but for infinite sequences $\mathbf{x}$ must be \emph{square-summable}. A signal with this property is sometimes referred to as a \emph{finite- energy signal}. For a simple example of the converse, any periodic signal which is not identically zero is \emph{not} square-summable.
We define the \emph{power} \index{power!of a signal} of a infinite-length signal as the limit of the ratio of energy over time, taking the limit over the entire number of samples:
Clearly, signals whose energy is finite have zero total power since their energy dilutes to zero over an infinite time duration. Conversely, unstable sequences such as diverging exponential sequences (i.e.\ $e^{an}$ with $|a | > 1$) possess infinite power. Constant signals, on the other hand, whose energy is infinite, do have finite power; the same holds for periodic signals: although, formally, the limit in~(\ref{eq:dt:power}) is undefined for periodic sequences, we simply define their power to be their
\emph{average energy over a period}. Assuming that the period is $N$~samples, we have
At the beginning of the 19th century, French polymath Jean-Baptiste Joseph Fourier turned his attention to the way in which heat propagates inside a solid object. In the resulting memoir he introduced an analytical tool that was destined to become so successful, both in mathematics and in applied engineering, as to acquire the eponymous name of \textit{Fourier Transform}. Fourier's intuition was that many functions (and, remarkably, even discontinuous functions) could be represented as linear combinations of simple, harmonically-related sinusoidal components.
The systematization of the mathematical aspects of Fourier's theory took over a century to complete, and produced a remarkable body of work collectively known under the name of Harmonic Analysis. At the same time, the practical applicability of the Fourier transform gained immediate recognition across all scientific domains --- the sole drawback being the need to carry out cumbersome numerical calculation in order to work out the necessary coefficients. It is no surprise, then, that Fourier analysis enjoyed a renewed, explosive success as soon as electronic computation became a commodity: the well-known Fast Fourier Transform algorithm (interestingly, originally sketched by Gauss in~1805) is arguably the fundamental ingredient at the heart of today's personal communication devices.
To understand why Fourier Analysis plays such a central role in so many disciplines, and in signal processing in particular, let's consider the physical processes behind most of the interesting phenomena that we want to model or describe. Signals are time-varying quantities: you can imagine, for instance, the voltage level at the output of a microphone or the measured level of the tide at a particular location; in all cases, the variation of a signal over time implies that a transfer of energy is taking place somewhere. Now, a time-varying value which only increases monotonically over time is clearly a physical impossibility in the long run: fuses will blow, wires will melt and engines will explode. Oscillations, on the other hand, are nature's (and man's) way of keeping things in motion indefinitely without incurring a meltdown; from Maxwell's wave equation to the mechanics of the vocal cords, from the motion of an engine to the ebb and flow of the tide, oscillatory behavior is always the recurring theme. Sinusoidal oscillations are the purest form of such a constrained motion and, details aside, Fourier's lasting contribution was to show that one could express any reasonably well-behaved phenomenon as the combined output of a number of harmonically related sinusoidal sources.
In this chapter we will describe some key properties and results of Fourier analysis as applied to discrete-time signals. We have already mentioned in the previous chapter that, by using a vector space framework for signal processing, the Fourier transform can be described as a change of basis. This guiding principle will prove extremely useful as we navigate the subtle differences that exists between the different flavors of the transform and as we interpret their properties.
%Sinusoids have another remarkable property which justifies their ubiquitous presence. Indeed, \emph{any linear time-invariant transformation of a sinusoid is a sinusoid at the same frequency}: we express this by saying that sinusoidal oscillations are eigenfunctions of linear time-invariant systems. This is a formidable tool for the analysis and design of signal processing structures, as we will see in much detail in the context of discrete-time filters.
\section{Chapter Outline}
-The Fourier transform of a signal is an alternative representation of the data in the signal. While a signal lives in the \emph{time domain}\index{time domain},\footnote{\emph{Discrete}-time, of course.} its Fourier representation lives in the \emph{frequency domain}\index{frequency domain}. We can move back and forth at will from one domain to the other using the direct and inverse Fourier operators, since these operators are invertible.
+The Fourier transform of a signal is an alternative representation of the data in the signal. While a signal lives in the \emph{time domain}\index{time domain},\footnote{\emph{Discrete}\/ time, of course.} its Fourier representation lives in the \emph{frequency domain}\index{frequency domain}. We can move back and forth at will from one domain to the other using the direct and inverse Fourier operators, since these operators are invertible.
In this chapter we describe three types of Fourier transform which apply to the three main classes of signals that we have seen so far:
\begin{itemize}
-\item the Discrete Fourier Transform (DFT), which maps length-$N$ signals into a set of $N$ discrete frequency components; the transform is a change of basis in the underlying vector space $\mathbb{C}^N$.
-\item the Discrete Fourier Series (DFS), which maps $N$-periodic sequences into a set of $N$ discrete frequency components; this also is a change of basis in $\mathbb{C}^N$.
-\item the Discrete-Time Fourier Transform (DTFT), which maps infinite sequences into the space of $2\pi$-periodic function of a real-valued argument; this transform can also be interpreted structurally as a change of basis, as we will see in detail.
+\item the \textbf{Discrete Fourier Transform (DFT)}, which maps length-$N$ signals into a set of $N$ discrete frequency components; the transform is a change of basis in the underlying vector space $\mathbb{C}^N$.
+\item the \textbf{Discrete Fourier Series (DFS)}, which maps $N$-periodic sequences into a set of $N$ discrete frequency components; this also is a change of basis in $\mathbb{C}^N$.
+\item the \textbf{Discrete-Time Fourier Transform (DTFT)}, which maps infinite sequences into the space of $2\pi$-periodic function of a real-valued argument; this transform can also be interpreted structurally as a change of basis, as we will see in detail.
\end{itemize}
-The frequency representation of a signal (given by a set of coefficients in the case of the DFT and DFS and by a frequency distribution in the case of the DTFT) is called the \emph{spectrum}.
+The frequency representation of a signal (given by a set of coefficients in the case of the DFT and DFS and by a function the case of the DTFT) is called the \emph{spectrum}.
Regardless of the underlying signal space, the Fourier transform decomposes a discrete-time signal into a superposition of suitably scaled discrete-time oscillatory components. We will express the prototypical oscillation (i.e., the basic ingredient of all transforms) as a \textit{discrete-time complex exponential}, namely a sequence of the form\index{complex exponential}
\begin{equation}
x[n] = A\,e^{j(\omega n + \phi)} %= A \bigl[ \cos(\omega n + \phi) + j\sin(\omega n + \phi) \bigr]
\end{equation}
where $A \in \mathbb{R}$ is the amplitude, $\phi$ is the phase offset (in radians) and $\omega \in \mathbb{R}$ is the oscillation's frequency, also expressed in radians. Note that, although it is convenient to think of the index $n$ as a measure of ``time'', such time is a-dimensional and therefore the units for the frequency are simply radians (and not, say, radians per second).
As we have already explored in Section~\ref{sec:dt:frequency}, using simple algebra, we can easily see that
\begin{equation}
x[n+1] = x[n]\,e^{j\omega};
\end{equation}
that is, we can imagine the sequence $x[n]$ as the sequential output of a ``complex exponential generating machine'' that, at each step, takes the previous sample and multiplies it by $e^{j\omega}$, a pure phase factor\footnote{A phase factor is a complex number whose magnitude is equal to one. As such, it lives on the unit circle and it can be uniquely identified by its phase only.}. Recall now that, given a point $z$ on the complex plane, the position of $ze^{j\theta}$ can be found by rotating $z$ around the origin by an angle $\theta$ as shown in Figure~\ref{fig:fa:rotation}; the rotation is counterclockwise if $\theta$ is positive and clockwise if $\theta$ is negative. With this, $x[n]$ can be plotted on the complex plane as a series of point on a circle of radius $|A|$ centered on the origin; each point is at an angular distance of $\omega$ from the previous one. Two examples of complex exponential sequences are shown for a few values of $n$ in Figure~\ref{fig:fa:cexp_pn}, using positive and negative frequencies and different phase offsets. As you can see, the complex exponential perfectly captures the concept of a point rotating in circles, i.e. the most fundamental type of oscillatory movement.
The choice of a complex-valued signal as the prototypical oscillation may appear needlessly esoteric at first and, in fact, if we only consider real-valued signals we could derive the Fourier transform using only standard, real-valued trigonometric functions. There are however some major advantages in using complex exponentials:
\itempar{The math is simpler:} simply put, by using complex exponentials, trigonometry boils down to algebra. For instance, how many times have you struggled with the correct angle-sum formula? You remember that $\cos(\alpha + \beta)$ will be equal to a sum of products of sines and cosines but in what order and with what sign? Using complex algebra, however,
No need to remember formulas means fewer chances of making mistakes.
\itempar{The notation is more compact:} the prototypical oscillation encodes the movement of a rotating point and its description always possesses two degrees of freedom, magnitude and phase. We can choose to encode this information using the real-valued vertical and horizontal coordinates of the point on a Cartesian plane (that is, using a cosine and a sine); or we can interpret the point's position as a point on the complex plane using a complex exponential. While the two representations are equivalent, the latter is obviously more compact and, to see how cumbersome the notation becomes if we insist on using only real-valued quantities, please refer to Appendix~A, where some real-valued transforms are worked out.
The same holds for the \textit{values} of the Fourier transform; we will soon see that each Fourier coefficient encodes ``how much'' an oscillation of appropriate frequency contributes to the signal. While it may seem disorienting at first to have a complex-valued transform of a real-valued signal, note that each oscillatory contribution will be scaled in amplitude and offset in phase. Again, these two degrees of freedom can be conveniently packaged into a complex-valued Fourier coefficient!
\itempar{In the digital world, signals can be complex:} indeed, why not? Digital signals are just sequences of numbers that will be processed by general-purpose computational units, and therefore these sequences can certainly be complex-valued. While the interfaces to and from the physical world will necessarily handle real values only, \textit{internally} a DSP system will often be easier to design if complex-valued operations are allowed; this is particularly true in the case of communication systems. By starting off with complex exponentials as the prototypical oscillation, we are already equipped with a more versatile tool for frequency analysis.
\caption{Same point, many aliases: (a) adding multiples of $2\pi$ to a pure phase term does not change its position; (b) a positive (counterclockwise) displacement by $\theta$ is equivalent to a negative (clockwise) displacement by ($2\pi - \theta)$}\label{fig:fa:point_alias}
In discrete time, things are a bit different with respect to the properties of classic sinusoidal functions of a real variable. First of all, perhaps surprisingly, not all complex exponential sequences are periodic in $n$. Without loss of generality, consider a sequence with zero phase offset $x[n] = e^{j\omega n}$; for $x[n]$ to be periodic, there must exist an integer $N$ so that
\[
x[n] = x[n + N], \quad \forall n\in\mathbb{Z}.
\]
The above expression is equivalent to $e^{j\omega N} = 1$, that is, there must exist an integer $M$ such that
\[
\omega N = 2\pi \, M.
\]
Periodicity therefore requires the frequency to be of the form
or, in other words, in discrete time the only periodic oscillations are those \textit{whose frequency is a rational multiple of $2\pi$}.
Incidentally, \marginpar{LATER} note that the energy of the infinite-length signal $x[n] = A\,e^{j(\omega n + \phi)}$ is infinite while the power is equal to $|A|^2$ irrespective of frequency and phase.
Another peculiarity of discrete-time oscillations is that there is a limit on ``how fast'' we can go. To understand why, let's start by recalling that a pure phase term is always $2\pi$-periodic, in the sense that
This inherent phase ambiguity is called aliasing and stems from the simple fact that a point on the unit circle has an infinite number of possible ``names'', as shown in Figure~\ref{fig:fa:point_alias}-(a) --- etymologically, ``alias'' is Latin for ``otherwise''. Applied to discrete-time complex exponential sequences, this clearly implies an upper limit on the rotational speed of a point around the unit circle, since adding multiples of $2\pi$ to the value of the frequency will not change the values of the samples:
consequently, we can limit the range of distinct angular speeds to a representative interval of size $2\pi$. To choose the most suitable interval, consider what happens if we gradually increase the frequency of a discrete-time complex exponential starting from zero:
\begin{itemize}
\item for $0 \leq \omega < \pi$ we have a counterclockwise motion with increasing angular speed, i.e., we cover the full circle in fewer and fewer steps.
\item for $\omega = \pi$ we have the maximum possible forward speed of a discrete-time complex exponential; this corresponds to a sequence whose values are alternating between $+1$ and $-1$, which represents the maximum displacement attainable by successive points on the unit circle (antipodal points); we cover the full circle in 2 steps.
\item for $\pi < \omega < 2\pi$ at each step the point on the unit circle moves by more than $\pi$, as shown in Figure~\ref{fig:fa:cexp_alias}. Such a large counterclockwise motion is more ``economically'' explained by a \textit{clockwise} motion by an angle $2\pi- \omega < \pi$, i.e., the motion is better described by a \textit{negative} frequency whose magnitude is less than $\pi$.
\item for $\omega \geq 2\pi$ we can subtract a suitable multiple of $2\pi$ to $\omega$ until we fall into one of the three preceding cases.
\end{itemize}
The reference interval of choice for angular frequencies is therefore $[-\pi, \pi]$.
Note that, as we increase $\omega$ beyond $\pi$, we obtain a perceived reversal of direction and a \textit{decreasing} angular speed. This aliasing phenomenon is well known in cinematography where it is called the \textit{wagonwheel effect}; you can experience it in full by watching an old western movie: if a stagecoach enters the scene, you will see that its multi-spoked wheels seem to spin alternately backwards and forward as the speed of the vehicle changes. For a more detailed discussion of this optical illusion, see Appendix~B.
\caption{Complex-exponential sequence at angular speed $\omega = 2\pi - \theta$, with $\theta$ small: (a) at each step, the point's displacement is larger than $\pi$; (b) the movement is more ``economical'' if one assumes a negative frequency $\omega' = -\theta$.}\label{fig:fa:cexp_alias}
Consider $\mathbb{C}^N$, the space of complex-valued signals of finite length $N$; what are all the possible sinusoidal signals that span a \textit{whole} number of periods over the $[0, N-1]$ interval? We will presently show that:
\begin{itemize}
\item there are exactly $N$ such sinusoids
\item their frequencies are all harmonically related, i.e. they are all multiples of the fundamental frequency $2\pi/N$:
\begin{equation}\label{eq:fa:fund_freqs}
\omega_k = \frac{2\pi}{N}k, \quad k = 0, 1, \ldots, N-1;
\end{equation}
\item the set of $N$ length-$N$ complex exponentials at frequencies $\omega_k$ form a set of orthogonal vectors and therefore a basis for $\mathbb{C}^N$
\end{itemize}
With this basis, we are be able to express \textit{any} signal in $\mathbb{C}^N$ as a linear combination of $N$ harmonically-related sinusoids; the set of $N$ coefficients in the linear combination are called the \textit{Discrete Fourier Transform} of the signal, which can be easily computed algorithmically for any input data vector.
The fundamental oscillatory signal, the discrete-time complex-exponential $e^{j\omega n}$, is equal to $1$ for $n=0$; therefore, if the signal is to span a whole number of periods over $N$ points, we must have
\[
e^{j\omega N} = 1.
\]
In the complex field, the equation $z^N = 1$ has $N$ distinct solutions, given by the $N$ roots of unity
and so the $N$ possible frequencies that fulfill the $N$-periodicity requirements are those given in~(\ref{eq:fa:fund_freqs}). We can now use these frequencies to define a set $\bigl\{ \mathbf{w}_k \bigr\}_k$ containing $N$ signals of length $N$, where
\begin{equation}
w_k[n] = e^{j\frac{2\pi}{N}kn} \quad n, k = 0, 1, \dots, N-1.
\end{equation}
The real and imaginary parts of $\mathbf{w}_k$ for $N = 32$ and for some values of $k$ are plotted in Figures~\ref{fig:fa:basis_vector_0} to~\ref{fig:fa:basis_vector_31}; note how $\mathbf{w}_k = \mathbf{w}_{N-k}^*$.
\displaystyle \frac{1-e^{j\frac{2\pi}{N}(h-k)N}}{1-e^{j\frac{2\pi}{N}(h-k)}} = 0 & \mbox{ for } k \neq h
\end{cases}.
\end{align}
Because of orthogonality, as explained in section~\ref{sec:vs:bases}, the vectors form a basis for $\mathbb{C}^N$, called the \textit{Fourier basis} for the space of finite-length signals. In compact form, we can express the orthogonality of the Fourier vectors as:
Clearly the basis is not orthonormal; while it could be normalized by multiplying each vector by $1/\sqrt{N}$, in signal processing practice it is customary to keep the normalization factor explicit in the change of basis formulas. We too will follow this convention, that exists primarily for computational reasons.
In section~\ref{sec:vs:dt_space} we have illustrated how we can efficiently perform an orthonormal change of basis in $\mathbb{C}^N$; the Discrete Fourier Transform is such a transformation, allowing us to move from the time domain, represented by the canonical basis $\bigl\{ \boldsymbol{\delta}_k \bigr\}_k$, to the frequency domain, spanned by the Fourier basis $\bigl\{ \mathbf{w}_k \bigr\}_k$. Here, since the Fourier basis is orthogonal but not orthonormal, we simply need to slightly adjust the formulas in section~\ref{sec:vs:dt_space} to take into account the required normalization factors.
Given a vector $\mathbf{x} \in \mathbb{C}^N$ and the Fourier basis $\bigl\{ \mathbf{w}_k \bigr\}_k$, %the \textit{synthesis formula} allows us to
we can always express $\mathbf{x}$ as the following linear combination of basis vectors:
Using the orthogonality of the Fourier basis in~(\ref{eq:fa:orthogonalityCompact}), and taking the left inner product of the left- and right-hand sides of~(\ref{eq:fa:synthesis}) with each $\mathbf{w}_k$, we can see that the $N$~complex scalars $X[k]$, called the \textit{Fourier coefficients}, can be obtained simply as
\begin{equation}\label{eq:fa:fouriercoefficient}
X[k] = \langle \mathbf{w}_k, \mathbf{x} \rangle
\end{equation}
The coefficients capture the similarity between $\mathbf{x}$ and each of the basis vectors via an inner product; they are themselves a vector $\mathbf{X} \in \mathbb{C}^N$ and they can be computed all at once via the following matrix-vector multiplication, also known as the \textit{analysis formula}:
\[
\mathbf{X = Wx};
\]
the matrix $\mathbf{W}$ is built by stacking the Hermitian transposes of the basis vectors as
where, for convenience, we have introduced the scalar $W_N = e^{-j\frac{2\pi}{N}}$. Again using~(\ref{eq:fa:orthogonalityCompact}), it is immediate to show that $\mathbf{W}^H\mathbf{W} = \mathbf{I}\;N$, where $\mathbf{I}$ is the identity matrix. Since the change of basis is invertible, both $\mathbf{x}$ and $\mathbf{X}$ represent the same information, albeit from two different ``points of view'': $\mathbf{x}$ lives in the time domain, while $\mathbf{X}$ lives in the frequency domain.
In order to ``go back'' we can use the synthesis formula in matrix form, taking into account the explicit normalization factor:
\begin{equation} \label{eq:fa:idft_matrix}
\mathbf{x} = \frac{1}{N}\mathbf{W}^H\mathbf{X}.
\end{equation}
Some examples of Fourier matrices for low-dimensional spaces are:
\item the elements in the first row and the first columns are all equal to one since $W_N^0 = 1$ for all $N$;
\item powers of $W_N$ can be computed modulo $N$ because of the ``aliasing'' property of complex exponentials: $W^n_N = W_N^{n \mod N}$;
\item the matrices display a regular internal structure, which allows for extremely efficient computation methods such as the Fast Fourier Transform (FFT) algorithm (see Section~\ref{sec:fa:fft}).
\end{itemize}
\begin{comment}
In summary, the Discrete Fourier Transform is a change of basis in $\mathbb{C}^N$; from the canonical basis (the time domain) we move to the Fourier basis (the frequency domain). The elements of the vector in the frequency domain are obtained via the inner products
or, more compactly, via~(\ref{eq:fa:dft_matrix}). The elements $X_k$ are referred to as the \emph{spectrum} of the signal. The transform is perfectly invertible and, to move back to the time domain, we use the linear combination
The analysis and synthesis formula of the previous section can be written out explicitly in terms of the samples in the original vector and the vector of Fourier coefficients. This formulation highlights the algorithmic nature of the DFT and provides a straightforward way to implement the transform numerically.
The DFT coefficients can be computed using the following formula:
\begin{equation}\label{eq:fa:dft}
X[k] = \sum_{n=0}^{N-1}x[n]e^{-j\frac{2\pi}{N}nk}, \qquad \quad k = 0, \ldots, N-1
\end{equation}
while the inverse DFT is computed from the Fourier coefficients as
\begin{equation}\label{eq:fa:idft}
x[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}nk}, \qquad \quad n = 0, \ldots, N-1.
\end{equation}
The explicit formulas allows us to appreciate the highly structured form of the summations, which is exploited in the various efficient algorithmic implementations available in most numerical packages.
+The magnitude of the Fourier coefficients represent the energy distribution of the signal in the frequency domain, that is, ``how much'' of the signal is represented by an oscillation of frequency $(2\pi/N)k$. Parseval's theorem, as we saw in Section~\ref{sec:vs:parseval}, guarantees that energy is preserved across a orthogonal change of basis;
+%: $\|\mathbf{x}\|^2 = \sum|\alpha_k|^2$
+in the case of the DFT, this means that the norms of the time-domain and the frequency-domain coefficients satisfy the relationship
+The DFT, being a change of basis in $\mathbb{C}^N$, maps a length-$N$ vector $\mathbf{x}$ onto another length-$N$ vector $\mathbf{X}$. The graphical representation of $\mathbf{x}$ is intuitively clear as a succession of values in time; for the complex-valued DFT vector we choose a graphical representation that splits the information into \textit{magnitude} and \textit{phase}. The magnitude, as we have seen, represents the energy contribution of each basis vector while the phase represents the relative initial alignment of each oscillatory component necessary to recovers the time-domain signal, as we will show in detail in Section~\ref{sec:fa:reco}.
+The horizontal axis for the plot corresponds to the DFT coefficients' index $k$, which uniquely identifies the underlying frequency of the Fourier basis vectors, $\omega_k = (2\pi/N)k$; as we move along the horizontal axis from left to right, therefore, the coefficients correspond to increasingly faster frequencies up to a midpoint, after which the frequencies start to decrease again as shown in Figure~\ref{fig:fa:dftread}. The first half of the Fourier coefficients correspond to counterclockwise rotations, while the second half to clockwise rotations, as we saw in Section~\ref{sec:fa:cexp}.
+\itempar{Structure and symmetries.} The Fourier basis for $\mathbb{C}^N$ contains both clockwise and counterclockwise oscillations, corresponding to indices smaller or greater than $N/2$ respectively; since frequencies can always be taken modulo $2\pi$, the frequency corresponding to the index $k$ for $k > N/2$ is the same as a \textit{negative} frequency with index $k-N$:
+With this, we could always rearrange the plot of the Fourier coefficients around $k=0$, so that frequencies from zero to the fastest counterclockwise frequency lie on the right semiaxis while the corresponding clockwise frequencies lie to the left, as in Figure~\ref{fig:fa:dftsymm}; this would place slow-frequency coefficients together and push high frequencies to the edges of the plot.
+
+In the particular (and very common) case of a real-valued signal $\mathbf{x} \in \mathbb{R}^N$, it is easy to verify that the DFT is Hermitian-symmetric:
+this implies that the magnitude DFT of a real signal is symmetric, while its phase is antisymmetric. Because of the finite size of the data vector we must distinguish two cases:
+\begin{enumerate}
+ \item when $N$ is odd, $X[0]$ must be real-valued, and to each counterclockwise frequency corresponds a conjugate clockwise coefficient, as in the top row of Figure~\ref{fig:fa:dftsymmex}; the magnitude plot can be made symmetric around zero by shifting the last $(N-1)/2$ coefficients to the left semiaxis.
+ \item when $N$ is even, $X[0]$ and $X[N-1]$ are both real-valued, and every other intermediate frequency occurs in a counterclockwise-clockwise pair, as in the bottom row of Figure~\ref{fig:fa:dftsymmex}; the magnitude plot can be made symmetric around zero by shifting the last $N/2$ coefficients to the left semiaxis (thus duplicating the highest-frequency coefficient).
+\end{enumerate}
+Because of these symmetries, in practical DFT magnitude plots we only need to show the coefficients for $k=0,\ldots, \lfloor N/2 \rfloor$ since the remaining values are implicitly determined.
+ \caption{Magnitude DFT plots for odd-length (top) and even-length sequences (bottom); colors indicate equal-magnitude coefficients. The plots on the left show the standard DFT representation whereas those on the right show the coefficients rearranged to be symmetric around zero.}\label{fig:fa:dftsymmex}
+\end{figure}
+
+\clearpage
\subsection{DFT of Elementary Signals}
In the following examples we will compute the DFT of some elementary signals in $\mathbb{C}^{N}$, illustrating the examples for $N=64$; the associated plots will show the DFT coefficients in magnitude and phase.
\itempar{Impulse.} The DFT of the discrete-time delta is the constant signal $\mathbf{1}$ since
\[
\sum_{n=0}^{N-1}\delta[n]e^{-j\frac{2\pi}{N}nk} = e^{-j\frac{2\pi}{N}nk}\bigl|_{n=0} = 1 \quad \forall k.
\]
-For the shifted delta $\boldsymbol{\delta}_m = \mathcal{D}_m\{\boldsymbol{\delta}\}$, the $k$-th coefficient of the DFT is
+For the shifted delta $\boldsymbol{\delta}_m = \mathcal{D}_m\,\boldsymbol{\delta}$, the $k$-th coefficient of the DFT is
the DFT of an element of the canonical basis in time is the conjugate of the corresponding element of the Fourier basis. Note that, in the DFT of an impulse, all the coefficients have unit magnitude, that is, the most ``concentrated'' signal in time has nonzero Fourier coefficients at every frequency. This inverse relationship between time and frequency supports is a general property of Fourier analysis and will reappear frequently in the rest of the book.
\itempar{Rectangular signal.} Consider the step signal $\mathbf{x}$ defined by
\begin{equation}
x[n] = \begin{cases}
1 & \mbox{for $0 \leq n < M$} \\
0 & \mbox{for $M \leq n < N$,}
\end{cases} \label{eq:fa:step}
\end{equation}
shown in Figure~\ref{fig:fa:ex4} for $M=5$ and $N=64$. We can express the signal as
In the derivation above, we have manipulated the expression for $X[k]$ into a product of a real-valued term (which captures the magnitude) and a pure phase term; this allows us to easily plot the DFT as in Figure~\ref{fig:fa:ex4}. Note that, while the phase grows linearly with $k$, it is customary in phase plots to ``wrap'' its value over a $2\pi$-wide interval; in signal processing the interval of reference is $[-\pi, \pi]$.
\itempar{Constant signal.} By setting $M=N$ in~(\ref{eq:fa:stepDFT}) we obtain the DFT of the constant signal $\mathbf{x=1}$; the sum, as in~(\ref{eq:fa:orthogonality}), shows the orthogonality of the roots of unity and we have
\begin{equation} \label{eq:fa:unitDFT1}
X[k] = \sum_{n=0}^{M-1}e^{-j\frac{2\pi}{N}nk}
= \begin{cases}
N & \mbox{ for } k=0 \\
0 & \mbox{ otherwise}
\end{cases}
\end{equation}
so that
\begin{equation} \label{eq:fa:unitDFT2}
\DFT{\boldsymbol{1}} = N\boldsymbol{\delta}.
\end{equation}
The above result, up to a normalization factor, is the dual of what we obtained when we computed the DFT of the delta signal; in this case the time-domain signal with the largest support yields a transform that has a single nonzero coefficient.
\itempar{Harmonic sinusoids.} The fundamental frequency for $\mathbb{C}^{N}$ is $2\pi/N$ and a complex exponential at a multiple of this frequency will coincide with a Fourier basis vector. Given the signal $\mathbf{x}$ defined by
\[
x[n] = e^{j\frac{2\pi}{N}mn}
\]
it is clearly $\mathbf{x} = \mathbf{w}_m$ and its DFT coefficients can be easily computed using~(\ref{eq:fa:orthogonalityCompact}) as
up to a normalization factor, this is the dual of~(\ref{eq:fa:DFTdelta}). Note that the result in the previous section, for the constant signal $\mathbf{x=1}$, is just a particular case of the above relationship when $m=0$.
We can now easily compute the DFT of standard trigonometric functions in which the frequency is harmonically related to the fundamental frequency for the space. Consider for instance the signal $\mathbf{x}$ defined by
\[
x[n] = \cos \left(\frac{\pi}{8} n \right), \qquad n=0,1,2,\ldots, 63.
\]
With a simple manipulation we can write:
\begin{align*}
\cos \left(4\frac{2\pi}{64} n \right) &= \frac{1}{2}\left( e^{j4\frac{2\pi}{64}n} + e^{-j4\frac{2\pi}{64}n}\right) \\
where we have used the fact that we can take all frequency indexes modulo 64 because of the aliasing property of complex exponentials; we can therefore express $\mathbf{x}$ as
Since the signal is the sum of two Fourier basis vectors, orthogonality implies that only two of the inner products in~(\ref{eq:fa:fouriercoefficient}) will be nonzero:
The DFT of the signal is plotted in Figure~\ref{fig:fa:ex1}; the spectrum shows how the entire frequency content of the signal is concentrated over two single frequencies. Since the original signal is real-valued, the DFT component at $k=60$ ensures that the imaginary parts in the reconstruction formula cancel out; this symmetry is a general property of the Fourier transform that we will examine in more detail later.
Consider now a slight variation of the previous signal obtained by introducing a phase offset:
so that the resulting DFT coefficients are all zero except for
\begin{align*}
X[4] &= 32\,e^{j2\pi/3} \\
X[60] &= 32\,e^{-j2\pi/3}.
\end{align*}
The resulting DFT is plotted in Figure~\ref{fig:fa:ex2}; the magnitude does not change but the phase offset is reflected in nonzero values for phase at $k=4, 60$.
Consider now a sinusoid whose frequency is {\em not} a multiple of the fundamental frequency for the space, such as
\[
x[n] = \cos \left(\frac{\pi}{5} n \right), \qquad n=0,1,2,\ldots, 63.
\]
In this case we cannot decompose the signal into a sum of basis vectors and we must therefore explicitly compute all the DFT coefficients. We could do this algebraically and work out the resulting geometric sums as we did for the step signal. More conveniently, however, since the DFT is at heart a \textit{numerical} algorithm, we can just use a standard numerical package (Numpy, Matlab, Octave) and use the built-in \texttt{fft()} function. The resulting DFT is shown in Figure~\ref{fig:fa:ex3}; note how here \textit{all} the DFT coefficients are nonzero. While the magnitude is larger for frequencies close to that of the original signal ($6\pi/64 < \pi/5 < 7\pi/64$), to reconstruct $\mathbf{x}$ exactly, we need a contribution from each one of the basis vectors.
+\subsection{The DFT as a Synthesis Tool}\label{sec:fa:reco}
An instructive way to explain the physical interpretation of the DFT is to start with the synthesis formula~(\ref{eq:fa:idft}). The expression tells us that we can build \textit{any} signal in $\mathbb{C}^N$ as a suitable combination of $N$ complex sinusoids with harmonically related frequencies; the magnitude and initial phase of each oscillation are given by the complex-valued Fourier coefficients. Consider a ``complex exponential generator'' as in Figure~\ref{fig:fa:gen}; the oscillator works at a frequency $(2\pi/N)k$ and it can be tuned in magnitude and phase. The DFT synthesis formula can be represented graphically as in Figure~\ref{fig:fa:idft}, that is, as a bank of $N$ oscillators working in parallel; to reproduce any signal $\mathbf{x}$ from its DFT $\mathbf{X}$:
\begin{itemize}
\item set the amplitude $A_k$ of the $k$-th generator to $\bigl|X[k]\bigr|$, i.e.\ to the magnitude of the $k$-th DFT coefficient;
\item set the phase $\phi_{k}$ of the $k$-th generator to $\measuredangle{X[k]}$, i.e.\ to the phase of the $k$-th DFT coefficient;
\item start all the generators at the same time and sum their outputs for $N$ cycles
\item divide the result by $N$.
\end{itemize}
This ``machine'' shows that each Fourier coefficient captures ``how much of'' and ``how in phase'' an oscillation at frequency $2\pi/k$ is contained in $\mathbf{x}$; this is consistent with the fact that each $X[k]$ is computed as the inner product between $\mathbf{x}$ and $\mathbf{w}_k$, and that the inner product is a measure of similarity.
consequently, the square magnitude of a DFT coefficient $|X[k]|^2$ is proportional, via to a scale factor $N$, to the energy of $\mathbf{x}$ at the frequency $(2\pi/N)k$: the magnitude of the DFT therefore shows how the global energy of the original signal is distributed in the frequency domain. The phase of each DFT coefficient specifies the initial phase of each oscillator in the reconstruction formula, i.e. the \emph{relative alignment} of each complex exponential at the onset of the signal. While this does not affect the energy distribution in frequency, it does have a significant impact on the \textit{shape} of the signal in the time domain as we will see in more detail shortly.
On a historical note, far from being a pure academic exercise, the synthesis of signals via a bank of oscillators is a topic that attracted great interest and efforts long before the advent of digital computers. Figure~\ref{fig:fa:tidemachine} shows a mechanical implementation of the inverse DFT as we just described; in this case the device was designed in order to anticipate the evolution of the sea tide. The original design of these tide prediction machines dates back to Lord Kelvin in the 1860s.
-The DFT, being a change of basis in $\mathbb{C}^N$, maps a length-$N$ vector $\mathbf{x}$ onto another length-$N$ vector $\mathbf{X}$. The graphical representation of $\mathbf{x}$ is intuitively clear as a succession of values in time; for the DFT vector, as we said, we choose a graphical representation that splits the information into magnitude and phase. The (discrete) horizontal axis for the plot corresponds to the DFT coefficients' index $k$ which uniquely identifies the underlying frequency of the Fourier vector $\omega_k = (2\pi/N)k$; as we move along the horizontal axis from left to right, therefore, we will encounter coefficients that correspond to increasingly faster frequencies up to the midpoint $\lfloor k/2 \rfloor$, after which the frequencies will start to decrease again as shown in Figure~\ref{fig:fa:dftread}. The first half of the Fourier coefficients correspond to counterclockwise rotations, while the second half to clockwise rotations, as we saw in Section~\ref{sec:fa:cexp}.
-
-\itempar{Symmetries.} When the original signal $\mathbf{x}$ is real-valued, the DFT is hermitian-symmetric:
-\[
- X[k] = X^*[N-k];
-\]
-this implies that the magnitude DFT of a real signal is symmetric.
-
-
\itempar{Labeling the frequency axis.}
Each DFT coefficient is the inner product between the ``input'' signal and a sinusoidal signal with frequency $\omega_k = 2\pi k /N$, that is, each DFT coefficients indicates how much of the input signal is
Back to the theoretical side, consider again the DFT reconstruction formula in~(\ref{eq:fa:idft}) (or, equivalently, the ``machine'' in Figure~\ref{fig:fa:gen}); normally, the expression is defined for $0 \leq n < N$ but, if the index $n$ is outside of this interval, we can always write $n = mN + i$ with $m \in \mathbb{Z}$ and $i= n \mod N$. With this,
In other words, due to the aliasing property of the complex exponential, the inverse DFT formula generates an infinite, \emph{periodic} sequence of period $N$; this should not come as a surprise, given the $N$-periodic nature of the basis vectors for $\mathbb{C}^N$. Similarly, the DFT analysis formula remains valid if the frequency index $k$ is allowed to take values outside the $[0, N-1]$ interval and the resulting sequence of DFT coefficients is also an $N$-periodic sequence.
The natural Fourier representation of periodic signals is called the Discrete Fourier Series (DFS) and its explicit analysis and synthesis formulas are identical to~(\ref{eq:fa:dft}) and~(\ref{eq:fa:idft}), modified only with respect to the range of the indexes, which now span $\mathbb{Z}$; the DFS represents a change of basis in the space of periodic sequences $\tilde{\mathbb{C}}^N$. Since there is no mathematical difference between the DFT and the DFS, it is important to remember that even in the space of finite-length signals everything is implicitly $N$-periodic.
\subsection{Circular shifts revisited}
In Section~\ref{sec:dt:operators} we stated that circular shifts are the ``natural'' way to interpret how the delay operator applies to finite-length signals; considering the inherent periodicity of the DFT/DFS, the reason should now be clear. Indeed, the delay operator is always well-defined for a periodic signal $\mathbf{\tilde{x}}$ and, given its DFS $\mathbf{\tilde{X}}$, the $k$-th DFS coefficient of the sequence shifted by $m$ is easily computed as
in other words, a delay by $m$ samples in the time domain becomes a linear phase shift by $-2\pi m/N$ in the frequency domain.
With a finite-length signal $\mathbf{x}$, for which time shifts are not well defined, we can still always compute the DFT, multiply the DFT coefficients by a linear phase shift and compute the inverse DFT. The result is always well defined and, by invoking the mathematical equivalence between DFT and DFS, it is straightforward to show that
which justifies the circular interpretation for shifts of finite-length signals.
\subsection{DFT of multiple periods}
All the information carried by a $N$-periodic discrete-time signal is contained in $N$ consecutive samples and therefore its complete frequency representation is provided by the DFS, which coincides with the DFT of one period. This intuitive fact is confirmed if we try to compute the DFT of $L$ consecutive periods:
The above results shows that the DFT of $L$ periods is obtained simply by multiplying the DFT coefficients of one period by $L$ and appending $L-1$ zeros after each one of them.
\subsection{Pushing the DFS to the limit}
In the next section we will derive a complete frequency representation for aperiodic, infinite-length signals, which is still missing; but right now, the DFS can help us gain some initial intuition if we imagine such signals as the limit of periodic sequences when the period grows to infinity. Consider an aperiodic, infinite-length and absolutely summable sequence $\mathbf{x}$; given any integer $N > 0$ we can always build an $N$-periodic sequence $\mathbf{\tilde{x}}_N$, with \index{periodization}
the convergence of the sum for all $n$ is guaranteed by the absolute summability of $\mathbf{x}$ (see also Example~\ref{exa:dt:periodization}). As $N$ grows larger, the copies in the periodization will be spaced more and more apart and, in the limit,
in the above, we have used the definition of $\mathbf{\tilde{x}}_N$ and exploited the fact that $e^{-j(2\pi /N)nk} = e^{-j(2\pi /N)(n+pN)k}$. Now, for every value of $p$ in the outer sum, the argument of the inner sum varies between $pN$ and $pN + N - 1$ so that the double sum can be simplified to
it is immediate to see that, for every value of the period $N$, the DFS coefficients of $\mathbf{\tilde{x}}_N$ are given by regularly spaced samples of $X(\omega)$ computed at multiples of $2\pi/N$:
Figure~\ref{fig:fa:dsf2dtft} shows some examples for different values of $N$. As $N$ grows large, the set of samples will grow denser in the $[0, 2\pi]$ interval; since, in the limit, $\mathbf{\tilde{x}}_N$ tends to $\mathbf{x}$, it appears that $X(\omega)$ would indeed be a suitable frequency-domain representation for $\mathbf{x}$, as we will show momentarily.
\caption{top row: original infinite-length, absolutely summable signal (left) and the function $X(\omega)$ defined in~(\ref{eq:fa:dtftAnteLitteram}); rows 2-5: periodized signal $\mathbf{\tilde{x}}_N$ and its DFS for increasing values of $N$; all DFS values coincide with samples of $X(\omega)$.}\label{fig:fa:dsf2dtft}
Comment on the result: you should point out two major problems.
\end{enumerate}
As it appears, interpolating with a zero-order hold introduces in-band distortion in the region $|f| < 1/2$ and out-of-band spurious components at higher frequencies. Both problems could however be fixed by a well-designed continuous-time filter $G(f)$ applied to the ZOH interpolation.
\begin{enumerate}
\setcounter{enumi}{2}
\item Sketch the frequency response $G(f)$
\item Propose two solutions (one in the continuous-time omain, and another in the discrete-time domain) to eliminate or attenuate the in-band distortion due to the zero-order hold. Discuss the advantages and disadvantages of each.
i.e. a signal obtained from the discrete-time sequence using a zero-centered zero-order hold with interpolation period $T_s = 1$s. Let $X_0(f)$ be the Fourier transform of $x_0(t)$.
\begin{enumerate}
\item Express $X_0(f)$ in terms of $X(e^{j\omega})$.
\item Compare $X_0(f)$ to $X(f)$, where $X(f)$ is the spectrum of the continuous-time signal obtained using an ideal sinc interpolator with $T_s=1$:
Comment on the result: you should point out two major problems.
\item The signal $x(t)$ can be obtained back from the zero-order hold interpolation via a continuous-time filtering operation:
\[
x(t)=x_0(t)\ast g(t).
\]
Sketch the frequency response of the filter $g(t)$.
\item Propose two solutions (one in the continuous-time domain, and another in the discrete-time domain) to eliminate or attenuate the distortion due to the zero-order hold. Discuss the advantages and disadvantages of each.
\end{enumerate}
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
\begin{enumerate}
\item
\begin{align*}
X_0(f) &= \int_{-\infty}^{\infty} x_0(t)e^{-j2\pi f t} dt\\
&= \int_{-\infty}^{\infty} \sum_{n=-\infty}^{\infty}x[n]\mbox{rect}\,(t-n)e^{-j2 \pi f t} dt \\
&= \sum_{n=-\infty}^{\infty} x[n] \int_{-\infty}^{\infty} \mbox{rect}\,(t-n)e^{-j2\pi f t} dt \\
&= \sum_{n=-\infty}^{\infty}x[n] e^{-j2\pi f n} \, \int_{-1/2}^{1/2} e^{-j2\pi f \tau} d\tau \\
\item To understand the effects of the zero-order hold consider for instance a discrete-time signal with a triangular
spectrum, as shown in the left panel below. We know that sinc interpolation will exactly preserve the shape of the spectrum and return a continuous-time signal that is strictly bandlimited to the $[-F_s/2, F_s/2]$ interval (with $F_s = 1/T_s = 1$), that is:
Conversely, the spectrum of the continuous-time signal interpolated by the zero-order hold is the product of $X(e^{j2\pi f})$, which is periodic with period $F_s = 1$ Hz, and the term $\sinc(f)$, whose first spectral null is for $f=1$ Hz. Here are the two terms, and their product, in magnitude:
- \dspFunc{x \dspPeriodize \dspTri{0}{1} x \dspSinc{0}{2} mul abs}
+ \dspFunc[linecolor=ColorCF]{x \dspPeriodize \dspTri{0}{1} x \dspSinc{0}{2} mul abs}
\dspCustomTicks[axis=x]{0 0 2 1 4 2}
\end{dspPlot}
\end{center}
There are two main problems in the zero-order hold interpolation as compared to the sinc interpolation:
\begin{itemize}
\item The zero-order hold interpolation is NOT bandlimited: the $2\pi$-periodic replicas of the digital spectrum leak through in
the continuous-time signal as high frequency components. This is due to the sidelobes of the interpolation function in the frequency domain (rect in time $\leftrightarrow$ sinc in frequency) and it represents an undesirable high-frequency content which is typical of all local interpolation schemes.
\item There is a distortion in the main portion of the spectrum (that between$-F_s/2$ and $F_s/2 = 0.5$ Hz) due to the non-flat frequency response of the interpolation function. It can be seen if we zoom in the baseband portion:
\item A first solution is to compensate for the distortion introduced by $G(f)$ in the discrete-time domain. This is
equivalent to pre-filtering $x[n]$ with a discrete-time filter of magnitude $1/G(e^{j2\pi f})$ . The advantages of this method is that digital filters such as this one are relatively easy to design and that the filtering can be done in the discrete-time domain. The disadvantage is that this approach does not eliminate or attenuate the high frequency leakage outside the baseband.
In continuous time, one could cascade the interpolator with an analog lowpass filter to eliminate the leakage. The disadvantage is that it is hard to implement an analog lowpass which can also compensate for the in-band distortion introduced by $G(f)$; such a filter will also introduce unavoidable phase distortion (no analog filter has linear phase).
Consider the local interpolation scheme of the previous exercise but assume that the characteristic of the
interpolator is the following:
\[
I(t) = \begin{cases}
1 - 2|t| & \mbox{ for } |t| \leq 1/2 \\
0 & \mbox{ otherwise}
\end{cases}
\]
This is a triangular characteristic with the same support as the zero-order hold. If we pick an interpolation interval $T_s$ and interpolate a discrete-time signal $x[n]$ with $I(t)$, we obtain the continuous-time signal:
Assume that the spectrum of $x[n]$ between $-\pi$ and $\pi$ is
\[
X(e^{j\omega}) =
\begin{cases}
1 & \mbox{ for } |\omega| \leq 2\pi/3 \\
0 & \mbox{ otherwise}
\end{cases}
\]
(with the obvious $2\pi$-periodicity over the entire frequency axis).
\begin{enumerate}
\item Compute and sketch the Fourier transform $I(f)$ of the interpolating function $I(t)$. (Recall that the triangular function can be expressed as the convolution of $\mbox{rect}(t/2)$ with itself).
\item Sketch the Fourier transform $X(f)$ of the interpolated signal $x(t)$; in particular, clearly mark the Nyquist frequency $f_N = 1/(2T_s)$.
\item The use of $I(t)$ instead of a sinc interpolator introduces two types of errors: briefly describe them.
\item To eliminate the error \emph{in the baseband} $[-f_N, f_N]$ we can pre-filter the signal $x[n]$ with a filter $h[n]$ \emph{before} interpolating with
$I(t)$. Write the frequency response of the discrete-time filter $H(e^{j\omega})$.
An alternative way of describing the sampling operation relies on the concept of \textit{modulation by a pulse train}. Given a sampling interval $T_s$, a continuous-time pulse train $p(t)$ is an infinite collection of equally spaced Dirac deltas:
\[
p(t) = \sum_{k=-\infty}^{\infty}\delta(t-kT_s).
\]
The pulse train is the used to modulate a continuous-time signal:
\[
x_s(t) = p(t)\,x(t).
\]
Intuitively, $x_s(t)$ represents a ``hybrid'' signal where the nonzero values are only those of the discrete time samples that would be obtained by raw-sampling $x(t)$ with period $T_s$; however, instead of representing the samples a countable sequence (i.e. with a different mathematical object) we are still using a continuous-time signal that is nonzero only over infinitesimally short instants centered on the sampling times. Using Dirac deltas allows us to embed the instantaneous sampling values in the signal.
Note that the Fourier Transform of the pulse train is
(where, as per usual, $F_s = 1/T_s$). This result is a bit tricky to show, but the intuition is that a periodic set of pulses in time produces a periodic set of pulses in frequency and that the spacing between pulses frequency is inversely proportional to the spacing between pulses in time.
Derive the Fourier transform of $x_s(t)$ and show that if $x(t)$ is bandlimited to $F_s/2$, where $F_s = 1/T_s$, then we can reconstruct $x(t)$ from $x_s(t)$.
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
By using the modulation theorem, the product in time becomes a convolution in frequency:
In other words, the spectrum of the delta-modulated signal is the periodization (with period $F_s=1/T_s$) of the original spectrum. If the latter is bandlimited to $F_s/2$ there will be no overlap between copies in the periodization and therefore $x(t)$ can be obtained simply by lowpass filtering $x_s(t)$ in the continuous-time domain.
\item What is the bandwidth of the signal? What is the minimum sampling period in order to satisfy the sampling theorem?
\item Take a sampling period $T_s = 1/(2f_0)$; clearly, with this sampling period, there will be aliasing. Plot the DTFT of the discrete-time signal $x_a[n] = x_c(nT_s)$.
\item Suggest a block diagram to reconstruct $x_c(t)$ from $x_a[n]$.
\item With such a scheme available, we can therefore exploit aliasing to reduce the sampling frequency necessary to sample a bandpass signal. In general, what is the minimum sampling frequency to be able to reconstruct, with the above strategy, a real-valued signal whose frequency support on the positive axis is $[f_0, f_1]$ (with the usual symmetry around zero, of course)?
\end{enumerate}
\end{exercise}
}\fi
\ifanswers{%
\begin{solution}{}
\begin{enumerate}
\item The highest nonzero frequency is $2f_0$ and therefore $x_c(t)$ is $4f_0$-bandlimited. The minimum sampling frequency that satisfies the sampling theorem is $F_s = 4f_0$. Note however that the support over which the (positive) spectrum is nonzero is the interval $[f_0, 2f_0]$ so that one could say that the total \emph{effective} bandwidth of the signal is only $2f_0$.
\item The digital spectrum will be the periodized continuous-time spectrum, rescaled to $[-\pi, \pi]$; the periodization after sampling at a frequency $F_a = 2f_0$, yields
for instance, as $\omega$ goes from zero to $\pi$, the nonzero contribution to the DTFT will be the term $X_c(\frac{\omega}{\pi}f_0 - 2f_0)$ where the argument goes from $-2f_0$ to $-f_0$. The spectrum is represented here with periodicity explicit:
\item filter in continuous time $x_p(t)$ with an ideal bandpass filter with (positive) passband equal to $[f_0, 2f_0]$ to obtain $x_c(t)$.
\end{itemize}
\item The effective \emph{positive} bandwidth of a signal whose spectrum is nonzero only over $[-f_1, -f_0] \cup [f_0, f_1]$ is $W = f_1 - f_0$. Obviously the sampling frequency must be at least equal to the \textit{total} effective bandwidth, so a first condition on the sampling frequency is $F_s \geq 2W.$ We can now distinguish two cases.
\begin{itemize}
\item[1)] assume $f_1$ is a multiple of the positive bandwidth, i.e.\ $f_1 = MW$ for some integer $M > 1$ (for $x[n]$ above, it was $M = 2$). Then the argument we made before can be easily generalized: if we pick $F_s = 2W$ and sample we have that
Since $f_0 = f_1 - W = (M-1)W$, this translates to
\begin{align*}
(2k+M-1)W & \leq f \leq (2k+M)W \\
(2k-M)W & \leq f \leq (2k-M+1)W
\end{align*}
which, once again, are non-overlapping intervals.
\item[2)] if $f_1$ is \emph{not} a multiple of $W$ the easiest thing to do is to decrease the lower frequency $f_0$ to a new {\em smaller} frequency $f_0'$ so that the new positive bandwidth $W' = f_1 - f_0'$ divides $f_1$ exactly. In other words we set a new lower frequency $f_0'$ so that it will be $f_1 = M(f_1-f_0')$ for some integer $M>1$; it is easy to see that
\[
M = \biggl\lfloor \frac{f_1}{f_1 - f_0} \biggr\rfloor.
\]
since this is the maximum number of copies of the $W$-wide spectrum which fit \emph{with no overlap} in the $[0, f_0]$ interval. If $W > f_0$ obviously we cannot hope to reduce the sampling frequency and we have to use normal sampling. This artificial change of frequency will leave a small empty ``gap'' in the new bandwidth $[f_0', f_1]$, but that's no problem. Now we use the previous result and sample at $F_s = 2(f_1 - f_0')$ with no overlaps. Since $f_1 - f_0' = f_1/M$, we have that, in conclusion, the minimum sampling frequency is
\begin{exercise}{Digital processing of continuous-time signals}
For your birthday, you receive an unexpected present: a $4$~MHz sampler, complete with anti-aliasing filter. This means you can safely sample signals up to a frequency of $2$ ~MHz; since this frequency is above the AM radio frequency band, you decide to hook up the sampler to your favorite signal processing system and build an entirely digital radio receiver. In this exercise we will explore how to do so.
To simplify things a bit, assume that the AM radio spectrum extends from $1$~Mhz to $1.2$~Mhz and that in this band you have ten channels side by side, each one of which occupies $20$~KHz.
\begin{enumerate}
\item Sketch the digital spectrum at the output of the A/D converter, and show the bands occupied by the channels, numbered from $1$ to $10$, with their beginning and end frequencies.
\end{enumerate}
The first thing that you need to do is to find a way to isolate
the channel you want to listen to and to eliminate the rest.
For this, you need a bandpass filter centered on the band of
interest. Of course, this filter must be {\em tunable} in
the sense that you must be able to change its spectral
location when you want to change station. An easy way to
obtain a tunable bandpass filter is by modulating a lowpass
filter with a sinusoidal oscillator whose frequency is
controllable by the user:
\begin{enumerate}\setcounter{enumi}{1}
\item As an example of a tunable filter, assume $h[n]$ is an ideal lowpass filter with cutoff frequency $\pi/8$. Plot the magnitude response of the filter $h_m[n] = \cos(\omega_m n)h[n]$, where $\omega_m = \pi/2$; $\omega_m$ is called the {\em tuning frequency}.
\item Specify the cutoff frequency of a lowpass filter which can be used to select one of the AM channels above.
\item Specify the tuning frequencies for channel~$1$, $5$ and $10$.
\end{enumerate}
Now that you know how to select a channel, all that is left to do is to demodulate the signal and feed it to an interpolator and then to a loudspeaker.
\begin{enumerate}\setcounter{enumi}{4}
\item Sketch the complete block diagram of the radio receiver, from the antenna going into the sampler to the final loudspeaker. Use only one sinusoidal oscillator. Do not forget the filter before the interpolator (specify its bandwidth).
\end{enumerate}
The whole receiver now works at a rate of $4$ MHz; since it outputs audio signals, this is clearly a waste.
\begin{enumerate}\setcounter{enumi}{5}
\item Which is the minimum interpolation frequency you can use? Modify the receiver's block diagram with the necessary elements to use a lower frequency interpolator.
Assume $x(t)$ is a continuous-time pure sinusoid at $10$ KHz. It is sampled with a sampler at $8$ KHz and then interpolated back to a continuous-time signal with an interpolator at $8$ KHz. What is the perceived frequency of the interpolated sinusoid?
We have seen that any discrete-time sequence can be sinc-interpolated into a continuous-time signal which is $F_s$-bandlimited; $F_s$ depends on the
interpolation interval $T_s$ via the relation $F_s = 1/T_s$.
Consider the continuous-time signal $x_c(t) = e^{-t/T_s}u(t)$ and the discrete-time sequence $x[n] = e^{-n}u[n]$. Clearly, $x_c(nT_s) = x[n]$; but, can we also say that $x_c(t)$ is the signal we obtain if we apply sinc interpolation to the sequence $x[n] = e^{-n}$ with interpolation interval $T_s$?
Consider a real continuous-time signal $x(t)$. All you know about the signal is that $x(t) = 0$ for $|t| > t_0$. Can you determine a sampling frequency $F_s$ so that when you sample $x(t)$, there is no aliasing? Explain.