Homec4science
Diffusion epfl (master) Browse

R2653/writing/sp4comm.multipub/110-multirate/10-mr-multirate.tex
a3c8ac432673master

Last Change
Raw File

10-mr-multirate.tex
View Options

\chapter{Multirate Signal Processing}
\label{ch:mr}
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Upsampling}
\label{sec:mr:up}\index{upsampling|(}
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
\begin{equation} \label{eq:mr:up}
(\mathcal{U}_N \mathbf{x})[n] = x_{N\uparrow}[n] = \
\begin{cases}
x[k] & \mbox{ for $n = kN,\;\; k \in \mathbb{Z}$} \\
0 & \mbox{ otherwise}
\end{cases} % \right.
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b]
\center
\begin{dspBlocks}{1.3}{0.6}
$x[n]$ & \BDupsmp{N} & $x_{N\uparrow}[n]$
\psset{arrows=->,linewidth=1pt}
\ncline{1,1}{1,2} \ncline{1,2}{1,3}
\end{dspBlocks}
\caption{Symbol for the upsampling operator.}\label{fig:mr:downSym}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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\:
\begin{align}
X_{N\uparrow}(z) &= \sum_{n=-\infty}^{\infty} x_{N\uparrow}[n]z^{-n} \nonumber \\
&= \sum_{k=-\infty}^{\infty} x[k]z^{-kN} = X(z^N)
\end{align}
and therefore
\begin{equation} %\label{}
X_{N\uparrow}(e^{j\omega}) = X(e^{j\omega N}).
\end{equation}
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
%
\def\dtSig{5 sub 0.3 mul RadtoDeg %
dup cos 0.3 mul exch %
0.5 mul sin 0.6 mul %
add 0.5 add }
%
\center
\begin{dspPlot}[height=\dspHeightCol,xout=true,sidegap=1]{-8,8}{-0.5,1.1}
\dspSignal[linecolor=ColorDT]{x \dtSig}
\end{dspPlot}
\begin{dspPlot}[height=\dspHeightCol,xout=true,sidegap=4]{-32,32}{-0.5,1.1}
\dspSignal[linecolor=ColorDT!65]{x dup cvi 4 mod 0 eq {4 div \dtSig} {pop 0} ifelse}
\dspSignal[linecolor=ColorDT,xunit=4,xmax=8,xmin=-8]{x \dtSig}
\end{dspPlot}
\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}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{1}}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-5,5}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{1}}
\pnode(-2,0){leftA}\pnode( 2,0){rightA}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X_{2\downarrow}(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x 2 mul \dspPeriodize \dspTri{0}{1}}
\pnode(-1,1.2){leftB}\pnode( 1,1.2){rightB}
\ncline[linewidth=2pt,linecolor=blue!40]{leftA}{rightA}
\ncline[linewidth=1pt,linecolor=blue!40,linestyle=dashed]{->}{leftA}{leftB}
\ncline[linewidth=1pt,linecolor=blue!40,linestyle=dashed]{->}{rightA}{rightB}
\end{dspPlot}
\caption{Upsampling by $2$.}\label{fig:mr:upA}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!htb]
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{1}}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-5,5}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{1}}
\pnode(-3,0){leftA}\pnode( 3,0){rightA}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X_{2\downarrow}(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x 3 mul \dspPeriodize \dspTri{0}{1}}
\pnode(-1,1.2){leftB}\pnode( 1,1.2){rightB}
\ncline[linewidth=2pt,linecolor=blue!40]{leftA}{rightA}
\ncline[linewidth=1pt,linecolor=blue!40,linestyle=dashed]{->}{leftA}{leftB}
\ncline[linewidth=1pt,linecolor=blue!40,linestyle=dashed]{->}{rightA}{rightB}
\end{dspPlot}
\caption{Upsampling by $3$.}\label{fig:mr:upB}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Upsampling and Interpolation}
\label{sec:mr:upfilt}
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}:
\itempar{Zero-Order Hold.}\index{zero-order hold!(discrete-time)}
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}.
\itempar{First-Order Hold.}\index{first-order hold!(discrete-time)}
In this discrete-time interpolation scheme, we obtain a piecewise linear interpolation after upsampling by $N$ by using
\begin{equation}\label{eq:mr:foi}
h_1 [n] =
\begin{cases}
\displaystyle 1 - \frac{|n|}{N} & \ |n| < N \\
0 & \mbox{ otherwise}
\end{cases}
\end{equation}
The impulse response is the familiar triangular function\footnote{
Once again, let us note that the triangle is the convolution of two rects, $h_1[n] = (1/N) \bigl(h_0[\cdot] \ast h_0[\cdot] \bigr)[n]$.}
shown in Figure~\ref{fig:mr:zoh}-(b). An example of the resulting interpolation is shown in the bottom panel of Figure~\ref{fig:mr:upinterp}.
\itempar{Sinc Interpolation.}\index{sinc interpolation!(discrete-time)}
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]$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!t]
\center
\begin{tabular}{cc}
\begin{dspPlot}[sidegap=0.5,xout=true,width=0.5\dspWidth,height=0.5\dspHeight]{-4,4}{0,1.1}
\dspSignal[linecolor=ColorDT]{x 0 ge {x 4 lt {1} {0} ifelse} {0} ifelse}
\end{dspPlot}
&
\begin{dspPlot}[sidegap=0.5,xout=true,width=0.5\dspWidth,height=0.5\dspHeight]{-4,4}{0,1.1}
\dspSignal[linecolor=ColorDT]{x abs 4 lt {1 x abs 4 div sub} { 0} ifelse}
\end{dspPlot}
\\ (a) & (b)
\end{tabular}
\caption{Discrete-time zero-order (a) and first-order (b) interpolators for $N = 4$.}\label{fig:mr:zoh}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[!b]
%
\def\dtSig{5 sub 0.3 mul RadtoDeg %
dup cos 0.3 mul exch %
0.5 mul sin 0.6 mul %
add 0.5 add }
\def\upFact{4 }
%
\center
\begin{dspPlot}[height=\dspHeightCol,xout=true,sidegap=4]{-32,32}{-0.5,1.1}
\dspSignal[linecolor=ColorDT!65]{x 4 div floor \dtSig}
\dspSignal[linecolor=ColorDT,xunit=4,xmax=8,xmin=-8]{x \dtSig}
\end{dspPlot}
\begin{dspPlot}[height=\dspHeightCol,xout=true,sidegap=4]{-32,32}{-0.5,1.1}
\dspSignal[linecolor=ColorDT!65]{%
x \upFact div floor % last upsampling point n
dup \upFact mul x exch sub \upFact div % intra-interval factor p
2 copy % n p n p
1 exch sub exch \dtSig mul % n p (1-p)f(n)
3 1 roll exch % A p n
1 add \dtSig % A p f(n+1)
mul add}
\dspSignal[linecolor=ColorDT,xunit=4,xmax=8,xmin=-8]{x \dtSig}
\end{dspPlot}
\caption{Upsampling by $4$ followed by interpolation: zero-order hold (top panel); linear interpolation (bottom panel).}\label{fig:mr:upinterp}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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$.
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Downsampling}\index{downsampling|(}
\label{sec:mr:down}
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.},
we can write
\begin{equation}
(\mathcal{S}_N \mathbf{x})[n] = x_{N\downarrow}[n] = x[nN]
\end{equation}
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b]
\center
\begin{dspBlocks}{1.3}{0.6}
$x[n]$ & \BDdwsmp{N} & $x_{N\downarrow}[n]$
\psset{arrows=->,linewidth=1pt}
\ncline{1,1}{1,2} \ncline{1,2}{1,3}
\end{dspBlocks}
\caption{Symbol for the downsampling operator.}\label{fig:mr:upSym}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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
\[
\mathcal{S}_2 \mathbf{x} = \ldots,\, x[-4],\, x[-2], \,x[0], \,x[2], \,x[4],\, \ldots
\]
On the other hand, if we delay $\mathbf{x}$ by one before subsampling we obtain
\[
\mathcal{S}_2 (\mathcal{D} \mathbf{x}) = \ldots ,\, x[-5],\, x[-3], \, x[1], \, x[3], \, x[5],\, \ldots
\]
so that clearly
\[
\mathcal{D} (\mathcal{S}_2 \mathbf{x}) \neq \mathcal{S}_2 (\mathcal{D} \mathbf{x}).
\]
The downsampling operator is {\em periodically time-varying\/} since, for all $k \in \mathbb{N}$,
\begin{equation}
\mathcal{S}_N (\mathcal{D}_{kN} \mathbf{x}) = \mathcal{D}_{kN}(\mathcal{S}_N \mathbf{x}).
\end{equation}
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
%
\def\dtSig{5 sub 0.3 mul RadtoDeg %
dup cos 0.3 mul exch %
0.5 mul sin 0.6 mul %
add 0.5 add }
%
\center
\begin{dspPlot}[height=\dspHeightCol,xout=true,sidegap=4]{-32,32}{-0.5,1.1}
\dspSignal[linecolor=ColorDT]{x \dtSig}
\end{dspPlot}
\begin{dspPlot}[height=\dspHeightCol,xout=true,sidegap=4]{-32,32}{-0.5,1.1}
\dspSignal[linecolor=ColorDT!20]{x \dtSig}
\dspSignal[linecolor=ColorDT,plotpoints=17]{x \dtSig}
%\dspSignal{x 4 mul \dtSig}
\end{dspPlot}
\begin{dspPlot}[height=\dspHeightCol,xout=true,sidegap=1]{-8,8}{-0.5,1.1}
\dspSignal[linecolor=ColorDT]{x 4 mul \dtSig}
\end{dspPlot}
\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}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Frequency-Domain Representation}
Let's consider first the \ztrans\ of a downsampled signal:
\begin{equation}
X_{N\downarrow}(z) = \sum_{n=-\infty}^{\infty} x[nN] z^{-n}.
\end{equation}
Our goal is to find a relationship between $X_{N\downarrow}(z)$ and $X(z)$; to do so let's introduce an ``auxiliary'' \ztrans\ $A(z)$ defined as
\begin{equation}
A(z) = \sum_{n=-\infty}^{\infty} x[nN] z^{-nN};
\end{equation}
clearly
\begin{equation} %\label{}
X_{N\downarrow} (z) = A(z^{1/N})
\end{equation}
so that now we need to express $A(z)$ in terms of $X(z)$, which is easier (although not entirely straightforward). The first step is to write
\begin{equation}
A(z) = \sum_{n=-\infty}^{\infty} \xi_N[n] x[n] z^{-n}
\end{equation}
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}
\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}):
\[
\boldsymbol{\xi}_N = \mbox{IDFS}_N\{\mathbf{1}\};
\]
explicitly:
\begin{equation}
\xi_N[n] = \frac{1}{N}\sum_{k=0}^{N-1} e^{j\frac{2\pi}{N}nk}.
\end{equation}
With this we can write
\begin{align}
A(z) &= \sum_{n=-\infty}^{\infty} \xi_N[n] x[n] z^{-n} \nonumber \\
&= \sum_{n=-\infty}^{\infty}\frac{1}{N} \sum_{k=0}^{N-1} e^{j\frac{2\pi}{N}nk} x[n] z^{-n} \nonumber \\
&= \frac{1}{N} \sum_{k=0}^{N-1} \sum_{n=-\infty}^{\infty} x[n] \big(e^{-j\frac{2\pi}{N}k} z \big)^{-n} \nonumber \\
&= \frac{1}{N} \sum_{k=0}^{N-1} X\big(e^{-j\frac{2\pi}{N}k} z \big)
\end{align}
so that finally:
\begin{equation}\label{eq:mr:dss}
X_{N\downarrow} (z) = \frac{1}{N} \sum_{k=0}^{N-1} X \big(e^{-j\frac{2\pi}{N}k} z^{\frac{1}{N}} \big)
\end{equation}
The Fourier transform of the downsampled signal is obtained by evaluating $X_{N\downarrow} (z)$ on the unit circle:
\begin{equation} %\label{}
X_{N\downarrow} (e^{j\omega}) = \frac{1}{N}\sum_{k=0}^{N-1} X\left(e^{j\frac{\omega - 2\pi k}{N}} \right).
\end{equation}
To understand the shape of the downsampled spectrum, let's first examine $A(z)$ on the unit circle:
\begin{equation} \label{eq:mr:nonscaled}
A(e^{j\omega}) = \frac{1}{N}\sum_{k=0}^{N-1} X\left(e^{j(\omega - \frac{2\pi}{N}k)} \right);
\end{equation}
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
\begin{align*}
X_{2\downarrow}(z) &= \frac{1}{2}\, \bigl[ X(z) + X(-z) \bigr] \\
X_{2\downarrow}(e^{j\omega}) &= \frac{1}{2} \left[X\bigl(e^{j\frac{\omega}{2}} \bigr) + X\bigl(e^{j(\frac{\omega}{2} - \pi)} \bigr) \right]
\end{align*}
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center
\begin{dspBlocks}{1.3}{0.6}
$x[n]$ & \BDfilter{LP$\{\pi/N\}$} & \BDdwsmp{N} & $y[n]$
\end{dspBlocks}
\psset{arrows=->,linewidth=1pt}
\ncline{1,1}{1,2} \ncline{1,2}{1,3}
\ncline{1,3}{1,4}
\caption{Anti-aliasing filter before downsampling.}\label{fig:mr:downfilt}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\center
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.5}}
\dspCustomTicks[axis=x]{0.5 $\omega_M$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.5}}
\dspCustomTicks[axis=x]{0.5 $\omega_M$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j(\omega - 2\pi k/N)})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.5}}
\dspFunc[linecolor=ColorDF,linestyle=dashed]{x 1 sub \dspPeriodize \dspTri{0}{0.5}}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$A(e^{j\omega})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{
x \dspPeriodize \dspTri{0}{0.5}
x 1 sub \dspPeriodize \dspTri{0}{0.5}
add 2 div}
\pnode(-.5,0){A}\pnode(0.5,0){B}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X_{2\downarrow}(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{
x 2 div \dspPeriodize \dspTri{0}{0.5}
x 2 div 1 sub \dspPeriodize \dspTri{0}{0.5}
add 2 div}
\pnode(-1,1.2){a}\pnode(1,1.2){b}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{A}{a}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{B}{b}
\end{dspPlot}
\caption{Downsampling by $2$; the highest frequency is $\pi/2$ and no aliasing occurs.}\label{fig:mr:exA}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\center
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.75}}
\dspCustomTicks[axis=x]{0.75 $\omega_M$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.75}}
\dspCustomTicks[axis=x]{0.75 $\omega_M$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j(\omega - 2\pi k/N)})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.75}}
\dspFunc[linecolor=ColorDF,linestyle=dashed]{x 1 sub \dspPeriodize \dspTri{0}{0.75}}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$A(e^{j\omega})$}]{-3,3}{0,1.2}
\dspCustomTicks[axis=x]{0.5 $\pi/2$}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x \dspPeriodize \dspTri{0}{0.75} 2 div}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x 1 sub \dspPeriodize \dspTri{0}{0.75} 2 div}
\dspFunc[linecolor=ColorDF]{
x \dspPeriodize \dspTri{0}{0.75}
x 1 sub \dspPeriodize \dspTri{0}{0.75}
add 2 div}
\pnode(-.5,0){A}\pnode(0.5,0){B}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X_{2\downarrow}(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{
x 2 div \dspPeriodize \dspTri{0}{0.75}
x 2 div 1 sub \dspPeriodize \dspTri{0}{0.75}
add 2 div}
\pnode(-1,1.2){a}\pnode(1,1.2){b}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{A}{a}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{B}{b}
\end{dspPlot}
\caption{Downsampling by $2$; the highest frequency is larger than $\pi/2$ (here, $\omega_M = 2\pi/3$) and aliasing occurs.}\label{fig:mr:exB}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\center
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.6667}}
\dspCustomTicks[axis=x]{0.6667 $\omega_M$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.6667}}
\dspCustomTicks[axis=x]{0.6667 $\omega_M$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j(\omega - 2\pi k/N)})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.6667}}
\dspFunc[linecolor=ColorDF,linestyle=dashed]{x 0.6667 sub \dspPeriodize \dspTri{0}{0.6667}}
\dspFunc[linecolor=ColorDF,linestyle=dotted]{x 1.3333 sub \dspPeriodize \dspTri{0}{0.6667}}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$A(e^{j\omega})$}]{-3,3}{0,1.2}
\dspCustomTicks[axis=x]{0.33333 $\pi/3$}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x \dspPeriodize \dspTri{0}{0.6667} 3 div}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x 0.6667 sub \dspPeriodize \dspTri{0}{0.6667} 3 div}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x 1.3333 sub \dspPeriodize \dspTri{0}{0.6667} 3 div}
\dspFunc[linecolor=ColorDF]{0.3333}
\pnode(-.3333,0){A}\pnode(0.3333,0){B}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X_{2\downarrow}(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{0.33333}
\pnode(-1,1.2){a}\pnode(1,1.2){b}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{A}{a}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{B}{b}
\end{dspPlot}
\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}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
%
\def\specHW{0.5 }
\def\specFunPar{abs 1 sub \specHW 2 mul div dup mul 1 exch sub }
\def\specFun{abs \specHW sub 1 \specHW sub div }
\def\periodSpec{%
1 sub 2 div dup floor sub 2 mul 1 sub dup dup % periodization
-\specHW lt {pop \specFunPar } {%
\specHW gt {\specFun } {pop 0} ifelse } %
ifelse }
%
\center
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=2,ylabel={$X(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \periodSpec}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \periodSpec}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j(\omega - 2\pi k/N)})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \periodSpec}
\dspFunc[linecolor=ColorDF,linestyle=dashed]{x 1 sub \periodSpec}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$A(e^{j\omega})$}]{-3,3}{0,1.2}
\dspCustomTicks[axis=x]{0.5 $\pi/2$}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x \periodSpec 2 div}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x 1 sub \periodSpec 2 div}
\dspFunc[linecolor=ColorDF]{
x \periodSpec
x 1 sub \periodSpec
add 2 div}
\pnode(-.5,0){A}\pnode(0.5,0){B}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X_{2\downarrow}(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{
x 2 div \periodSpec
x 2 div 1 sub \periodSpec
add 2 div}
\pnode(-1,1.2){a}\pnode(1,1.2){b}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{A}{a}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{B}{b}
\end{dspPlot}
\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}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
%
\def\dspTriTru#1#2#3{ #1 sub abs #3 div dup 0.5 gt {pop 0} {#3 mul #2 div 1 exch sub} ifelse }
%
\center
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTri{0}{0.75}}
\dspFunc[linecolor=ColorDFilt,linestyle=dashed]{x \dspRect{0}{1}}
\dspCustomTicks[axis=x]{0.75 $\omega_M$ 0.5 $\pi/2$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$H(e^{j\omega})X(e^{j\omega})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTriTru{0}{0.75}{1}}
\dspCustomTicks[axis=x]{0.75 $\omega_M$}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X(e^{j(\omega - 2\pi k/N)})$}]{-3,3}{0,1.2}
\dspFunc[linecolor=ColorDF]{x \dspPeriodize \dspTriTru{0}{0.75}{1}}
\dspFunc[linecolor=ColorDF,linestyle=dashed]{x 1 sub \dspPeriodize \dspTriTru{0}{0.75}{1}}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$A(e^{j\omega})$}]{-3,3}{0,1.2}
\dspCustomTicks[axis=x]{0.5 $\pi/2$}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x \dspPeriodize \dspTriTru{0}{0.75}{1} 2 div}
\dspFunc[linecolor=ColorDF!50,linestyle=dashed]{x 1 sub \dspPeriodize \dspTriTru{0}{0.75}{1} 2 div}
\dspFunc[linecolor=ColorDF]{
x \dspPeriodize \dspTriTru{0}{0.75}{1}
x 1 sub \dspPeriodize \dspTriTru{0}{0.75}{1}
add 2 div}
\pnode(-.5,0){A}\pnode(0.5,0){B}
\end{dspPlot}
\begin{dspPlot}[xtype=freq,height=\dspHeightCol,xticks=1,ylabel={$X_{2\downarrow}(e^{j\omega})$}]{-1,1}{0,1.2}
\dspFunc[linecolor=ColorDF]{
x 2 div \dspPeriodize \dspTriTru{0}{0.75}{1}
x 2 div 1 sub \dspPeriodize \dspTriTru{0}{0.75}{1}
add 2 div}
\pnode(-1,1.2){a}\pnode(1,1.2){b}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{A}{a}
\ncline[linewidth=1pt,linecolor=gray,linestyle=dashed]{->}{B}{b}
\end{dspPlot}
\caption{Downsampling by $2$ with preliminary anti-aliasing filtering.}\label{fig:mr:exE}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Fractional Sampling Rate Changes}
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
a simple example:
\begin{align*}
(\mathcal{S}_2 ( \mathcal{U}_2 \mathbf{x}))[n] &= x[n] \\
(\mathcal{U}_2 ( \mathcal{S}_2 \mathbf{x}))[n] &=
\begin{cases}
x[n] & \mbox{ for $n$ even} \\
0 & \mbox{ for $n$ odd}
\end{cases}.
\end{align*}
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}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\center
\begin{dspBlocks}{1}{0.6}
$x[n]$ & \BDupsmp{N} & \BDfilter{LP$\{\min(\pi/N,\pi/M)\}$} & \BDdwsmp{M} & $y[n]$
\psset{arrows=->,linewidth=1pt}
\ncline{1,1}{1,2} \ncline{1,2}{1,3}
\ncline{1,3}{1,4}\ncline{1,4}{1,5}
\end{dspBlocks}
\caption{Block diagram for a fractional sampling rate change.}\label{fig:mr:frac}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b]
\center
\begin{dspBlocks}{1}{0.6}
$x[n]$ & \BDupsmp{4} & \BDfilter{LP$\{\pi/4\}$} & \BDdwsmp{3} & $y[n]$ \\
\\
$x[n]$ & \BDupsmp{3} & \BDfilter{LP$\{\pi/4\}$} & \BDdwsmp{4} & $y[n]$ \\
\psset{arrows=->,linewidth=1pt}
\ncline{1,1}{1,2} \ncline{1,2}{1,3}
\ncline{1,3}{1,4}\ncline{1,4}{1,5}
\ncline{3,1}{3,2} \ncline{3,2}{3,3}
\ncline{3,3}{3,4}\ncline{3,4}{3,5}
\end{dspBlocks}
\caption{Block diagram for a fractional sampling rate change.}\label{fig:mr:fracex}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Practical resamplers}
\index{resampling}
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}
\begin{equation} \label{eq:mr:LagInterp}
l_x(n; t) = \sum_{k=-N}^{N} x[n - k] \,L_k^{(N)}(t)
\end{equation}
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}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% finite-length signal, 3 taps
\def\ta{.5 } \def\tb{2 } \def\tc{-1 }
%
%% the tap plotting string
\def\taps{-1 \ta 0 \tb 1 \tc }
\def\plotTaps{\dspTaps[linecolor=ColorDT]{\taps}}
%
%% Lagrange polynomials
\def\lpa{%
dup -1 div exch %
-1 add -2 div %
mul }
\def\lpb{%
dup 1 add -1 div exch %
-1 add %
mul }
\def\lpc{%
dup exch %
1 add 2 div %
mul }
%
\def\lagInterp{%
dup \lpa \ta mul exch %
dup \lpb \tb mul exch %
\lpc \tc mul %
add add}
%
\def\LagPoly#1#2#3#4{%
\dspFunc[linecolor=#2]{x \csname lp#3\endcsname}%
\dspText(! #1 2 sub 2 div 1.3){\color{#2} $L_{#4}^{1}(t)$}}
%
\def\sample#1#2#3{%
\psline[linewidth=\dspStemWidth,#1](#2,0)(! #2 dup #3)
\psdot[dotstyle=*,dotsize=\dspDotSize,#1](! #2 dup #3)}
%
\begin{figure}[t!]
\center
\begin{dspPlot}[sidegap=0]{-1.5,1.5}{-0.8,1.6}
\begin{dspClip}%
\LagPoly{1}{green}{a}{-1}%
\LagPoly{2}{blue}{b}{0}%
\LagPoly{3}{orange}{c}{1}%
\end{dspClip}
\end{dspPlot}
\caption{Quadratic (second-order) Lagrange interpolation polynomials $L_k^{(1)}(t)$.}\label{fig:mr:lagPoly}
\end{figure}
\begin{figure}[t!]
\begin{dspPlot}[sidegap=0,xticks=custom]{-1.5,1.5}{-2,3}
\begin{dspClip}%
\psset{linecolor=lightgray}
\dspFunc[linewidth=0.4pt]{x \lpa \ta mul}%
\dspFunc[linewidth=0.4pt]{x \lpb \tb mul}%
\dspFunc[linewidth=0.4pt]{x \lpc \tc mul}%
\dspFunc[linewidth=2pt,linecolor=ColorCT,xmin=-2,xmax=2]{x \lagInterp}%
\end{dspClip}
\plotTaps%
\sample{linecolor=ColorDT2}{0.3}{\lagInterp}%
\dspCustomTicks[axis=x]{-1 $n-1$ 0 $n$ 1 $n+1$ 0.3 $n+\tau$}%
\end{dspPlot}
\caption{Local quadratic Lagrange interpolation $l_x(n; t)$ around $n$ and approximation $x_\tau[n] = l_x(n; \tau)$ for $\tau=0.3$}\label{fig:mr:approx}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Once the value for $\tau$ is known, it can be plugged into~(\ref{eq:mr:LagInterp}) so that~(\ref{eq:mr:subApprox}) becomes
\begin{equation}\label{eq:mr:FIRapprox}
x_\tau[n] = \sum_{k=-N}^{N} x[n - k] \,L_k^{(N)}(\tau) = \left(\mathbf{d}_\tau \ast \mathbf{x}\right)[n]
\end{equation}
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$.}
\begin{align}
m(n) &= \mbox{nint}\left( n\beta \right) \, \in \mathbb{N} \\
|\tau(n)| & = \left| n\beta - m(n) \right| \leq \frac{1}{2} \label{eq:mr:taun}
\end{align}
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
\[
y[n] \approx \sum_{k=-N}^{N} x[m(n) - k] \,\hat{d}_{\tau(n)}[k]
\]
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]$.
\item $n=1$: $m(1) = \mbox{nint}(0.78) = 1$, $\tau(1) = \beta - m(1) = -0.22$ so $y[1] = x_{-0.22}[1]$
\item $n=2$: $m(2) = \mbox{nint}(2\beta) = \mbox{nint}(1.56) = 2$, $\tau(2) = -0.44$ so $y[2] = x_{-0.44}[2]$
\item $n=3$: $m(3) = \mbox{nint}(3\beta) = \mbox{nint}(2.34) = 2$ and therefore we must \textit{reuse} the previous interpolation anchor $x[2]$; $\tau(3) = 2.34 - 2 = 0.34$ so $y[3] = x_{0.34}[2]$
\item $n=4$: $m(4) = \mbox{nint}(3.12) = 3$, $\tau(4) = 0.12$ so $y[2] = x_{0.12}[3]$
\item $\ldots$
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[t]
\def\polyFun{%
dup 2.315 mul exch %
dup dup mul -3.74583333 mul exch %
dup dup dup mul mul 2.05 mul exch %
dup dup dup dup mul mul mul -0.45416667 mul exch %
dup dup dup dup mul mul mul mul 0.035 mul %
0.5 add add add add add %
3 mul }
%
\def\subs#1{%
\psline[linecolor=ColorDT2,linewidth=\dspStemWidth](#1,0)(! #1 dup \polyFun)
\psdot[linecolor=ColorDT2,dotstyle=*,dotsize=\dspDotSize](! #1 dup \polyFun)}
%
\center
\begin{dspPlot}[sidegap=0.2,yticks=none,xout=true]{0,4}{0,3}
\dspFunc[linecolor=ColorCT,linestyle=dashed,xmin=-.12,xmax=0]{x \polyFun}
\dspFunc[linecolor=ColorCT,linestyle=dashed,xmin=4,xmax=4.1]{x \polyFun}
\dspFunc[linecolor=ColorCT,xmin=0,xmax=4]{x \polyFun}
\multips(0,0)(0.78,0){6}{\psline[linewidth=0.5pt,linestyle=dashed,linecolor=lightgray](0,3)}
\subs{0}
\subs{0.78}
\subs{1.56}
\subs{2.34}
\subs{3.12}
\subs{3.9}
\dspSignal[linecolor=ColorDT]{x \polyFun}
\rput*[b](.89,.38){$\tau(1)$}\psline[arrows=->](1,.3)(.78,.3)
\rput*[b](1.78,.38){$\tau(2)$}\psline[arrows=->](2,.3)(1.56,.3)
\rput*[b](2.18,.38){$\tau(3)$}\psline[arrows=->](2,.3)(2.34,.3)
\rput*[b](3.18,.38){$\tau(4)$}\psline[arrows=->](3,.3)(3.12,.3)
\end{dspPlot}
\caption{On-the-fly fractional resampling with a 28\% frequency increase ($\beta = 0.78$).}\label{fig:mr:resampleup}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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]$.
\item $n=1$: $m(1) = \mbox{nint}(1.22) = 1$, $\tau(1) = \beta - m(1) = 0.22$ so $y[1] = x_{0.22}[1]$
\item $n=2$: $m(2) = \mbox{nint}(2\beta) = \mbox{nint}(2.44) = 2$, $\tau(2) = 0.44$ so $y[2] = x_{0.44}[2]$
\item $n=3$: $m(3) = \mbox{nint}(3\beta) = \mbox{nint}(3.66) = 4$ and therefore we must \textit{skip} $x[3]$; $\tau(3) = 3.66 - 4 = -0.34$ so $y[3] = x_{-0.34}[4]$
\item $n=4$: $m(4) = \mbox{nint}(3.12) = 3$, $\tau(4) = 0.12$ so $y[2] = x_{0.12}[3]$
\item $\ldots$
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[b]
\def\polyFun{%
dup 2.315 mul exch %
dup dup mul -3.74583333 mul exch %
dup dup dup mul mul 2.05 mul exch %
dup dup dup dup mul mul mul -0.45416667 mul exch %
dup dup dup dup mul mul mul mul 0.035 mul %
0.5 add add add add add %
3 mul }
%
\def\subs#1{%
\psline[linecolor=ColorDT2,linewidth=\dspStemWidth](#1,0)(! #1 dup \polyFun)
\psdot[linecolor=ColorDT2,dotstyle=*,dotsize=\dspDotSize](! #1 dup \polyFun)}
%
\center
\begin{dspPlot}[sidegap=0.2,yticks=none,xout=true]{0,4}{0,3}
\dspFunc[linecolor=ColorCT,linestyle=dashed,xmin=-.12,xmax=0]{x \polyFun}
\dspFunc[linecolor=ColorCT,linestyle=dashed,xmin=4,xmax=4.1]{x \polyFun}
\dspFunc[linecolor=ColorCT,xmin=0,xmax=4]{x \polyFun}
\multips(0,0)(1.22,0){4}{\psline[linewidth=0.5pt,linestyle=dashed,linecolor=lightgray](0,3)}
\subs{0}
\subs{1.22}
\subs{2.44}
\subs{3.66}
\dspSignal[linecolor=ColorDT]{x \polyFun}
\rput*[b](1.11,.38){$\tau(1)$} \psline[arrows=->](1,.3)(1.22,.3)
\rput*[b](2.22,.38){$\tau(2)$} \psline[arrows=->](2,.3)(2.44,.3)
\rput*[b](3.83,.38){$\tau(3)$} \psline[arrows=<-](3.66,.3)(4,.3)
\end{dspPlot}
\caption{On-the-fly fractional resampling with a 18\% frequency reduction ($\beta = 1.22$).}\label{fig:mr:resampledown}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
c4science · Help