# Efficient sinc interpolation of a list

I'm looking for an efficient way to construct a sinc interpolation i.e. an interpolating function of a list of numbers in the form of a Fourier series with a chosen maximum frequency. The input list is supposed to be an equidistant sample of some real smooth function (e.g. a function with a finite-support Fourier transform) and the interpolation should be a continuous function that is a sum of harmonic functions (cosine and sine functions) such that their coefficients are equal to the Fourier coefficients of the list for all frequencies up to the chosen max frequency, (see Wikipedia: Whittaker–Shannon interpolation formula). Increasing the max frequency beyond the band-limit of the original function should give a perfect reconstruction of it. I'd guess there is a predefined Mathematica function doing that but I couldn't find it.

One way would be to apply Interpolation on the list, then NFourierSeries on the resulting InterpolationFunction (choosing the desired max frequency), but this is obviously very inefficient. LowpassFilter should be doing something similar but the output is a list, not a function. One option for an efficient solution that I'm trying to implement is to compute the Fourier coefficients from the list using Fourier, construct a function that is a list of harmonic functions with the right frequencies (e.g. Table[Sin[2 Pi n #/T], {n, nmax}] & for sines) and then Dot-multiply the two lists. I found this excellent post that computes the correct factors: Numerical Fourier transform of a complicated function

• Have you seen this and this? Commented Jan 10, 2021 at 2:07
• Maybe this, too. Commented Jan 10, 2021 at 3:12
• @J.M.'sennui Thanks! I had seen the second one but not the first. However, these approaches wouldn't be as efficient as I need, for example they use Sum which is slow for large lists. Commented Jan 10, 2021 at 12:06
• Sinc interpolation is interpolation using the assumption that the spectrum of the function is zero above the sampling rate of the function. So if you sinc-interpolate, then do NFourierSeries, you should find that all the additional terms are zero. An efficient way of doing sinc-interpolation is to Fourier transform, zero-pad the spectrum, then inverse Fourier transform. Commented Jan 10, 2021 at 13:03
• @MichaelE2 Thanks! I should read this more carefully. It seems to go a bit deeper into the efficiency aspect. It should also be possible to make it run faster using Compile. Commented Jan 10, 2021 at 19:27

A Sinc interpolation can be done along the same same idea as a Lagrange interpolation: Assuming we have equally spaced data points, we seek a function ipol that is 1 at x==0 and zero at every other data argument. The sum of ipol[x-x[[i]]] dat[[i]] over i gives then an interpolation function. By this the highest frequency is given by half the sample rate.

To demonstrate a Sinc interpolation, we first need some arbitrary test data, assuming for simplicity, that the x values are 1,2,..:

n = 20;
dat = Table[Exp[-x^2] Sin[ 4 x], {x, -Pi, Pi, 2 Pi/n}];
ListLinePlot[dat, PlotRange -> All]


Then we define a scaled Sinc that is 1 at x==0 and zero at every other integer value:

ipol[x_] = Sinc[2 Pi  x/2];
Plot[ipol[x], {x, -5, 5}, PlotRange -> All]


Now, we place such an function at every data point and multiply it with the value of the data point. Then we plot it together with the original data points :

fun[x_] = Table[ipol[x - y], {y, Length[dat]}].dat;
Plot[fun[x], {x, 0, 20}, PlotRange -> All,
Epilog -> Point[Transpose@{Range[n + 1], dat}]]


• Sinc[2 Pi x/2] is just Sinc[Pi x], no? Commented Jan 10, 2021 at 11:24
• Thanks! I need to apply this operation on ~10000 lists of length ~1000 each, which is why I need an efficient implementation, but this simple solution may actually be sufficient for my purposes. I would first apply a LowpassFilter on each list with a reasonable max frequency and then construct the function list Table[ipol[x - y], {y, Length[dat]}] only once and Dot-multiply it with each of the filtered lists. Commented Jan 10, 2021 at 14:03
• For some reasons, I would still prefer it if the output was a Fourier series instead of a sum of sinc-functions (for example I would like to be able to construct periodic interpolations as with the option Padding -> "Periodic" of Interpolation, which would be simpler using sines and cosines). But of course your solution gives the true sinc-interpolation, the truncated Fourier series is not precisely the same and it's not a true interpolation either. Commented Jan 10, 2021 at 14:04
• @Daniel Lichtblau. Hi Daniel, I think if you truncate the Fourier transform, the resulting back transform goes no more through the original points. You will then have a smoothing operation and no interpolation. Commented Jan 10, 2021 at 20:56
• @DanielHuber Actually my post wasn't so clear on this: I'm interested in both the smoothing operation and the interpolation as a special case (when the number of frequencies kept is equal to the size of the original list). Commented Jan 11, 2021 at 0:05

If we want to use FFT to get an expansion in circular functions we can proceed as follows:

FFT can be look at as a base change with the new base functions: Exp[2Pi I x f]. This is not too simple, because the FFT returns the coefficients like: First DC component, then coefficients of functions with increasing frequencies and pos. exponents. Up to the highest frequency (there may be one or two). Then down again in frequency with negative exponents.

Again we need same test data:

n = 20;
dat = Table[Exp[-x^2] Sin[4 x], {x, -Pi, Pi, 2 Pi/(n - 1)}];
ListLinePlot[dat, PlotRange -> All]


Then we calculate the FFT and if we want to smooth, truncate it:

fft = Fourier[dat];
smooth = 0;
If[smooth > 0,
If[EvenQ[n], fft[[2 + n/2 - smooth ;; n/2 + smooth]] = 0,
fft[[(n + 1)/2 - smooth ;; (n + 1)/2 + smooth]] = 0]];


We now have the coefficients of the new base and can assemble the searched for function and plot the result:

(*fourier synthese*)
tab = Table[
Exp[ -2 Pi I If[i < (n + 1)/2, (i - 1), -(n + 1 - i)] (x - 1)/
n], {i, n}] ;
fun[x_] = (tab.fft +
If[EvenQ[n],
fft[[1 + n/2]] (E^(-I \[Pi] (-1 + x)) - E^(I \[Pi] (-1 + x)))/2,
0])/Sqrt[n];
Plot[fun[x], {x, 0, n}, PlotRange -> All,
Epilog -> Point[Transpose[{Range[n], dat}]]]


No truncation:

With truncation (smooth=4):

As you see, truncation is not without drawbacks, it increases the wiggles. This is because use a rectangular window in the frequency domain. To improve, we would use a smoother window like Hamming, von Hann, ...

The function fun is now in exponential form. If you want it expanded in sin and cos:

fun[x_] = fun[x] // ComplexExpand // Chop