Fourier's fit coefficients

Question

lately i am been working fitting a fourier series function to a periodic signal for retrieve the amplitude and the phase of each component via least squares, so i modified the code of this file for it:

import math
import numpy as np
#period of the signal
per=1.0
w = 2.0*np.pi/per
#number of fourier components.
nf = 5
fp = open("file.cat","r")
# m1 is the number of unknown coefficients.
m1 = 2*nf + 1
# Create empty matrices.
x = np.zeros((m1,m1))
y = np.zeros((m1,1))
xi = [0.0]*m1

# Read (time, value) from each line of the file.
for line in fp:
    t = float(line.split()[0])
    yi = float(line.split()[1])
    xi[0] = 1.0
    for k in range(1,nf+1):
        xi[2*k-1] = np.sin(k*w*t)
        xi[2*k] = np.cos(k*w*t)
    for j in range(m1):
        for k in range(m1):
            x[j,k] += xi[j]*xi[k]
        y[j] += yi*xi[j]
fp.close()
# Copy to big matrices.
X = np.mat( x.copy() )
Y = np.mat( y.copy() )
# Invert X and multiply by Y to get coefficients.
A = X.I*Y
A0 = A[0]
# Solution is A0 + Sum[ Amp*sin(k*wt + phi) ]
print "a[0] = %f" % A[0]
for k in range(1,nf+1):
    amp = math.sqrt(A[2*k-1]**2 + A[2*k]**2)
    phs = math.atan2(A[2*k],A[2*k-1])
    print "amp[%d] = %f phi = %f" % (k, amp, phs)

but the plot show this (without the points, of course):

and it should show something like this:

somebody can tell me how can i compute the phase and the amplitude in another simpler way? a guide maybe, i will be very grateful. cheers!

PD. I will attach the FILE that i used, just because :)

EDITED

The error was with a index :(

First, I defined the vector with the values:

amp = np.array([np.sqrt((A[2*k-1])**2 + (A[2*k])**2) for k in range(1,nf+1)])
phs = np.array([math.atan2(A[2*k],A[2*k-1]) for k in range(1,nf+1)])

and then, to build the signal, I defined:

def term(t): return np.array([amp[k]*np.sin(k*w*t + phs[k]) for k in range(len(amp))])
Signal = np.array([A0+sum(term(phase[i])) for i in range(len(mag))])

but within the np.sin(), k should be k+1, because the index start in 0 ·__·

def term(t): return np.array([amp[k]*np.sin((k+1)*w*t + phs[k]) for k in range(len(amp))])
plt.plot(phase,Signal,'r-',lw=3)

and that is all.

Thanks Marco Tompitak for the help!!

Answer 1

You're specifying the wrong period for the signal:

#period of the signal
per=0.178556

This gives you the resulting Fourier fit, indeed with a maximum period of ~0.17. The problem is that this number specifies the longest period that is present in your Fourier series. The function only has components with perior 0.17 or shorter. Apparently you are expecting a fit with period ~1, so it can never approximate that properly. You should specify per=1.0. There's nothing wrong with the algorithm; a quick writeup of a similar algorithm in Mathematica gives the same output and plausible results: