# Whittaker-Shannon ($\mathrm{sinc}$) interpolation for a finite number of samples

Given an infinite number of samples $(N)$, a higher (or lower) number of samples $(cN)$ can be derived using sinc interpolation followed by sampling. How can this be applied to finite length signals?

With $\mathrm{sinc}$ interpolation, one can derive a continuous-time signal as:

$$y(t) = \sum^{\infty}_{n=-\infty} y[n]\mathrm{sinc}\left({t\over T}-n\right)$$

• For a finite number of sample points, should (can) we consider the $x[n]$ in the picture as

$$y[n] = \begin{cases} x[n], & \text{if } n \in [0, N-1] \\0, & \text{otherwise} \end{cases}$$

• Or should $y[n]$ be considered as a periodic version of $x[n]$? (This link briefly addresses the same. The original stated form cannot be directly used with periodic signals)

$$y[n] = x\left[n\pmod N\right]$$

In the first consideration, outside the region $[0,\ N-1]$, if I understand correctly, the Gibb's phenomenon would result in a ringing effect. Would this completely invalidate any values predicted outside the non-zero region or is it only that the degree of inconsistency is high? (More specifically for points close to but just outside the boundary in the interpolated continuous-time signal)

I am interested to know whether the addition of zeros would pollute the input set of points during the interpolation stage.

• Finite length signals have infinite bandwidth, so cannot be sampled according to the Nyquist criterion. The fact that I've never yet seen an infinite duration signal is neither here nor there. ;-) – Peter K. Jun 24 '16 at 15:26
• The sinc weights tend to zero at times far away from the interpolated sample, so eventually you can just ignore them. The fact that digital audio works is proof that the error can be negligible. You may need to do some tests for your particular case, especially if the number of samples is small. – MBaz Jun 24 '16 at 17:11
• @Fat32 Yes that was one of my questions. The other one was what should I choose as $x[n]$ in the original sinc-interpolation equation in order to resample inside the original sampling region: (i) the samples padded with zeros on both sides infinitely (similar to what you have considered) or (ii) the samples taken repetitively as a periodic sequence? Should edit the question to make it clearer? – Television Jun 24 '16 at 17:16
• @Peter I am not sure if I understood what you mean, but my intention is not to sample the finite length signal directly, but sample it after interpolation which would possibly result in an infinite duration signal. – Television Jun 24 '16 at 17:17
• @MBaz If I was to apply the samples affixed with zeros on either side then this behavior of $sinc(x)$ functions would ensure that the reconstructed signal dies to zero eventually giving a large error (wrt to the original signal), so I agree with your point that the magnitude of error really depends on the number of samples provided and the original signal itself. How about the second question with regard to samples inside the sampling period (refer to my comment above) – Television Jun 24 '16 at 17:25

An alternative explanation of what the function $$g(u)=\sum_{m\in\mathbb{Z}}\operatorname{sinc}(u-mN)$$ is.

First, note that $$g(u)=\sum_{m\in\mathbb{Z}}\operatorname{sinc}(u-mN)=\text{sinc}(u) \ \circledast \ \sum_{m\in\mathbb{Z}}\delta(u-mN)$$.

Let's calculate the Fourier Transform of $$g(u)$$:

$$G(f) \triangleq \mathscr{F}\{g(u)\}=\operatorname{rect}(f)\cdot\frac{1}{N}\sum_{m\in\mathbb{Z}}\delta(f-\tfrac{m}{N})$$

Here, we have used the convolution theorem. As we see, the spectrum is a sampled version of the rect function. It's clear that it is discrete, since $$g(u)$$ is periodic. We can now distinguish two cases: $$N$$ odd and $$N$$ even.

For the $$N$$ odd case the result is clear:

\begin{align}G(f)&=\frac{1}{N}\sum_{m=-(N-1)/2}^{(N-1)/2}\delta(f-\tfrac{m}{N})\\ g(u)&=\frac{1}{N}\sum_{m=-(N-1)/2}^{(N-1)/2}e^{j2\pi \frac{mu}{N}} \end{align}

which is just a time-scaled Dirichlet kernel. Let $$t=2\pi\frac{u}{N}$$, then

$$g(\tfrac{N}{2\pi}t)=\sum_{m=-((N-1)/2)}^{(N-1)/2}e^{jmt}=\frac{1}{N}\frac{\sin((N-1)/2+1/2)t)}{\sin(t/2)}=\frac{\sin(tN/2)}{N\sin(t/2)}.$$

Now doing backsubstituting to $$u$$, we get to $$g(u)=\frac{\sin(\pi u)}{N\sin(\pi u/N)}$$

which confirms the result from Olli.

For the case of even $$N$$, the treatment is a bit more difficult: The Diracs appear exactly at the discontinuity of the $$\operatorname{rect}()$$ function. But, treating the $$\operatorname{rect}()$$ function as e.g. the limit of a raised cosine where the rolloff goes to zero, we can argue that the rect has value of $$\tfrac12$$ at the discontinuity. Hence, in frequency domain $$G(f)$$ can be expressed as

\begin{align}G(f)&=\sum_{m=-(N/2-1)}^{N/2-1}\delta(f- \tfrac{m}{N})+\tfrac{1}{2}(\delta(f-\tfrac{1}{2})+\delta(f+\tfrac{1}{2}))\\ g(u)&=\sum_{m=-(N/2-1)}^{N/2-1}e^{j2\pi\frac{mu}{N}}+\cos(2\pi\tfrac{1}{2}u) \end{align}

Doing the same substitution as above, we get to

\begin{align}g(u)&=\frac{1}{N}\frac{\sin(\pi\frac{N-1}{N}u)}{\sin(\pi u/N)}+\frac{\cos(\pi u)}{N}\\ &=\frac{1}{N}\frac{\sin(\pi u - \pi u/N)+\cos(\pi u)\sin(\pi u/N)}{\sin(\pi u/N)}\end{align}

Now, using the trigonometric identity $$\sin(\pi u - \pi u/N)=\sin(\pi u)\cos(\pi u/N)-\cos(\pi u)\sin(\pi u/N)$$, some parts in the numerator cancel and we get

\begin{align}g(u)&=\frac{\sin(\pi u)\cos(\pi u/N)}{N\sin(\pi u/N)}\\&=\frac{\sin(\pi u)}{N\tan(\pi u/N)}\end{align}

which again confirms the result from Olli.

• Max, you're on the right track, too. check out Olli's answer and, maybe, between the two of you, you will get to a consistent and rigorous result. – robert bristow-johnson Mar 11 '17 at 23:15
• @robertbristow-johnson I've derived the result more rigorously and could confirm the expressions from Olli. – Maximilian Matthé Mar 12 '17 at 0:01
• you might win yourself some rep, Max. use the geometric series: $$\sum\limits_{n=0}^{N-1} a^n \ = \ \frac{a^N -1}{a-1}$$ and a shitload of trig identities. – robert bristow-johnson Mar 12 '17 at 2:20
• oh, i see you did. – robert bristow-johnson Mar 12 '17 at 2:24
• it's unlikely i will change my mind, but i will give this 24 hours and then award the bounty to you, Max. i will make the case for the $\tfrac12$ weighting of the dirac impulses at the edges of $\operatorname{rect}$ function. but your intuition is correct. – robert bristow-johnson Mar 12 '17 at 2:51

For finite length signals (and finite time DSP), one can approximate Sinc interpolation by using a windowed Sinc interpolation kernel, with a window of finite length. The width and shape (Von Hann, etc.) of the window determine the quality of the interpolation approximation. The window widths commonly chosen are far narrower than the periodicity of the periodic-Sinc or Dirichlet kernel, so that difference does not matter.

Sinc interpolation is usually invalid not only outside the non-zero signal region, but near both edges as well, as pure Whittaker-Shannon reconstruction is only valid for infinite length strictly bandlimited signals, which the edges of a rectangular window do not approximate well.

• Thanks, that's what I was looking for. Adding what @MBaz pointed out above, the absolute error would depend on the distance from the boundaries of the sampling period. Just one clarification with what you have mentioned about the width of the window chosen: what relation would it bear with the initial number of samples provided? Also, where does the periodicity of the periodic-sinc kernel come into the picture when choosing the window length? – Television Jun 24 '16 at 17:49

The continuous function (of Robert's bounty) that can be used to interpolate any $N$-periodic signal by convolving with it any $N$ consequtive uniform samples of the signal:

$$g(u) = \sum_{m=-\infty}^{\infty} \operatorname{sinc}\left(u-mN\right)$$

could be called the "$N$-periodic sinc" (Fig. 1).

Figure 1: The $N$-periodic sinc $g(u)$ with $N=6.$

The $N$-periodic sinc can only consist of those zero-phase complex exponentials that are harmonics of frequency $\frac{2 \pi}{N}$ (which has period $N$) and which are at most of frequency $\pi$. Those complex exponentials are of equal amplitude except for when $N$ is even in which case there is a positive and a negative Nyquist frequency present. The amplitudes of the frequency $\pm\pi$ complex exponentials need to be halved.

$$\begin{array}{l}g(u)&= \left\{\begin{array}{ll}\displaystyle\sum_{m=-(N-1)/2}^{(N-1)/2}\frac{e^{2\pi m i u / N}}{N}&\text{if }N\text{ is odd,}\\ \displaystyle\sum_{m=-(N-2)/2}^{(N-2)/2}\frac{e^{2\pi m i u / N}}{N} + \displaystyle\frac{\cos(\pi u)}{N}&\text{if }N\text{ is even.}\end{array}\right.\\ &= \left\{\begin{array}{ll}\displaystyle\frac{1}{N}+\displaystyle\sum_{m=1}^{(N-1)/2}\frac{2\cos(2\pi m u / N)}{N}&\text{if }N\text{ is odd,}\\ \displaystyle\frac{1}{N}+\displaystyle\sum_{m=1}^{(N-2)/2}\frac{2\cos(2\pi m u / N)}{N} + \displaystyle\frac{\cos(\pi u)}{N}&\text{if }N\text{ is even.}\end{array}\right.\end{array}$$

I don't have the math skills to derive the above properly, but I have tested it enough to be sure that it's right.

Except for at $u=0$ (and also there if we take the limit), sinc equals $\frac{\sin(\pi u)}{\pi u}.$ Similarly, and again, I don't have a proper derivation, but examining what we need to divide $\sin(\pi u)$ by to get $g(u)$:

$$g(u) = \left\{\begin{array}{ll}1&\text{if }u/N\text{ is integer},\\ \displaystyle\frac{\sin(\pi u)}{N\sin(\pi u/N)}&\text{if (}u/N\text{ is not integer) and (}N\text{ is odd),}\\ \displaystyle\frac{\sin(\pi u)}{N\tan(\pi u/N)}&\text{if (}u/N\text{ is not integer) and (}N\text{ is even).}\end{array}\right.$$

Also the limit of the second and the third case is $1$ as $u$ approaches an integer multiple of $N.$

It is not surprising that:

$$\lim_{N\to\infty}g(u) = \operatorname{sinc}(u).$$

• My derivations below agree to your results, derived via the convolution theorem and back-transform to the time domain. – Maximilian Matthé Mar 12 '17 at 0:01
• Olli, i think Max is getting the booty. – robert bristow-johnson Mar 12 '17 at 2:28
• Good, since he gave the derivations too. – Olli Niemitalo Mar 12 '17 at 5:31

This is not an answer, but I am developing the question a bit more.

So, in any case, we have

$$y(t) = \sum_{n=-\infty}^{\infty} y[n] \, \operatorname{sinc}\left(\frac{t - nT}{T}\right)$$

where

$$\operatorname{sinc}(u) \triangleq \begin{cases} \frac{\sin(\pi u)}{\pi u}, & \text{if } u \ne 0 \\1, & \text{if } u = 0 \end{cases}$$

All terms are bandlimited to a maximum frequency of $\frac{1}{2T}$, so the summation is bandlimited to the same bandlimit. And, in any case, we have

$$y(t) \Bigg|_{t = nT} = y[n]$$

so the reconstruction works out exactly at the sampling instances.

In the zero extended case,

$$y[n] = \begin{cases} x[n], & \text{if } 0 \le n < N \\0, & \text{otherwise} \end{cases}$$

it's easy:

$$y(t) = \sum_{n=0}^{N-1} x[n] \, \operatorname{sinc}\left(\frac{t - nT}{T}\right)$$

But in the periodic case,

$$y[n+N] = y[n] \qquad \forall n$$

$$y[n] = x\left[n\pmod N\right]$$

what is $y(t)$?

\begin{align} y(t) &= \sum_{n=-\infty}^{\infty} y[n] \, \operatorname{sinc}\left(\frac{t - nT}{T}\right) \\ &= \sum_{m=-\infty}^{\infty} \sum_{n=0}^{N-1} y[n+mN] \, \operatorname{sinc}\left(\frac{t - (n+mN)T}{T}\right) \\ &= \sum_{m=-\infty}^{\infty} \sum_{n=0}^{N-1} y[n] \, \operatorname{sinc}\left(\frac{t - (n+mN)T}{T}\right) \\ &= \sum_{m=-\infty}^{\infty} \sum_{n=0}^{N-1} x[n] \, \operatorname{sinc}\left(\frac{t - (n+mN)T}{T}\right) \\ &= \sum_{n=0}^{N-1} x[n] \, \sum_{m=-\infty}^{\infty} \operatorname{sinc}\left(\frac{t-nT-mNT}{T}\right) \\ \end{align}

Substituting $u \triangleq t-nT$

$$y(t) = \sum_{n=0}^{N-1} x[n] \, g(t-nT)$$

where

$$g(u) = \sum_{m=-\infty}^{\infty} \operatorname{sinc}\left(\frac{u-mNT}{T}\right)$$

Clearly $g(u)$ is periodic with period $NT$.

$$g(u+NT) = g(u) \qquad \forall u$$

What is the closed-form expression for $g(u)$ in terms of $u$, $N$, and $T$?

That's what the bounty is for.