Computing Fourier Series to represent data points

I wish to compute a function (a Fourier series) that passes through some set of given points.

Similar to what is going on here https://gofigure.impara.ai/ , but I wish not to animate it. I merely want the function so that I can draw the shape myself. I have read lots of math stuff describing it and code that animates it, but I am struggling with my implementation.

My current code is as follows [should be able to run in a python notebook alone]

import matplotlib.pyplot as plt 
import scipy
import cmath 
import math
import numpy as np
from matplotlib import collections  as mc
import pylab as pl


midpoints = [[305.0, 244.5], [293.5, 237.0], [274.5, 367.5], [270.5, 373.5], [229.5, 391.0], [216.0, 396.0], [302.0, 269.0], [295.0, 271.0], [60.5, 146.5], [54.0, 153.0], [52.0, 167.0], [54.0, 153.0], [52.0, 167.0], [45.0, 178.0], [75.0, 76.5], [68.5, 98.5], [75.0, 76.5], [97.0, 58.5], [283.5, 357.5], [274.5, 367.5], [309.0, 255.0], [305.0, 244.5], [309.0, 255.0], [302.0, 269.0], [299.5, 291.5], [300.0, 297.0], [300.0, 297.0], [295.0, 309.5], [62.5, 105.0], [61.0, 118.5], [62.5, 105.0], [68.5, 98.5], [58.0, 139.0], [60.5, 146.5], [241.0, 111.5], [252.0, 124.5], [256.0, 132.5], [252.0, 124.5], [283.0, 356.0], [283.5, 357.5], [300.0, 290.5], [299.5, 291.5], [300.0, 290.5], [296.0, 280.5], [296.0, 280.5], [295.0, 271.0], [158.0, 387.0], [177.0, 396.5], [197.5, 402.0], [192.5, 403.0], [189.5, 400.5], [192.5, 403.0], [197.5, 402.0], [202.5, 401.0], [214.0, 395.5], [216.0, 396.0], [202.5, 401.0], [214.0, 395.5], [233.5, 375.0], [229.5, 391.0], [233.5, 375.0], [249.0, 372.5], [282.5, 340.0], [284.5, 328.0], [284.5, 328.0], [295.0, 309.5], [45.0, 178.0], [49.5, 189.0], [57.0, 198.0], [49.5, 189.0], [238.5, 108.5], [241.0, 111.5], [162.0, 57.5], [170.0, 59.0], [239.5, 204.5], [239.0, 200.0], [293.5, 237.0], [291.0, 227.5], [265.0, 229.5], [291.0, 227.5], [239.0, 189.5], [245.5, 178.0], [239.0, 189.5], [241.0, 193.5], [241.0, 193.5], [239.0, 200.0], [55.0, 119.0], [61.0, 118.5], [53.5, 134.0], [58.0, 139.0], [50.0, 129.0], [55.0, 119.0], [50.0, 129.0], [53.5, 134.0], [107.0, 46.0], [119.5, 50.5], [97.0, 54.0], [97.0, 58.5], [107.0, 46.0], [97.0, 54.0], [150.5, 377.0], [158.0, 387.0], [257.5, 367.5], [270.5, 373.5], [249.0, 372.5], [257.5, 367.5], [280.0, 349.5], [282.5, 340.0], [280.0, 349.5], [283.0, 356.0], [239.5, 90.0], [238.5, 98.0], [238.5, 108.5], [238.5, 98.0], [130.0, 49.0], [119.5, 50.5], [189.0, 65.0], [191.0, 64.5], [189.0, 65.0], [177.0, 62.5], [170.0, 59.0], [177.0, 62.5], [256.0, 132.5], [257.5, 139.5], [128.0, 361.5], [127.5, 360.0], [136.5, 382.5], [131.5, 378.5], [126.5, 370.0], [131.5, 378.5], [128.0, 361.5], [126.5, 370.0], [105.5, 343.5], [101.0, 324.5], [105.5, 343.5], [121.5, 347.5], [126.0, 353.0], [127.5, 360.0], [121.5, 347.5], [126.0, 353.0], [191.0, 64.5], [198.5, 72.0], [237.5, 83.5], [239.5, 90.0], [145.5, 49.0], [138.5, 49.0], [159.0, 57.0], [162.0, 57.5], [145.5, 49.0], [159.0, 57.0], [265.0, 229.5], [254.5, 220.0], [253.0, 216.5], [254.5, 220.0], [253.0, 216.5], [248.0, 208.5], [248.0, 208.5], [239.5, 204.5], [245.0, 173.5], [245.5, 178.0], [250.0, 158.0], [245.0, 173.5], [257.5, 139.5], [250.0, 158.0], [177.0, 396.5], [181.0, 395.5], [181.0, 395.5], [189.5, 400.5], [147.0, 377.0], [150.5, 377.0], [140.5, 381.5], [147.0, 377.0], [140.5, 381.5], [136.5, 382.5], [92.5, 313.5], [101.0, 324.5], [99.5, 290.0], [92.5, 313.5], [98.0, 271.0], [99.5, 290.0], [134.5, 47.5], [130.0, 49.0], [134.5, 47.5], [138.5, 49.0], [73.0, 222.5], [71.0, 214.5], [107.5, 246.0], [110.5, 248.0], [104.0, 266.5], [98.0, 271.0], [69.0, 209.5], [71.0, 214.5], [226.5, 87.0], [237.5, 83.5], [226.5, 87.0], [205.5, 94.5], [205.5, 94.5], [205.5, 77.5], [198.5, 72.0], [205.5, 77.5], [96.5, 244.5], [107.5, 246.0], [109.0, 265.5], [104.0, 266.5], [108.0, 263.5], [110.5, 248.0], [109.0, 265.5], [108.0, 263.5], [65.0, 205.5], [69.0, 209.5], [57.0, 198.0], [56.0, 201.5], [65.0, 205.5], [56.0, 201.5], [90.0, 240.0], [96.5, 244.5], [83.0, 224.0], [73.0, 222.5], [90.5, 231.5], [83.0, 224.0], [90.0, 240.0], [90.5, 231.5]]

x_list = [ p[0] for p in midpoints ]
y_list = [ p[1] for p in midpoints ]
complexmdpts = [ [p[0]+1j*p[1]] for p in midpoints ]

plt.scatter(x_list, y_list, s=50, marker="x", color='y') 

coefs = np.fft.fftshift(scipy.fft.fft(complexmdpts))
n = len(coefs)
print("coeffs[{}]:\n{}".format(n, coefs[:5]))

#todo: sort coeffs?

# function in terms of t to trace out curve
def f(t):
    ftx=0
    fty=0
    for i in range(-int(n/2), int(n/2)+1):
        ftx+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).real.tolist()[0]
        fty+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).imag.tolist()[0]
    return [ftx/n, fty/n]
    

lines = [] # store computed lines segments to approximate function

pft = f(0) # compute first point

t_list = np.linspace(0, 2*math.pi, n) # compute list of dts to use when drawing

for t in t_list: 
    cft = f(t)
#     print("f({}): {} \n".format(t, cft))
    lines.append([cft,pft])
    pft = cft

#draw f(t) approximation
lc = mc.LineCollection(lines)
fig, ax = pl.subplots()
ax.add_collection(lc)
ax.autoscale()
ax.margins(0.1)

And here is my output:

Points I wish to outline

My function approximation

I am not convinced that the fast fourier transform is used correctly. From what I read the fft is what I need, and I shift it because the scipy fft returns the array shifted, and I think that the rest of my code is right, assuming the coefficients are correct, which is why I am suspicious of the coefficients.

Is there a step between the transform and the coefficients that I am missing? or is my function evaluation given the coefficients incorrect? or am I missing something else?

There are actually a few issues that need to be resolved to get the expected outline. Let's go over each of those issues.

Fourier coefficient computation

Since you are computing the FFT of the 2D array complexmdpts (each element being an array of 1 complex value), the default behavior of scipy.fft is to compute the FFT along the last axis. In this case it means that you are in fact computing n FFTs of length 1 and given that FFTs of length 1 are identities, the entire computation returns the original array. One solution would be to specify axis=0 explicitly:

scipy.fft(complexmdpts,axis=0) 

You could alternatively flatten the 2D array into a 1D array given that has fundamentally only one dimension, but that would result in more changes to the rest of your code that is based on this 2D structure.

Coefficient shifting

There is a common confusion when trying to interpret the FFT coefficients. The coefficients in the higher indices correspond to negative frequencies after the wraparound that occurs above the Nyquist frequency. To visualize these coefficients on a plot where the negative frequencies appear at their intuitive location, the coefficients in the higher indices are swapped using fftshift with the coefficients in the lower indices. This way the coefficients corresponding to negative frequencies appear before the coefficients corresponding to positive frequencies.

In your specific case, the computations in f(t) associate negative frequencies (whenever i is negative, so is the argument to cmath.exp) with coefficients at negative indices. Since python array indexing conveniently uses elements counting back from the end of the array for negative indices, you are correctly using coefficients corresponding to negative frequencies when indexing with negative indices. All this to say that you do not need to swap the coefficients with fftshift, and the coefficients are obtained directly with:

coefs = scipy.fft(complexmdpts,axis=0)

Sampling domain

You've specified t to range from 0 to 2*math.pi. Given your implementation of f(t) the domain actually ranges from 0 to n (for example, when i==1 the argument to cmath.exp goes from 0 to 1j*2*math.pi as t goes from 0 to n). To get the curve to span the entire domain, you should update your values of t accordingly with:

t_list = np.linspace(0, n, n)

Points ordering

Finally, by plotting the original series of points using a scatter plot you are hiding what a line plot would look like:

Trying to get a smooth function to interpolating this result in a plot with similar lines jumping from one side of the outline to the other. If you want to capture these jumps, then I guess you're done. But I suspect you likely want to capture the outline of the area. In that case you'd have to sort your original points so that they follow the outline. One approach would be to iteratively append the closest point to the previously selected point. This will give something that more closely resemble the contour:

For demonstration purposes (i.e. no claim that this is the best way to do it), you could do this sorting with something like the following:

def sort_points(points):
  # pick a point
  reference_point = points[0]
  sorted = [reference_point]
  remaining_points = range(1,len(points))
  for i in range(1,len(points)):

    # find the closest point to reference_point, 
    mindiff = np.sum(np.square(np.array(points[remaining_points[0]])-reference_point))
    idx = 0
    # loop over all the other remaining points
    for j in range(1,len(remaining_points)):
      diff = np.sum(np.square(np.array(points[remaining_points[j]])-reference_point))
      if diff < mindiff:
        mindiff = diff
        idx = j        
    # found the closest: update the selected point, and add it to the list of sorted points
    reference_point = points[remaining_points[idx]]
    sorted.append(reference_point )
    remaining_points = np.delete(remaining_points, idx)
  return sorted    

For reference, here's the complete code:

import matplotlib.pyplot as plt 
import scipy
import cmath 
import math
import numpy as np
from matplotlib import collections  as mc
import pylab as pl

def sort_points(points):
  # pick a point
  reference_point = points[0]
  sorted = [reference_point]
  remaining_points = range(1,len(points))
  for i in range(1,len(points)):

    # find the closest point to reference_point, 
    mindiff = np.sum(np.square(np.array(points[remaining_points[0]])-reference_point))
    idx = 0
    # loop over all the other remaining points
    for j in range(1,len(remaining_points)):
      diff = np.sum(np.square(np.array(points[remaining_points[j]])-reference_point))
      if diff < mindiff:
        mindiff = diff
        idx = j        
    # found the closest: update the selected point, and add it to the list of sorted points
    reference_point = points[remaining_points[idx]]
    sorted.append(reference_point )
    remaining_points = np.delete(remaining_points, idx)
  return sorted    


midpoints = [[305.0, 244.5], [293.5, 237.0], [274.5, 367.5], [270.5, 373.5], [229.5, 391.0], [216.0, 396.0], [302.0, 269.0], [295.0, 271.0], [60.5, 146.5], [54.0, 153.0], [52.0, 167.0], [54.0, 153.0], [52.0, 167.0], [45.0, 178.0], [75.0, 76.5], [68.5, 98.5], [75.0, 76.5], [97.0, 58.5], [283.5, 357.5], [274.5, 367.5], [309.0, 255.0], [305.0, 244.5], [309.0, 255.0], [302.0, 269.0], [299.5, 291.5], [300.0, 297.0], [300.0, 297.0], [295.0, 309.5], [62.5, 105.0], [61.0, 118.5], [62.5, 105.0], [68.5, 98.5], [58.0, 139.0], [60.5, 146.5], [241.0, 111.5], [252.0, 124.5], [256.0, 132.5], [252.0, 124.5], [283.0, 356.0], [283.5, 357.5], [300.0, 290.5], [299.5, 291.5], [300.0, 290.5], [296.0, 280.5], [296.0, 280.5], [295.0, 271.0], [158.0, 387.0], [177.0, 396.5], [197.5, 402.0], [192.5, 403.0], [189.5, 400.5], [192.5, 403.0], [197.5, 402.0], [202.5, 401.0], [214.0, 395.5], [216.0, 396.0], [202.5, 401.0], [214.0, 395.5], [233.5, 375.0], [229.5, 391.0], [233.5, 375.0], [249.0, 372.5], [282.5, 340.0], [284.5, 328.0], [284.5, 328.0], [295.0, 309.5], [45.0, 178.0], [49.5, 189.0], [57.0, 198.0], [49.5, 189.0], [238.5, 108.5], [241.0, 111.5], [162.0, 57.5], [170.0, 59.0], [239.5, 204.5], [239.0, 200.0], [293.5, 237.0], [291.0, 227.5], [265.0, 229.5], [291.0, 227.5], [239.0, 189.5], [245.5, 178.0], [239.0, 189.5], [241.0, 193.5], [241.0, 193.5], [239.0, 200.0], [55.0, 119.0], [61.0, 118.5], [53.5, 134.0], [58.0, 139.0], [50.0, 129.0], [55.0, 119.0], [50.0, 129.0], [53.5, 134.0], [107.0, 46.0], [119.5, 50.5], [97.0, 54.0], [97.0, 58.5], [107.0, 46.0], [97.0, 54.0], [150.5, 377.0], [158.0, 387.0], [257.5, 367.5], [270.5, 373.5], [249.0, 372.5], [257.5, 367.5], [280.0, 349.5], [282.5, 340.0], [280.0, 349.5], [283.0, 356.0], [239.5, 90.0], [238.5, 98.0], [238.5, 108.5], [238.5, 98.0], [130.0, 49.0], [119.5, 50.5], [189.0, 65.0], [191.0, 64.5], [189.0, 65.0], [177.0, 62.5], [170.0, 59.0], [177.0, 62.5], [256.0, 132.5], [257.5, 139.5], [128.0, 361.5], [127.5, 360.0], [136.5, 382.5], [131.5, 378.5], [126.5, 370.0], [131.5, 378.5], [128.0, 361.5], [126.5, 370.0], [105.5, 343.5], [101.0, 324.5], [105.5, 343.5], [121.5, 347.5], [126.0, 353.0], [127.5, 360.0], [121.5, 347.5], [126.0, 353.0], [191.0, 64.5], [198.5, 72.0], [237.5, 83.5], [239.5, 90.0], [145.5, 49.0], [138.5, 49.0], [159.0, 57.0], [162.0, 57.5], [145.5, 49.0], [159.0, 57.0], [265.0, 229.5], [254.5, 220.0], [253.0, 216.5], [254.5, 220.0], [253.0, 216.5], [248.0, 208.5], [248.0, 208.5], [239.5, 204.5], [245.0, 173.5], [245.5, 178.0], [250.0, 158.0], [245.0, 173.5], [257.5, 139.5], [250.0, 158.0], [177.0, 396.5], [181.0, 395.5], [181.0, 395.5], [189.5, 400.5], [147.0, 377.0], [150.5, 377.0], [140.5, 381.5], [147.0, 377.0], [140.5, 381.5], [136.5, 382.5], [92.5, 313.5], [101.0, 324.5], [99.5, 290.0], [92.5, 313.5], [98.0, 271.0], [99.5, 290.0], [134.5, 47.5], [130.0, 49.0], [134.5, 47.5], [138.5, 49.0], [73.0, 222.5], [71.0, 214.5], [107.5, 246.0], [110.5, 248.0], [104.0, 266.5], [98.0, 271.0], [69.0, 209.5], [71.0, 214.5], [226.5, 87.0], [237.5, 83.5], [226.5, 87.0], [205.5, 94.5], [205.5, 94.5], [205.5, 77.5], [198.5, 72.0], [205.5, 77.5], [96.5, 244.5], [107.5, 246.0], [109.0, 265.5], [104.0, 266.5], [108.0, 263.5], [110.5, 248.0], [109.0, 265.5], [108.0, 263.5], [65.0, 205.5], [69.0, 209.5], [57.0, 198.0], [56.0, 201.5], [65.0, 205.5], [56.0, 201.5], [90.0, 240.0], [96.5, 244.5], [83.0, 224.0], [73.0, 222.5], [90.5, 231.5], [83.0, 224.0], [90.0, 240.0], [90.5, 231.5]]
midpoints = sort_points(midpoints)

x_list = [ p[0] for p in midpoints ]
y_list = [ p[1] for p in midpoints ]
complexmdpts = [ [p[0]+1j*p[1]] for p in midpoints ]

plt.scatter(x_list, y_list, s=50, marker="x", color='y') 

coefs = scipy.fft(complexmdpts,axis=0)
n = len(coefs)
print("coeffs[{}]:\n{}".format(n, coefs[:5]))

# function in terms of t to trace out curve
m = n
def f(t):
    ftx=0
    fty=0
    for i in range(-int(m/2), int(m/2)+1):
        ftx+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).real.tolist()[0]
        fty+=(coefs[i]*cmath.exp(1j*2*math.pi*i*t/n)).imag.tolist()[0]
    return [ftx/n, fty/n]
    

lines = [] # store computed lines segments to approximate function

pft = f(0) # compute first point

t_list = np.linspace(0, n, n) # compute list of dts to use when drawing

for t in t_list: 
    cft = f(t)
    lines.append([cft,pft])
    pft = cft

#draw f(t) approximation
lc = mc.LineCollection(lines)
fig, ax = pl.subplots()
ax.add_collection(lc)
ax.autoscale()
ax.margins(0.1)