# 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, 2016 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, 2016 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? Jun 24, 2016 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. Jun 24, 2016 at 17:17
• Similar ideas at dsp.stackexchange.com/questions/72433. Nov 29, 2021 at 18:15

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. Mar 11, 2017 at 23:15
• @robertbristow-johnson I've derived the result more rigorously and could confirm the expressions from Olli. Mar 12, 2017 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. Mar 12, 2017 at 2:20
• oh, i see you did. Mar 12, 2017 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. Mar 12, 2017 at 2:51

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. Mar 12, 2017 at 0:01
• Olli, i think Max is getting the booty. Mar 12, 2017 at 2:28
• Good, since he gave the derivations too. Mar 12, 2017 at 5:31

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? Jun 24, 2016 at 17:49

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.

• And let's also assume that $x[n]$, $y[n]$, and $y(t)$ are all real. Mar 11, 2017 at 4:55
• And, I think it won't kill generality to assume $T=1$. Mar 11, 2017 at 4:57
• Oh yes, back in 2001 on comp.dsp... groups.google.com/forum/#!msg/comp.dsp/RKY-gXf_c2U/jS2eq20-QdwJ Mar 12, 2017 at 16:36
• How is $y(t)$ a recovery of $y[n]$? Is it missing a windowing to make $y(t)=0$ outside $y[n]$'s support? Jan 9 at 15:46
• Because @OverLordGoldDragon, $y(t)$ is just and only the bandlimited reconstruction of the $y[n]$ samples, however they are defined. If all of those $y[n]$ are equal to zero (outside of $0\le n<N$), it's not because $y(t)$ was multiplied by some window in the continuous-time domain. $y[n]$ might be zero, but $y(t)$ is not zero except at those integer values of $t=nT$. $y(t)$ is not time-limited. Jan 10 at 8:25