Fitting sinusoidal data in Python

Question

I am trying to fit experimental data

with a function of the form:

A * np.sin(w*t + p) * np.exp(-g*t) + c

However, the fitted curve (the line in the following image) is not accurate:

If I leave out the exponential decay part, it works and I get a sinus function that is not decaying:

The function that I use is from this thread:

def fit_sin(tt, yy):
    '''Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"'''
    tt = np.array(tt)
    yy = np.array(yy)
    ff = np.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(np.fft.fft(yy))
    guess_freq = abs(ff[np.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = np.std(yy) * 2.**0.5
    guess_offset = np.mean(yy)
    guess_phase = 0.
    guess_damping = 0.5
    guess = np.array([guess_amp, 2.*np.pi*guess_freq, guess_phase, guess_offset, guess_damping])

    def sinfunc(t, A, w, p, g, c):  return A * np.sin(w*t + p) * np.exp(-g*t) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, g, c = popt
    f = w/(2.*np.pi)
    fitfunc = lambda t: A * np.sin(w*t + p) * np.exp(-g*t) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "damping": g, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": np.max(pcov), "rawres": (guess,popt,pcov)}

res = fit_sin(x, y)
x_fit = np.linspace(np.min(x), np.max(x), len(x))

plt.plot(x, y, label='Data', linewidth=line_width)
plt.plot(x_fit, res["fitfunc"](x_fit), label='Fit Curve', linewidth=line_width)
plt.show()

I am not sure if I implemented the code incorrectly or if the function is not able to describe my data correctly. I appreciate your help!

You can load the txt file from here:

GitHub

and manipulate the data like this to compare it with the post:

file = 'A2320_Data.txt'
column = 17

data = np.loadtxt(file, float)

start = 270
end = 36000

time_scale = 3600

x = []
y = []
for i in range(len(data)):
    if start < data[i][0] < end:
        x.append(data[i][0]/time_scale)
        y.append(data[i][column])

x = np.array(x)
y = np.array(y)

plt.plot(x, y, label='Pyro Oscillations', linewidth=line_width)

Answer 1

Your fitted curve will look like this

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
def sinfunc(t, A, w, p, g, c):  return A * np.sin(w*t + p) * np.exp(-g*t) + c

tt = np.linspace(0, 10, 1000)
yy = sinfunc(tt, -1, 10, 2, 0.3, 2)
plt.plot(tt, yy)

g stretches the envelope horizontally, c moves the center vertically, w determines the oscillation frequency, A stretches the envelope vertically.

So it can't accurately model the data you have.

Also, you will not be able to reliably fit w, to determine the oscilation frequency it is better to try an FFT.

Of course you could adjust the function to look like your data by adding a few more parameters, e.g.

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
def sinfunc(t, A, w, p, g, c1, c2, c3):  return A * np.sin(w*t + p) * (np.exp(-g*t) - c1) + c2 * np.exp(-g*t) + c3

tt = np.linspace(0, 10, 1000)
yy = sinfunc(tt, -1, 20, 2, 0.5, 1, 1.5, 1)
plt.plot(tt, yy)

But you will still have to give a good guess for the frequency.