How do I fit a sine curve to my data with pylab and numpy?

by

Here is a parameter-free fitting function fit_sin() that does not require manual guess of frequency:

import numpy, scipy.optimize

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 = numpy.array(tt)
    yy = numpy.array(yy)
    ff = numpy.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(numpy.fft.fft(yy))
    guess_freq = abs(ff[numpy.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = numpy.std(yy) * 2.**0.5
    guess_offset = numpy.mean(yy)
    guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])

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

The initial frequency guess is given by the peak frequency in the frequency domain using FFT. The fitting result is almost perfect assuming there is only one dominant frequency (other than the zero frequency peak).

import pylab as plt

N, amp, omega, phase, offset, noise = 500, 1., 2., .5, 4., 3
#N, amp, omega, phase, offset, noise = 50, 1., .4, .5, 4., .2
#N, amp, omega, phase, offset, noise = 200, 1., 20, .5, 4., 1
tt = numpy.linspace(0, 10, N)
tt2 = numpy.linspace(0, 10, 10*N)
yy = amp*numpy.sin(omega*tt + phase) + offset
yynoise = yy + noise*(numpy.random.random(len(tt))-0.5)

res = fit_sin(tt, yynoise)
print( "Amplitude=%(amp)s, Angular freq.=%(omega)s, phase=%(phase)s, offset=%(offset)s, Max. Cov.=%(maxcov)s" % res )

plt.plot(tt, yy, "-k", label="y", linewidth=2)
plt.plot(tt, yynoise, "ok", label="y with noise")
plt.plot(tt2, res["fitfunc"](tt2), "r-", label="y fit curve", linewidth=2)
plt.legend(loc="best")
plt.show()

The result is good even with high noise:

Amplitude=1.00660540618, Angular freq.=2.03370472482, phase=0.360276844224, offset=3.95747467506, Max. Cov.=0.0122923578658

With noise
Low frequency
High frequency

More Related Contents:

  1. Linear regression with matplotlib / numpy
  2. How to find the exact intersection of a curve (as np.array) with y==0?
  3. How to do exponential and logarithmic curve fitting in Python? I found only polynomial fitting
  4. FutureWarning: elementwise comparison failed; returning scalar, but in the future will perform elementwise comparison
  5. Matplotlib: save plot to numpy array
  6. How to apply piecewise linear fit in Python?
  7. find length of sequences of identical values in a numpy array (run length encoding)
  8. Matplotlib log scale tick label number formatting
  9. How can I attach a pyplot function to a figure instance?
  10. scipy: savefig without frames, axes, only content
  11. How to remove the space between subplots in matplotlib.pyplot?
  12. Plotting a decision boundary separating 2 classes using Matplotlib’s pyplot
  13. No module named when using PyInstaller
  14. How to plot empirical cdf (ecdf)
  15. How can I release memory after creating matplotlib figures
  16. Matplotlib: TypeError: can’t multiply sequence by non-int of type ‘numpy.float64’
  17. How to convert a .wav file to a spectrogram in python3
  18. figure of imshow() is too small
  19. How to ‘turn off’ blurry effect of imshow() in matplotlib?
  20. Creating a Confidence Ellipses in a sccatterplot using matplotlib
  21. How do I calculate r-squared using Python and Numpy?
  22. How to generate audio from a numpy array?
  23. How to overplot a line on a scatter plot in python?
  24. scipy.curve_fit() returns multiple lines
  25. Python: find contour lines from matplotlib.pyplot.contour()
  26. Change values on matplotlib imshow() graph axis
  27. setting spacing between grouped bar plots in matplotlib
  28. pylab.hist(data, normed=1). Normalization seems to work incorrect
  29. How to create a legend for 3D bar in matplotlib?
  30. 3D Plots over non-rectangular domain

Leave a Comment