Interpolation is the procedure with which we convert a discrete-time sequence $x[n]$ to a continuous-time function $x(t)$. 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 the 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-length signal 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 $x(t)$ 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[(b)] the interpolation does not ``go through'' the data points; indeed, it seems intuitive that the interpolation should respect the original dataset;
\item[(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[(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 ``minimalistic'' 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 $x[n]$, once we have chosen an interpolation interval $T_s$, the first requirement for the interpolating function $x(t)$ is that its values at multiples of $T_s$ be equal to the corresponding points of the discrete-time sequence, i.e.
\[
x(t)\Bigr|_{t = nT_s} = x[n];
\]
with our choice of $T_s =1$ the above condition becomes simply $x(n)= x[n]$. 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 $x(t)$ to be in $C^{\infty}$, where $C^{M}$ is the class of functions for which all derivatives up to order $M$ exist and are continuous.
One 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 interpolator's $2N+1$ coefficients 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 write $p(n)$ as a linear combination of basis vector for the space. This is not the first time we encounter the need for polynomial basis: in Chapter~\ref{ch:vs} we briefly encountered the Legendre polynomial basis in the context of approximation, while the interpolator in~(\ref{eq:is:straightPI}) is clearly expressed as a linear combination of the $2N+1${\em monomial}\/ basis vectors $\{1, t, t^2, \ldots t^{2N}\}$. For the task at hand, however, 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(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~(\ref{eq:is:lagInterp}) is indeed equal to $x[n]$ for $2N+1$ integer values of $t$, 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 two serious drawbacks. Firstly, although the Lagrangian basis provides a way to quickly compute the necessary coefficients, the set of basis vectors is completely different for every different length of the input data. From an engineering point of view it is obvious that we would like a more universal interpolation machine rather than having to redesign our interpolator as a function of problem size. The second difficulty, that is common to all types of polynomial fitting, is that the method becomes numerically ill-conditioned as the polynomial degree grows large.
A way out of these problems is to relax at least partially the smoothness constraint; this leads to the simpler interpolation schemes that we will look at in the next Section.
If we relax the smoothness constraint, we can design much simpler interpolation schemes. 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 significative:
\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.
\end{itemize}
In summary, the ZOH therefore creates a continuous-time signal by stitching together delayed and scaled versions of the \textit{same} interpolation function, also known as the interpolation kernel, independently of the amount of data to interpolate. The resulting ``interpolation machine'' is both conceptually and practically a much simpler device. Of course the price to pay for this simplicity is a discontinuous signal, but we can easily extend the local interpolation paradigm to improve the result.
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.
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 later.
The tradeoff seems however rather 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.
\end{itemize}
Interestingly, however, as the size of the data set increases to infinity, the two methods become one and the same! This is perhaps one of the most remarkable mathematical results in discrete-time signal processing, and one that will lead naturally to the sampling theorem.
\subsection{Sinc Interpolation}
\label{sec:is:sincinterp}
Consider again the expression for the Lagrange polynomial (\ref{eq:is:lagPoly}) and let's try to determine what happens when the size of the data set to interpolate (and therefore the degree of the Lagrange basis vectors) goes to infinity:
\[
\lim_{N\rightarrow\infty} L^{(N)}_n(t)=\prod_{ {\scriptstyle k =-\infty\atop\scriptstyle k\neq n }}^{\infty} \frac{t - k}{n - k} =\ ?
\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}. In other words, the maximally smooth interpolation for an infinite-length sequence can be obtained using a kernel-based method where the kernel is the sinc function:
Obviously, since the sinc is a two-sided, infinite support function, sinc interpolation falls in the same unattainable category as ideal filters: something that we can try to approximate but that we will never be able to exactly implement in practice. Nevertheless, 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.
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 (in the figure, $T_s =1$). Note how the interpolation property of the sinc that we introduced in Section~\ref{sec:is:SincProp} guarantees that $x(nT_s)= x[n]$.
\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}
\itempar{Spectral Properties of the Sinc Interpolation.}
The interpolation kernel $\sinc(t/T_s)$ is $F_s$-bandlimited, with $F_s =1/T_s$. Therefore the sinc interpolation of a discrete-time sequence, being a linear combination of bandlimited functions, is itself bandlimited. If the DTFT $X(e^{j\omega})$ of the discrete-time sequence is well defined, the spectrum of the interpolated signal can be obtained as follows:
(1/F_s) \, X(e^{j2\pi f/F_s}) &\mbox{ for $ |f| \leq F_s/2$}\\
0 &\mbox{ otherwise}
\end{cases}\label{eq:is:interpSpec}
\end{align}
In other words, the continuous-time spectrum is just a scaled and stretched version of the DTFT of the discrete-time sequence 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).\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.} Figure~\ref{fig:is:interpSpeed} illustrates the concept graphically.
\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}