Fast Python implementation of Whittaker–Shannon / sinc / bandlimited interpolation.

sinc_interpolation.py

import finufft
import numpy as np
import scipy.fft as sfft
from numpy.typing import NDArray


def sinc_interpolation(x: NDArray, s: NDArray, u: NDArray) -> NDArray:
            """Whittaker–Shannon or sinc or bandlimited interpolation.

            Args:
                x (NDArray): signal to be interpolated, can be 1D or 2D
                s (NDArray): time points of x (*s* for *samples*) 
                u (NDArray): time points of y (*u* for *upsampled*)

            Returns:
                NDArray: interpolated signal at time points *u*

            Reference:
                This code is based on https://gist.github.com/endolith/1297227
                and the comments therein.

            TODO:
                * implement FFT based interpolation for speed up

            """
            sinc_ = np.sinc((u - s[:, None])/(s[1]-s[0]))

            return np.dot(x, sinc_)

# only needing len(u) and len(u)/len(s) but this way the signature is similar
# FFT version is much (!) faster
def sinc_interpolation_fft(x: np.ndarray, s: np.ndarray, u: np.ndarray) -> np.ndarray:
    """
    Fast Fourier Transform (FFT) based sinc or bandlimited interpolation.
    
    Args:
        x (np.ndarray): signal to be interpolated, can be 1D or 2D
        s (np.ndarray): time points of x (*s* for *samples*) 
        u (np.ndarray): time points of y (*u* for *upsampled*)
        
    Returns:
        np.ndarray: interpolated signal at time points *u*
    """
    num_output = len(u)

    # Compute the FFT of the input signal
    X = sfft.rfft(x)

    # Create a new array for the zero-padded frequency spectrum
    X_padded = np.zeros(num_output // 2 + 1, dtype=complex)

    # Copy the original frequency spectrum into the zero-padded array
    X_padded[:X.shape[0]] = X

    # Compute the inverse FFT of the zero-padded frequency spectrum
    x_interpolated = sfft.irfft(X_padded, n=num_output)

    return x_interpolated * (num_output / len(s))


def ndf_t_interpolate(x: np.ndarray, s: np.ndarray, u: np.ndarray) -> np.ndarray:
    """
    Interpolate uniformly sampled data (x at times s) to arbitrary u
    via direct evaluation of the inverse DFT (type-2 nonuniform DFT).

    Args:
        x (np.ndarray): signal to be interpolated, can be 1D or 2D
        s (np.ndarray): time points of x (*s* for *samples*) 
        u (np.ndarray): time points of y (*u* for *upsampled*)
        
    Returns:
        np.ndarray: interpolated signal at time points *u*
    """
    N = x.shape[-1]
    T = s[1] - s[0]
    
    # compute FFT and corresponding frequencies
    X = sfft.fft(x)
    freqs = sfft.fftfreq(N, T)
    
    # build phase matrix and compute inverse transform
    phase = np.exp(2j * np.pi * np.outer(u - s[0], freqs))
    y = phase @ X / N
    
    return y.real if np.isrealobj(x) else y


def finufft_interpolate(
    x: np.ndarray,
    s: np.ndarray,
    u: np.ndarray,
    eps: float = 1e-6
) -> np.ndarray:
    """
    Efficiently interpolate uniformly sampled data (x at times s) to arbitrary u
    via FINUFFT type-2 NUFFT.

    Note:
        FINUFFT type-2 NUFFT is from uniform to nonuniform points.
        With type-1 NUFFT, nonuniform to uniform can be achieved.

    Args:
        x (np.ndarray): signal to be interpolated, can be 1D or 2D
        s (np.ndarray): time points of x (*s* for *samples*) 
        u (np.ndarray): time points of y (*u* for *upsampled*)
        
    Returns:
        np.ndarray: interpolated signal at time points *u*
    """
    N = x.shape[-1]
    T = s[1] - s[0]

    # FFT for uniform-frequency coefficients (no fftshift needed)
    f = np.fft.fftshift(sfft.fft(x))

    # Map u into [-π, π) interval for FINUFFT
    u_scaled = 2 * np.pi * (u - s[0]) / (N * T)

    # FINUFFT type-2 NUFFT: from uniform modes f to nonuniform points u_scaled
    c = finufft.nufft1d2(u_scaled, f, eps=eps, isign=1)

    # Normalize and return real part if input was real
    result = c / N
    return result.real if np.isrealobj(x) else result

Quick comparison between time and frequency domain filtering on a 1d (random) signal.
Please note that both axes are scaled logarithmically.

Benchmark and graph created with perfplot.

The FFT method only allows uniform resampling, whereas the time-domain one allows for arbitrary steps in u no?

@nils-werner : in principle you are right. However, one can adapt the functions...