1
$\begingroup$

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)$$

Reading the data in:

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:

enter image description here

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[0],pfit[1],pfit[2],pfit[3],pfit[4],pfit[5],pfit[6],pfit[7],pfit[8],pfit[8],pfit[10],pfit[11],pfit[12],pfit[13],pfit[14],pfit[15],pfit[16],pfit[17],pfit[18],pfit[19],pfit[20],pfit[21])
plt.scatter(time,model,marker='+')
plt.plot(time,signal,c='r')
plt.show()

Yielding:

enter image description here

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):

enter image description here

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:

enter image description here

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.

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

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: enter image description here

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[0] )

        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] = 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()

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – zabop Feb 5 at 11:33
  • 1
    $\begingroup$ @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. $\endgroup$ – Cedron Dawg Feb 5 at 14:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.