Trigonometric interpolation (quadratic Wikipedia algorithm, supports nonuniform input/output domains)

triginterp.py

import numpy as np
from .triginterp_inner import triginterp_inner


def triginterp(x, y, xi):
    """Trigonometric interpolation with zero-padding (lol, not periodic). Easy to change."""
    P = np.zeros_like(xi, dtype=np.float64)
    triginterp_inner(
        x.astype(np.float64), y.astype(np.float64), xi.astype(np.float64), P
    )
    P[(xi > x.max()) | (xi < x.min())] = 0
    return P

triginterp_inner.pyx

#!python
#cython: language_level=3
cimport cython
from libc.math cimport sin, pi

cdef double halfpi = pi / 2.0

# https://en.wikipedia.org/wiki/Trigonometric_interpolation#Implementation

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def triginterp_inner(double[:] x, double[:] y, double[:] xi, double[:] P):
    cdef Py_ssize_t N = len(x)
    cdef double Nd = <double> N
    cdef Py_ssize_t M = len(xi)
    # print(N, M, Nd)

    cdef double h = 2. / Nd
    cdef double scale = (x[1] - x[0]) / h
    cdef double halfpi_scaled = halfpi / scale

    cdef Py_ssize_t j, k
    cdef double dxix, tau_kj, xij
    for j in range(M):
        xij = xi[j]
        for k in range(N):
            dxix = halfpi_scaled * (xij - x[k])
            if dxix == 0.0:
                P[j] += y[k]
            else:
                # needs tan in denom in even case if not zero padding
                tau_kj = sin(Nd * dxix) / (Nd * sin(dxix))
                P[j] += y[k] * tau_kj