Implementing Euler's Form of Trigonometric Interpolation

at the moment I am struggling to correctly implement Euler's from of trigonometric interpolation. To guide you through the work I have already done and information I have gathered, I will pair code with the mathematical definitions. The code will be written in python and will make use of numpy functions. Please note, that I won't use the Fast Fourier Transformation.

Firstly, all of the experiments will be performed on the following dataset (x,y):

[('0.0', '0.0'),
 ('0.6283185307179586', '0.6427876096865393'),
 ('1.2566370614359172', '0.984807753012208'),
 ('1.8849555921538759', '0.8660254037844387'),
 ('2.5132741228718345', '0.3420201433256689'),
 ('3.141592653589793', '-0.34202014332566866'),
 ('3.7699111843077517', '-0.8660254037844384'),
 ('4.39822971502571', '-0.9848077530122081'),
 ('5.026548245743669', '-0.6427876096865396'),
 ('5.654866776461628', '-2.4492935982947064e-16')]

we define the Discrete Fourier Transform

So my implementation for this function is as follows

import numpy as np #consider numpy as imported from here on

def F_n(Y):
    n = len(Y)
    Y_hat = []
    for k in range(len(Y)):
        transformed_k = 1/n * sum([y_l * np.exp(-2 * np.pi * 1j* k * l/n) for l, y_l in enumerate(Y) ])
        Y_hat.append(transformed_k)
    return Y_hat

So comparing the resulting coefficients with the ones using np.fft.fft it looks like the coefficients are almost correct. They only differ by being shifted by a digit, i.e.

# F_n(y)
['(-1.33907057366955e-17+0j)',
 '(0.14283712054380923-0.439607454395199j)',
 '(-0.048591754799448425+0.06688081278992913j)',
 '(-0.039133572999081954+0.028432205056635337j)',
 '(-0.036913281031968816+0.01199385205986717j)',
 '(-0.036397023426620205-2.0058074207055733e-17j)',
 '(-0.03691328103196878-0.011993852059867215j)',
 '(-0.03913357299908168-0.028432205056635646j)',
 '(-0.04859175479944824-0.06688081278992904j)',
 '(0.1428371205438091+0.439607454395199j)']
# np.fft.fft(y)
['(-1.1102230246251565e-16+0j)',
 '(1.428371205438092-4.39607454395199j)',
 '(-0.4859175479944836+0.6688081278992911j)',
 '(-0.3913357299908192+0.2843220505663533j)',
 '(-0.36913281031968803+0.11993852059867194j)',
 '(-0.36397023426620184-1.1102230246251565e-16j)',
 '(-0.36913281031968803-0.11993852059867194j)',
 '(-0.3913357299908196-0.2843220505663534j)',
 '(-0.4859175479944836-0.6688081278992911j)',
 '(1.4283712054380922+4.39607454395199j)']

Now I want to implement the trigonometric interpolation. Here is the definition

I use this implementation for the theorem

def trig_interpolation(Y_hat, x_range, depth=1000):
    n = len(Y_hat)

    get_summand = lambda c_j,l,x: c_j*np.exp(2 * np.pi * 1j * l*x)

    y_intp = []
    x_intp = list((i/depth)*x_range for i in range(depth))
    if n%2==0:
        K = n//2 
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K+1,K+1), Y_hat)]))
    else:
        K = n//2+1
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K,K+1), Y_hat)]))

    return x_intp, y_intp

x_range = max(x)-min(x)
x_intp, y_intp = trig_interpolation(np.fft.fft(y), x_range)

where x_range denotes the range of the x values in the original set of data points, depth denotes the resolution of the interpolation and in line 4 the get_summand function represents the term $$\exp(2\pi j x)$$ to make the summing process a little easier to read.

When running my code on the coefficients given by numpy's fft, I get

plt.plot(x_intp,np.real(y_intp))
plt.plot(x,y, 'o')

Though, the points are aligned with the interpolation using numpy's fft, I expect the curve to look differently, behaving closer to an actual sine curve.

Using my implementation calculating the Fourier coefficients gives me a curve which is to be expected wrong.

---

I am asking to point out my mistakes, such that I can implement the trigonometric interpolation correctly according to the mathematical characterization I presented.

There are many separate problems with what you have done.

The most important is that the discrete Fourier transform does NOT produce the order of frequencies that you require. The order is roughly (see https://numpy.org/doc/stable/reference/routines.fft.html ):

freq[0] .... positive frequencies ... negative frequencies (mapped from very positive ones)

So you have to shift your frequencies. There is a routine numpy.fft.fftshift but it doesn't quite shift in line with the numbering that you want, so I have used numpy.roll here instead to rotate the array.

You need to divide your FFT by n to get the coefficients in the form that you require.

Your x_range is wrong - it is based on an implied x_max, not the last element of the array.

K is wrong in one case.

You want x/x_range in your sum if x is not going between 0 and 1.

import numpy as np
import matplotlib.pyplot as plt

n = 10
x = np.linspace( 0.0, 2 * np.pi, n, endpoint=False )
y = np.sin( x )

def trig_interpolation(Y_hat, x_range, depth=1000):
    n = len(Y_hat)

    get_summand = lambda c_j, l, x: c_j * np.exp(2 * np.pi * 1j * l * x / x_range )

    x_intp = [(i / depth) * x_range for i in range(depth)]                     
    y_intp = []
    K = n // 2                                                       # Fix K for the odd/even cases
    if n%2==0:
        Y_hat = np.roll( Y_hat, K - 1 )                              # Rotate
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K+1,K+1), Y_hat)]))
    else:
        Y_hat = np.roll( Y_hat, K )                                  # Rotate
        for x in x_intp:
            y_intp.append(sum([get_summand(c_j,l,x) for l,c_j in zip(range(-K,K+1), Y_hat)]))
    return x_intp, y_intp


x_range = n * ( x[1] - x[0] )                                        # Fix the range: it ISN'T defined by last element of x
x_intp, y_intp = trig_interpolation( np.fft.fft(y) / n, x_range)     # Divide the FFT by n

plt.plot(x_intp,np.real(y_intp))
plt.plot(x,y, 'o')
plt.show()