Calculating Fourier coefficients in Python using fft and correctly matching the fundamental frequency

Question

I have a curve in xyz coordinates. I want to fourier expand this curve as:

Curve fourier expansion.

I would like to calculate the coefficients of the expansion x_c,m,x_s,m, ... by using an fft in python. I am having trouble understanding how to match the fundamental frequency to θ so that the fft returns the correct coefficients. Note that I do not wish to use any expansion, I want this expansion specifically.

Here is my current code:

import numpy as np
from scipy import interpolate
from scipy.fft import rfft, rfftfreq
from math import pi
import matplotlib.pyplot as plt

def fourier_series(a0,a,b,order,x):
    return a0 + np.sum([a[n]*np.cos(n*x) + b[n]*np.sin(n*x) for n in range(order)])

order = 25
theta = np.linspace(0,2*pi,100)
curve = [[np.cos(x), 2*np.sin(x), np.cos(2*x) + np.sin(2*x)] for x in theta] #this could be any closed curve.

xArr, yArr, zArr = np.transpose(curve)

# Calculate the fft
freq = rfftfreq(len(xArr))
freq_series_cos = [n/len(xArr) for n in range(1,order+1)]
freq_series_sin = [n/len(xArr) for n in range(order+1)]
curvesFourier = []

#interpolate the fft and calculate the right frequencies to match CurveXYZFourier
for x in [xArr, yArr, zArr]:
    xf = rfft(x)/len(x)
    fft_0 = pi*xf[0].real
    fft_cos = pi*xf.real/2  #find the cosine coefficients
    fft_sin = -pi*xf.imag/2  #find the sine coefficients
    
    fft_cos_interpolated = interpolate.CubicSpline(freq,fft_cos) #interpolate the coefficients to pick the right frequencies 
    fft_sin_interpolated = interpolate.CubicSpline(freq,fft_sin) #interpolate the coefficients to pick the right frequencies
    
    b  = fft_sin_interpolated(freq_series_sin) 
    a0 = fft_0 
    a  = fft_cos_interpolated(freq_series_cos)
    
    x_fourier = [fourier_series(a0,a,b,order, i) for i in theta]
    curvesFourier.append(x_fourier)


ax = plt.axes(projection='3d')

ax.plot3D(xArr, yArr, zArr, "k--")
ax.scatter3D(curvesFourier[0], curvesFourier[1], curvesFourier[2])
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()

I am interpolating the frequencies. However, I expect there is a way of matching the fundamental frequency to what is intended, so that this step is not required. I have also tried adapting the solution given by gg349 in https://stackoverflow.com/questions/4258106/how-to-calculate-a-fourier-series-in-numpy with no success.

Thanks for the help!

Answer 1

Turns out I was overcomplicating the issue. Here is a solution:

import numpy as np
from scipy import interpolate
from scipy.fft import rfft
from math import pi
import matplotlib.pyplot as plt
   
def fourier_series(a0,a,b,order,x):
    return a0 + np.sum([a[n-1]*np.cos(n*x) for n in range(1,order)]) +np.sum([b[n]*np.sin(n*x) for n in range(order)])

order = 50
theta = np.linspace(0,2*pi,100)
curve = [[1 + 0.02*np.cos(x)+ 0.03*np.cos(2*x), 1 + 0.02*np.sin(x), np.cos(2*x) + np.sin(2*x)] for x in theta] #any closed curve

xArr, yArr, zArr = np.transpose(curve)

'''
# Compute the arclength parameterization of the curve, only needed for interpolating the data
dL = np.sqrt(np.sum(np.diff(curve, axis=0)**2, axis=1))
L = np.concatenate(([0], np.cumsum(dL)))
L = L*2*pi/L[-1]

#To find orders above len(x)/2 you must interpolate the data and choose the sampling so that len(n_samples) >= 2*order
x_interpol = interpolate.CubicSpline(L,xArr)
y_interpol = interpolate.CubicSpline(L,yArr)
z_interpol = interpolate.CubicSpline(L,zArr)
n_samples = np.linspace(0,2*pi,1000)
'''

coeffs = []
curvesFourier = []
curvesFourier2 = []

for x in [xArr, yArr, zArr]:
    #xf=rfft(x(n_samples))/len(x(n_samples)) #use this instead if you wish to find any order
    xf=rfft(x)/len(x)                        #using the real array produces more accurate results but the maximum order is limited to len(x)/2
    fft_0 = xf[0].real
    fft_cos = 2*xf[1:].real  #find the cosine coefficients
    fft_sin = -2*xf.imag  #find the sine coefficients
    x_fourier = [fourier_series(fft_0,fft_cos,fft_sin,order, i) for i in theta]
    curvesFourier.append(x_fourier)
    
ax = plt.axes(projection='3d')

ax.plot3D(xArr, yArr, zArr, "k--")
ax.plot3D(curvesFourier[0], curvesFourier[1], curvesFourier[2],"g-")

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")

plt.show()