Fourier series coefficients from fft vs least-squares fitting

gistfile1.txt

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
import abel


def cosine(r):
    return 10*np.cos(2*np.pi*2*r/r[-1]) +\
            5*np.cos(2*np.pi*8*r/r[-1]) +\
            2*np.cos(2*np.pi*11*r/r[-1])


def series(An, r):
    sum = np.zeros_like(r)
    for n, an in enumerate(An):
        sum += an*np.cos(2*np.pi*n*r/r[-1])
    return sum


def residual(An, r, signal):
    return signal - series(An, r)


# signal
r = np.linspace(0, 1, 201)
signal = cosine(r)

n = r.size

# Fourier transform --------------------
# duplicate function to make symmetric
signalx2 = np.append(signal[::-1], signal)

fourier_signal = np.fft.rfft(signalx2)/n
freq = np.fft.rfftfreq(signalx2.size, d=r[1]-r[0])

a = fourier_signal[0:-1].real

# Least-squares -------------------
# fitting Fourier series
Nu = 20
An = np.arange(Nu)

res = least_squares(residual, An, args=(r, signal))
An = res.x

# re-align fft coefficients
fA = np.zeros_like(An)
fA[0] = 0
print("coefficients")
print(" n    lsq-An   fft-An")
print("{:2d} {:8.2f} {:8.2f}".format(0, An[0], fA[0]))
for i in np.arange(1, Nu):
    fA[i] = a[np.abs(freq-i).argmin()]
    print("{:2d} {:8.2f} {:8.2f}".format(i, An[i], fA[i]))


# Plots ---------------
plt.plot(r, signal, label="signal")
plt.plot(r, series(An, r), label="series")
plt.plot(r, series(fA, r), label="fft")
plt.xlabel(r'$r$')
plt.ylabel(r'amplitude')
plt.title(r'$10\cos(2\pi\,2r)+5\cos(2\pi\,8r)+2\cos(2\pi\,11r)$')
plt.legend()

plt.savefig("fftvslsq.png", dpi=75)
plt.show()