The Whittaker-Shannon interpolation formula (inferring missing values from a sampled signal using sinc or Lanczos function)

In an earlier blog post we had talked about Lagrange polynomial interpolation which is when you need to find a polynomial that passes through a set of points. This time we will be seeing a similar problem that is encountered in signal processing.

In discrete signal processing, a continuous signal is represented by a discrete signal consisting of values taken at regular intervals in time, a process called sampling. Below we have an example of a continuous signal that is then sampled at intervals of 0.9 (arbitrarily chosen here).





Now it may be the case that you would like to reconstruct the original signal from the sampled signal. You may think that there is an infinite number of signals that can fit the same samples but because the discrete values are sampled at regular intervals, there is actually only one signal that can fit it, provided that the sampling rate is greater than twice the highest frequency of the original signal. If this condition is met then we can construct the original signal using the Whittaker-Shannon interpolation formula.

The process is similar to that of the Lagrange polynomial interpolation. In polynomial interpolation, the trick is to find a set of component polynomials each of which passes through one of the given set of points. In order to be able to combine them all together into a single polynomial that passes through all of the given set of points, the component polynomials also need to equal zero wherever there is a different point. This makes sure that adding the polynomials together will not cause an interference between them at the given set of points.

In Whittaker-Shannon interpolation, the components are the sinc function instead of polynomials. Sinc functions are periodic functions that are as follows:

$$
\text{sinc}(x) = \begin{cases}
1 & \text{for } x = 0 \\
\frac{sin(x)}{x} & \text{otherwise} \\
\end{cases}
$$



See how sinc is equal to 1 at x = 0 and equal to zero at regular intervals around it? This is how we'll take advantage of it as a component. The first thing we need to do is to synchronise the intervals at which it equals zero to the intervals of the sampled signal. If the sample interval, or period, is 0.9, then the synchronised sinc function would be:

$$
y = sinc\left(\frac{\pi x}{0.9}\right)
$$



Next we need to shift the center peak of the synchronised sinc over one of the samples. Let's pick sample number 3 where sample number 0 is the first sample:

$$
y = sinc\left(\frac{\pi (x - 3 \times 0.9)}{0.9}\right)
$$



Finally we need to adjust the height of the peak to be equal to that of that sample value. Since the peak has a value of one and all the other important points on the sinc function are equal to zero, then if the all we have to do is multiply the sinc function by the sample value:

$$
y = -0.3303 \times sinc\left(\frac{\pi (x - 3 \times 0.9)}{0.9}\right)
$$



Great! Now we just need to do the same for all the other samples and add up all the sinc functions into a single function:

$$
y =
0 \times sinc\left(\frac{\pi (x - 0 \times 0.9)}{0.9}\right) +
0.7628 \times sinc\left(\frac{\pi (x - 1 \times 0.9)}{0.9}\right) +
-0.4309 \times sinc\left(\frac{\pi (x - 2 \times 0.9)}{0.9}\right) +
-0.3303 \times sinc\left(\frac{\pi (x - 3 \times 0.9)}{0.9}\right) +
\dots
$$



Note that the point made at the top about how there is only one signal that passes through the given samples only applied when there is an infinite supply of such samples.



Now, the sinc function is only useful when you want to make use of all points together in order to infer any missing point. In practice, this is very computationally expensive and it would be better to instead use the points thatre close to the missing point rather than all the points. To do this, instead of the sinc function we use a truncated sinc function called a Lanczos filter:

$$
\text{lanczos}(x, a) = \begin{cases}
sinc(\pi x) \times sinc(\frac{\pi x}{a}) & -a < x < a \\
0 & \text{otherwise} \\
\end{cases}
$$

where $a$ is some positive whole number such as 2 and is the amount of points to use to the left and right of the point to infer.



Now we can do the same thing as before but with Lanczos instead of sinc, taking care that the sinc within the lanczos function already multiplies $x$ by $\pi$:

$$
y =
0 \times \text{lanczos}\left(\frac{x - 0 \times 0.9}{0.9}, 2\right) +
0.7628 \times \text{lanczos}\left(\frac{x - 1 \times 0.9}{0.9}, 2\right) +
-0.4309 \times \text{lanczos}\left(\frac{x - 2 \times 0.9}{0.9}, 2\right) +
-0.3303 \times \text{lanczos}\left(\frac{x - 3 \times 0.9}{0.9}, 2\right) +
\dots
$$



The Lanczos filter can be used in things like image enlargement. By treating each pixel in an image as a sample of points from a continuous field, values in between two existing pixels can be deduced by interpolation by reproducing the continuous field using the Lanczos interpolation and then picking missing values from between two pixels.

In general, for samples $s$ and a sample period of $T$,

$$y = \sum_i s[i] \times sinc\left(\frac{\pi (x - i \times T)}{T}\right)$$
$$y = \sum_i s[i] \times lanczos\left(\frac{x - i \times T}{T}\right)$$