Interpolation is the procedure by which we convert a discrete-time sequence $\mathbf{x}$ to a continuous-time function $\mathbf{x}_c$. At the core of the interpolation procedure, as we have mentioned, is the association of a {\em physical}\/ time duration $T_s$ to the intervals between samples in the discrete-time sequence; next, we need to ``fill the gaps'' between sample values. To develop the intuition, let's consider a simple example where we want to interpolate a finite set of values as in Figure~\ref{fig:is:diffint}-(a); Figures~\ref{fig:is:diffint}-(b), (c) and (d) show three possible interpolations of the given dataset to a continuous-time signal $\mathbf{x}_c$ using $T_s = 1$.
Although we haven't really formalized the requirements for the interpolation process, it is rather intuitive that all three proposed interpolations are not entirely satisfactory:
\begin{itemize}
\item in (b) the interpolation does not go through the data points; indeed, it seems reasonable to require that the interpolation include the original data set exactly;
\item in (c), while the function goes through the data points, it does not look smooth; in general, we are right to be suspicious of non-smooth curves, since no natural phenomenon can create discontinuous jumps in amplitude (achievable only via infinite speed) or curvature (achievable only via infinite acceleration);
\item in (d), while the function goes through the data points and is smooth, it seems to ``wander around'' too much, as if interpolating a more complicated dataset; indeed, we want our interpolation strategy to be ``economical'' in using the information contained in the discrete-time signal.
\end{itemize}
In the following sections we will try to make our requirements more precise and see how we can arrive at a universal interpolation scheme that is both intuitively and mathematically sound.
%To formalize the problem, and without loss of generality, we will assume in the following that $T_s = 1$ and that the finite-length discrete-time signal is defined for $n = -N, \ldots, 0, \ldots, N$ for an odd length of $M = 2N+1$ points.
Given a discrete-time signal $\mathbf{x}$, once we have chosen an interpolation interval $T_s$, the first requirement for the interpolating function $\mathbf{x}_c$ is that its values at multiples of $T_s$ be equal to the corresponding points of the discrete-time sequence, i.e.
\[
x_c(nT_s) = x[n].
\]
%with our choice of $T_s = 1$ the above condition becomes simply $x(n) = x[n]$. T
The second requirement is that we want the interpolating function to be smooth and, mathematically, the smoothness of a function increases with the number of its continuous derivatives. For maximal smoothness, therefore, we require $\mathbf{x}_c \in C^{\infty}$, where $C^{M}$ is the class of functions for which all derivatives up to order $M$ exist and are continuous.
For now, let's consider the case where $\mathbf{x}$ is a finite-length signal; without loss of generality, assume the signal has odd lenght $2N+1$ and that the signal is centered in zero, as in Figure\ref{fig:is:diffint}-(a) before. To lighten the notation, let's also pick $T_s = 1$ for now.
The simplest choice for a maximally differentiable curve through a set of $2N+1$ data points is the unique polynomial interpolator of degree $2N$
which, like all polynomials, belongs to $C^{\infty}$. Computation of the $2N+1$ coefficients $a_k$ is a classic algebraic problem, first solved in the 17th century by Newton. Numerically, one way to arrive at the solution is to work out the system of $2N+1$ equations
which can be carried out algorithmically using Pascal's triangle, for instance. A more interesting approach is to consider the space of finite-degree polynomials over the interval $[-N, N]$ and choose a ``smart'' basis to express $\mathbf{x}_c$. This is not the first time we encounter the idea of using a polynomial basis: in Chapter~\ref{ch:vs}, for instance, we introduced the Legendre basis to solve a function approximation problem. The interpolator in~(\ref{eq:is:straightPI}) clearly expresses $\mathbf{x}_c$ as a linear combination of the $2N+1$ {\em monomial}\/ basis vectors $\{1, t, t^2, \ldots t^{2N}\}$ but, for the task at hand, a more appropriate basis is the set of {\em Lagrange polynomials}.
The Lagrange polynomial basis for the interval $I=[-N, N]$ is the family of $2N+1$ polynomials
\begin{equation} \label{eq:is:lagPoly}
L^{(N)}_n(t) = \mathop{\prod_{k = -N}}_{k\neq n}^{N}\frac{t - k}{n - k}, \qquad n = -N, \ldots, N
\end{equation}
each of which has degree $2N$. As an example, the family of five polynomials for $N=2$ is shown in Figure~\ref{fig:is:LagBasis}. A key property of the Lagrangian basis vector is that, for integer values of their argument, we have
\begin{equation} \label{eq:is:lagInterpProp}
L^{(N)}_n(m) = \left\{
\begin{array}{ll}
1 & \mbox{if $n=m$}\\
0 & \mbox{if $n \neq m$}
\end{array}\right. \qquad\qquad -N \leq n,m \leq N.
\end{equation}
Using the above result, it is easy to verify that the polynomial interpolator for a discrete-time signal of length $2N+1$ is simply
\begin{equation} \label{eq:is:lagInterp}
x_c(t) = \sum_{n = -N}^{N} x[n]L^{(N)}_n(t)
\end{equation}
that is, a linear combination of the Lagrangian basis vectors where the scalar coefficients are the discrete-time samples themselves. Indeed, since a polynomial of degree $M$ is uniquely determined by $M$ of its values and since $x_c(n)$, because of~(\ref{eq:is:lagInterpProp}), is equal to $x[n]$ for $n = -N, -N+1, \ldots, N$, the interpolator is indeed the unique solution to the problem.
\caption{Lagrange interpolation polynomials $L_n^{(2)}(t)$ for $n=-2,\ldots,2$. Note that $L_n^{(N)}(t)$ is zero for $t$ integer except for $t = n$, where it is $1$.}\label{fig:is:LagBasis}
As an example, Figure~\ref{fig:is:weightedLag} shows how the five polynomials $L_n^{(2)}(t)$ beautifully come together to interpolate a $5$-point discrete-time signal into the smooth curve in Figure~\ref{fig:is:finalLag}. A fundamental characteristic of polynomial interpolation is that it is {\it global}: for all values of $t$, \textit{all} the original discrete-time data points contribute to the instantaneous value of $x(t)$.
Although it elegantly solves the problem of finding the smoothest curve through a finite data set, polynomial interpolation suffers from a series of drawbacks. First of all, although the Lagrangian basis provides a way to quickly compute the interpolator's coefficients, if the number of data points changes, then the set of basis vectors needs to be redetermined from scratch. From an engineering point of view it is obvious that we would like a more universal ``interpolation machine'' that does not depend on the size of the input. Additionally, and this is a problem common to all types of polynomial fitting, the method becomes numerically ill-conditioned as the polynomial degree grows large. Finally, the method produces a continuous-time interpolator that does not admit a simple frequency-domain representation, since the resulting function diverges outside the $[-N, N]$ interval.
We will now introduce an interpolation method that solves all of these problems, but the price to pay will be the partial relaxation of the smoothness constraint.
The idea behind kernel-based interpolation is to build a continuous-time signal by summing together scaled copies of the same compact-support prototype function (the ``kernel''). Consider for instance a device that builds a continuous-time signal simply by keeping its output constant and equal to the last discrete-time input data point between interpolation intervals:
\begin{equation*}
x_0(t) = x[\, \lfloor t + 0.5 \rfloor \,], \qquad -N \leq t \leq N;
\end{equation*}
the result of this interpolation method, known as Zero-Order Hold (ZOH), is shown in Figure~\ref{fig:is:zoh} for the same 5-point signal as in Figure~\ref{fig:is:diffint}-(a). The above expression can be rewritten as
\begin{equation} \label{eq:is:zoh}
x_0(t) = \sum_{n = -N}^{N}x[n]\rect(t - n),
\end{equation}
which shows a remarkable similarity to~(\ref{eq:is:lagInterp}). The differences, however, are extremely significant:
\begin{itemize}
\item the continuous-time term in the sum (i.e the rect function) is no longer dependent on the length of the original data set
\item the dependence of the continuous-time term on interpolation interval (that is, on $n$) is only via a simple time shift
\item the value of the output at any instant $t$ is dependent only on one input data point: the interpolation is \textit{local} rather than global and can be performed in real time as the discrete-time data flows in.
\end{itemize}
The ZOH creates a continuous-time signal by stitching together delayed and scaled versions of a rectangular kernel, independently of the amount of data to interpolate: the resulting ``interpolation machine'' is both conceptually and practically a very simple device. The price to pay for this simplicity, however, is a ``jumpy'' output signal (that is, the interpolator is discontinuous).
Generally speaking, local kernel-based interpolators produce a signal via what we could call a ``mixed-domain convolution''
\begin{equation} \label{eq:is:interplocal}
x(t) = \sum_{n = -N}^{N}x[n]i(t - n);
\end{equation}
where the kernel $i(t)$ is a compact-support function fulfilling the properties
\begin{equation*}
\left\{\begin{array}{ll}
i(0) & = 1 \\
i(t) & = 0 \quad \mbox{for $t$ a nonzero integer}
\end{array}\right.
\end{equation*}
(compare this to~(\ref{eq:is:lagInterpProp})). By crafting a good kernel we can achieve an interpolated signal with better continuity properties and we will now explore a few options; note that, for an arbitrary interpolation interval $T_s$, the interpolation is obtained simply by using a scaled kernel:
The kernel used in the ZOH is, as we have seen, the simple rect function shown in Figure~\ref{fig:is:interpKernels}-(a):
\begin{equation}
i_0(t) = \begin{cases}
1 & \mbox{ if } |t| < 1/2 \\
0 & \mbox{ otherwise}
\end{cases}
\end{equation}
The kernel is discontinuous (that is, it belongs to $C^0$) and so is the interpolated signal. Only one value of the input is used to produce the current value of the output.
The kernel used in first-order interpolation is the triangle function shown in Figure~\ref{fig:is:interpKernels}-(b):
\begin{equation}
i_1(t) = \begin{cases}
1- |t| & \mbox{ if } |t| < 1 \\
0 & \mbox{ otherwise}
\end{cases}
\end{equation}
The kernel belongs to $C^1$ and the same holds for the output signal. An example of first-order interpolation applied to the usual five-point data set is shown in Figure~\ref{fig:is:firstOrdInterp}. The interpolation strategy is akin to ``connecting the dots'' between discrete-time samples so that, at any time $t$, \textit{two} discrete-time samples contribute to the output value. If implemented in real time, first-order interpolation would therefore introduce a delay of one interpolation period.
Although the first-order interpolator is only marginally more complex than the ZOH, the continuity of the output implies a much better quality of the result as can be appreciated visually for familiar signals as in Figure~\ref{fig:is:interpCompare}.
In third-order interpolation\footnote{Even-order kernels are not used in practice since they are not symmetric around a central interpolation instant.} a commonly used kernel is the cubic piecewise function shown in Figure~\ref{fig:is:interpKernels}-(c):
The kernel belongs to $C^2$ (first and second derivatives are continuous) and the same holds for the output signal. The kernel's support is four and, as a consequence, the values of the continuous-time signal depend on four neighboring discrete-time samples. An example of third-order interpolation applied to the usual five-point data set is shown in Figure~\ref{fig:is:thirdOrdInterp}.
Local interpolation schemes can be extended to higher order kernels and, in general, a kernel of order $2N+1$ is composed of polynomials of degree $2N+1$ over a support of length $2N+2$. The resulting interpolation belongs to $C^{2N}$ and the lack of continuity in the derivative of order $2N+1$ ultimately will cause undesired high frequency content in the output, as we will see momentarily.
\subsection{Ideal Sinc Interpolation}
\label{sec:is:sincinterp}
The tradeoff, so far, seems clear:
\begin{itemize}
\item global interpolation schemes such as polynomial interpolation provide maximum smoothness at the price of a complex procedure that depends on the length of the data set;
\item local kernel-based methods are simple to implement but ultimately lack in smoothness (which leads to unwanted spectral artifacts).
\end{itemize}
However, a small miracle is in store: in the limit, as the size of the data set grows to infinity, polynomial interpolation and kernel-based interpolation become one and the same, leading to what is perhaps one of the most remarkable mathematical results in discrete-time signal processing.
Consider again the expression for the Lagrange polynomial (\ref{eq:is:lagPoly}) for $T_s = 1$ and let's try to determine what happens when the size of the data set, and therefore the degree of the Lagrange polynomials, goes to infinity:
\lim_{N\rightarrow\infty} L^{(N)}_n(t) = \sinc\left(t - n \right)
\end{equation}
Remarkably, as $N$ goes to infinity, all Lagrange polynomials converge to simple time shifts of the same function and the function in question is the sinc; a graphical illustration of the convergence process for $L^{(N)}_0(t)$ is shown in Figure~\ref{fig:is:Lag2Sinc}. This means that the maximally smooth interpolation for an infinite-length sequence can be obtained using a kernel-based method where the kernel is the sinc function:
Figure~\ref{fig:is:sincInterp} shows how the scaled and time-shifted copies of the sinc kernel come together to interpolated a discrete-time sequence; note how the interpolation property of the sinc that we introduced in Section~\ref{sec:is:SincProp} guarantees that $x_c(nT_s) = x[n]$.
Obviously, since the sinc is a two-sided, infinite support function, sinc interpolation falls in the same abstract category as ideal filters: a paradigm that we can try to approximate but that cannot be exactly implement in practice. Still, the mathematical construction represents a fundamental cornerstone of the translation machinery between discrete and continuous time and will represent the starting point for the development of sampling theory.
\caption{Sinc interpolation: discrete-time signal (top); first four scaled copies of the sinc kernel for $n = 0,1,2,3$ (center); final interpolation (bottom).}\label{fig:is:sincInterp}
The resulting spectrum is therefore the product of three factors:
\begin{enumerate}
\item a scaling factor $T_s = 1/F_s$;
\item the term $X(e^{j2\pi f/F_s})$, which is the DTFT of the original sequence, rescaled so that $\omega = \pi$ in the discrete-time spectrum is mapped to $f = F_s/2$ in the continuous-time spectrum; please note that, since the DTFT is $2\pi$-periodic, then $X(e^{j2\pi f/F_s})$ is an $F_s$-periodic function;
\item $I(f/F_s)$, the Fourier transform of the rescaled kernel, which acts as a continuous-time filter.
\end{enumerate}
\itempar{Ideal sinc interpolation.} The ideal interpolation kernel $\sinc(t/T_s)$ is $F_s$-bandlimited, with $F_s = 1/T_s$ and therefore the sinc interpolation of a discrete-time sequence is itself bandlimited:
In other words, when using sinc interpolation, the continuous-time spectrum is just a scaled and stretched version of the DTFT between $-\pi$ and $\pi$. The duration of the interpolation interval $T_s$ is inversely proportional to the resulting bandwidth of the interpolated signal; intuitively, a slow interpolation ($T_s$ large) results in a spectrum concentrated around the low frequencies; conversely, a fast interpolation ($T_s$ small) results in a spread-out spectrum (more high frequencies are present); figure~\ref{fig:is:interpSpeed} illustrates the concept graphically.\footnote{
%
To find a simple everyday (yesterday's?) analogy, think of a $45$ rpm vinyl\index{vinyl} record played at either $33$~rpm (slow interpolation) or at $78$ rpm (fast interpolation) and remember the acoustic effect on the sounds.}
\caption{Spectral characteristics of a sinc-interpolated signal: DTFT of the original discrete-time signal $x[n]$ (top); Fourier transform of its sinc interpolation with interpolation interval $T_s = T_0$ (second panel); spectrum of a sinc interpolation using a slower interpolation rate ($T_s = 2T_0$): spectrum concentrates in the low frequencies around the origin (third panel); spectrum of a sinc interpolation using a faster interpolation rate ($T_s = T_0/2$): spectrum spreads out to higher frequencies (third panel); the spectral shape and total energy remain the same.}\label{fig:is:interpSpeed}
\itempar{Practical interpolation.} In practical interpolation schemes, where the kernel is a compact-suport function, the spectrum of the interpolation will be dependent on the Fourier transform of the kernel and will exhibit unwanted artifact with respect to the ideal interpolation scheme. Figure~\ref{fig:is:zohfreq}, for instance, shows the factors in~(\ref{eq:is:interpSpec}) and their product when the kernel is the zero-order hold. We can remark the following:
\begin{itemize}
\item the interpolator acts as a non-ideal lowpass filter that largely preserves the baseband copy of the periodic spectrum but only attenuates the remaining copies;
\item the interpolator distorts the baseband component since its response is not flat;
\item the spectral width of the baseband component is again inversely proportional to the interpolation period.
Higher-order interpolators ameliorate the situation; a first-order interpolator, for instance, whose kernel can be expressed as
\[
\mathbf{i}_1 = \mathbf{i}_0 \ast \mathbf{i}_0
\]
has a sharper characteristic since
\begin{equation}
I_1(f) = \sinc^2(2f).
\end{equation}
As a result, as shown in Figure~\ref{fig:is:zohfreq}, it appears that the smoothness property that made sense in the time domain also leads to a continuous-time spectrum that more closely mirrors its discrete-time counterpart by rejecting most of the out-of-band energy. It is clear that we would like the kernel to approximate as well as possible an ideal lowpass filter characteristic with cutoff $F_s/2$, which leads once again to the ideal sinc interpolator!