# How can I improve my fit of cosines to periodic data using Python?

I have a space-separated csv file containing a measurement. First column is the time of measurement, second column is the corresponding measured value, third column is the error. The file can be found here. I would like to fit the parameters $$a_i$$, $$f$$, $$\phi_n$$ of the function $$g$$ to the data, using Python:

$$g(x) = A_0 + \sum_{n=0}^N a_i \cos(2n\pi fx + \phi_n)$$

import numpy as np
data=np.genfromtxt('signal.data')

time=data[:,0]
signal=data[:,1]
signalerror=data[:,2]


Plot the data:

import matplotlib.pyplot as plt
plt.figure()
plt.plot(time,signal)
plt.scatter(time,signal,s=5)
plt.show()


Getting the result: Now lets calculate a preliminary frequency guess of the periodic signal:

from gatspy.periodic import LombScargleFast
dmag=0.000005
nyquist_factor=40

model = LombScargleFast().fit(time, signal, dmag)
periods, power = model.periodogram_auto(nyquist_factor)

model.optimizer.period_range=(0.2, 10)
period = model.best_period


We get the result: 0.5467448186001437

I define the function to fit as follows, for $$N=10$$:

def G(x, A_0,
A_1, phi_1,
A_2, phi_2,
A_3, phi_3,
A_4, phi_4,
A_5, phi_5,
A_6, phi_6,
A_7, phi_7,
A_8, phi_8,
A_9, phi_9,
A_10, phi_10,
freq):
return (A_0 + A_1 * np.sin(2 * np.pi * 1 * freq * x + phi_1) +
A_2 * np.sin(2 * np.pi * 2 * freq * x + phi_2) +
A_3 * np.sin(2 * np.pi * 3 * freq * x + phi_3) +
A_4 * np.sin(2 * np.pi * 4 * freq * x + phi_4) +
A_5 * np.sin(2 * np.pi * 5 * freq * x + phi_5) +
A_6 * np.sin(2 * np.pi * 6 * freq * x + phi_6) +
A_7 * np.sin(2 * np.pi * 7 * freq * x + phi_7) +
A_8 * np.sin(2 * np.pi * 8 * freq * x + phi_8) +
A_9 * np.sin(2 * np.pi * 9 * freq * x + phi_9) +
A_10 * np.sin(2 * np.pi * 10 * freq * x + phi_10))


Now we need a function which fits G:

def fitter(time, signal, signalerror, LSPfreq):

from scipy import optimize

pfit, pcov = optimize.curve_fit(lambda x, _A_0,
_A_1, _phi_1,
_A_2, _phi_2,
_A_3, _phi_3,
_A_4, _phi_4,
_A_5, _phi_5,
_A_6, _phi_6,
_A_7, _phi_7,
_A_8, _phi_8,
_A_9, _phi_9,
_A_10, _phi_10,
_freqfit:

G(x, _A_0, _A_1, _phi_1,
_A_2, _phi_2,
_A_3, _phi_3,
_A_4, _phi_4,
_A_5, _phi_5,
_A_6, _phi_6,
_A_7, _phi_7,
_A_8, _phi_8,
_A_9, _phi_9,
_A_10, _phi_10,
_freqfit),

time, signal, p0=[11,  2, 0, #p0 is the initial guess for numerical fitting
1, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
LSPfreq],

sigma=signalerror, absolute_sigma=True)

error = []  # DEFINE LIST TO CALC ERROR
for i in range(len(pfit)):
try:
error.append(np.absolute(pcov[i][i]) ** 0.5)  # CALCULATE SQUARE ROOT OF TRACE OF COVARIANCE MATRIX
except:
error.append(0.00)
perr_curvefit = np.array(error)

return pfit, perr_curvefit


Check what we got:

LSPfreq=1/period
pfit, perr_curvefit = fitter(time, signal, signalerror, LSPfreq)

plt.figure()
model=G(time,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit,pfit)
plt.scatter(time,model,marker='+')
plt.plot(time,signal,c='r')
plt.show()


Yielding: Which is clearly wrong. If I play with the initial guesses p0 in the definition of function fitter, I can get a slightly better result. Setting

p0=[11,  1, 0,
0.1, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
LSPfreq]


Gives us (zoomed in): Which is a bit better. High frequency components are still in, despite the amplitude of the high frequency compontents were guessed to be zero. The original p0 seems also more justified than the modified version based on visual inspection of the data.

I played around with different values for p0, and while changing p0 indeed changes the result, I do not get a line reasonably well fitted to the data.

Why does this model fitting method fail? How can I improve get a good fit?

The whole code can be found here.

This question was also posted here.

EDIT:

PyCharm gives a warning for the p0 part of the code: Expected type 'Union[None,int,float,complex]', got 'List[Union[int,float],Any]]' instead

which with I don't know how to deal with, but might be relevant.

• Please explain why the downvote? Upon explanation I can improve the Q. Feb 4 '19 at 11:47
• I can't explain the downvote. +1 to at least bring you back to 0.
– Peter K.
Feb 4 '19 at 18:05

Your equation is a little muddled. It'll be easier for folks here if you state it this way:

$$g(x) = A_0 + \sum_{k=1}^L a_k \cos( 2\pi f x k + \phi_k)$$

I can't explain the downvote.

What you are trying to do is find the coefficients and phases for a Fourier Series for your function. Why your optiminzing best fit doesn't work, I don't know and I don't care to look at. It will not work very well on your signal as a whole, but there should be a good solution for a segment of your signal that displays a repeating pattern.

The conventional way to find the parameters is to frame a whole number of cycles, take the DFT, and use the bin values to derive your parameters. By using a multiple number of cycles, the fundamental and the harmonics will then fall into the bins of that multiple.

This is the results of a ten term series fit onto a repeat pattern of five: If you take the number of term in your series up to the number of points (and use all the bins) you will get an exact fit, for that interval.

Here is the code:



import matplotlib.pyplot           as plt
import numpy                       as np

#=========================================================================
def main() :

#---- Obtain the Data

theData = np.genfromtxt( "Data.txt" )

theTime        = theData[:,0]
theSignal      = theData[:,1]
theSignalError = theData[:,2]

#---- Show it

#        plt.figure()
#        plt.plot( theTime, theSignal )
#        plt.scatter( theTime, theSignal, s=5 )
#        plt.show()

#---- Find the First Minimum and the Fifth Minimum

theFirst = FindNextMajorMin( theSignal, 4 )

theNext = theFirst

for theMinCount in range( 5 ) :
theNext = FindNextMajorMin( theSignal, theNext )

#---- Extract a Repeat Pattern of Five

theSegment = theSignal[theFirst:theNext]

N = len( theSegment )

#---- Take the Normalized Fourier Transform

theDFT = np.fft.rfft( theSegment ) / float( N )

print theDFT

print N

#---- Translate into OP Parameters

a = np.zeros( 10 )
p = np.zeros( 10 )

A0 = np.real( theDFT )

for i in range( 1, 10 ) :
a[i] = np.abs(   theDFT[5*i] ) * 2.0
p[i] = np.angle( theDFT[5*i] )

#---- Construct the Partial Sums

theSums = np.zeros( (N,10) )

theFactor = 5.0 * 2.0 * np.pi / float( N )

for n in range( N ) :
theSums[n,0] = A0

for n in range( N ) :
for i in range( 1, 10 ) :
v = a[i] * np.cos( theFactor * n * i + p[i] )
theSums[n,i] = theSums[n,i-1] + v

#---- Show the Values

plt.plot( theSegment )
plt.plot( theSums )
plt.show()

plt.plot( theSegment )
plt.plot( theSums[:,9] )
plt.show()

theDFT = 0

plt.plot( np.abs( theDFT ) )
plt.show()

#=========================================================================
def FindNextMajorMin( argSignal, argLook ) :

theLimit = len( argSignal ) - 6

for n in range( argLook + 6, theLimit ) :
s = argSignal[n]

theMinFlag = True

for d in range( 1, 6 ) :
if s > argSignal[n-d] or \
s > argSignal[n-d] :
theMinFlag = False
break

if theMinFlag : return n

return -1

#=========================================================================
main()


• Indeed works, as prescribed. If you want to add this answer to here: stackoverflow.com/questions/54487279/… that might be useful for folks over there as well, or I can add it with your permission. Thx, great answer. Feb 5 '19 at 11:33
• @zabop, It would probably be best for you to post a link to this answer. Repost If you feel like it, I'm not a member of plain StachExchange. Please notice that the rough frequency calculation I did is slightly off, but you still get good coefficients. You can tell by how the fit is a little skinny, but the middle fit shows they define the waveform, the timbre so to speak. You can synthesize them back at any frequency by varying your time parameter. If you are going to implement this somewhere, you may want to calculate the DFT bins individually since you aren't using all of them. Feb 5 '19 at 14:33