Spectral Interpolation In Python

Interpolation is a common technique used in data analysis and signal processing to estimate values between known data points. A popular method for interpolation is using polynomials, but this is actually a pretty poor approach. A much more powerful technique is spectral interpolation which uses Fourier transforms as its backbone. In this blog post, we will explore how to calculate spectral interpolating functions and how to implement that in Python.

Suppose we have a function fⱼ given on an equidistant grid

The Fourier coefficients of fⱼ can be defined as

where k∈[−π/h,π/h] . The inverse transform is then

As spectral interpolating function, we simply replace xⱼ in this expression by x. Once we have the interpolating function, we can do whatever we want: differentiate, integrate, whatever.

But in practice, for given fⱼ , it is not a good idea to solve the integral above. It is much easier to use a trick: use the Kronecker delta!

The Kronecker delta is defined as

If we substitute that in the first equation, we get the Fourier transform of the Kronecker delta:

If we wanted to interpolate the Kronecker, we would write

Let’s see how this looks like. We can use the sinc function from numpy. By doing so we don't have to care programmatically about what happens at x=0 . But note that np.sinc already includes a factor of π in its argument!

import numpy as np
import matplotlib.pyplot as plt

h = 0.5
x_j = np.arange(-3, 3+h, h)
f_j = np.zeros_like(x_j)
f_j[abs(x_j) < h/2] = 1

x = np.linspace(-3, 3, 201)
f = np.sinc(x/h) # Note: np.sinc already contains a factor of pi!

plt.plot(x_j, f_j, 'o')
plt.plot(x, f, lw=2)
plt.xlabel('x');

Well, this doesn’t look like a very good interpolant. It is so wiggly! That’s because the Kronecker delta is extremely “non-smooth”. For smoother functions, this approach gives much better results. Now, the Kronecker delta has the wonderful property, that you can express every discrete function in terms of Kronecker deltas:

And since the Fourier transform is linear, we can simply write equation 2 as

Then we apply the usual recipe to replace xⱼ by x and get as interpolating function

Let’s try this in Python and suppose we want to interpolate some Gaussian:

def spectral_interpolate(f_j, x_j, x, h):
    f = np.zeros_like(x)
    for x_n, f_n in zip(x_j, f_j):
        f += f_n * np.sinc((x - x_n)/h)
    return f

f_j = np.exp(-x_j**2)
f = spectral_interpolate(f_j, x_j, x, h)

plt.plot(x_j, f_j, 'o')
plt.plot(x, f);

Now, this looks much better, right? Let’s see what the real error is by comparing with the exact function that we interpolated:

error = np.abs(np.exp(-x**2) - f)

plt.semilogy(x, error)
plt.grid()

So the error is less than 10⁻⁵. That’s quite amazing for such a coarse grid!

If you want to know more about spectral methods, I highly recommend L. Trefethen’s book “Spectral Methods in MATLAB”. The derivation here also roughly follows that book.

In conclusion, spectral interpolation is a powerful technique for estimating values between known data points. The implementation of the Kronecker delta trick makes the calculation and implementation of spectral interpolating functions computationally more efficient.