In [1]:

    

# Enable autoreloading of modules
%load_ext autoreload
%autoreload 2


    

In [2]:

    

import numpy as np
import scipy.optimize
import pandas as pd

import holoviews as hv
from holoviews import opts
hv.extension('bokeh', 'matplotlib')


    

In [3]:

    

# Set default figure size on curves
opts.defaults(opts.Curve( height=400, width=900 ,show_grid=True))


    

In [18]:

    

def fit_sin(tt, yy):
    """
    Fit sin to the input time sequence, and return fitting parameters "amp",
    "omega", "phase", "offset", "freq", "period" and "fitfunc"
    
    @param tt: independent variable, raw. List/array.
    @param yy: dependent variable, raw. List/array.
    @return: Struct of results. Dict.
    
    Source: https://stackoverflow.com/questions/16716302/how-do-i-fit-a-sine-curve-to-my-data-with-pylab-and-numpy
    
    """
    
    tt = np.array(tt)
    yy = np.array(yy)
    ff = np.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = np.fft.fft(yy)
    Fyy_phase = np.angle(Fyy)
    Fyy_mag = abs(Fyy)
    peak_index = np.argmax(Fyy_mag[1:])+1
    guess_freq = abs(ff[peak_index])  # exclude the zero freq "peak"
    guess_phase = Fyy_phase[peak_index]
    guess_amp = np.std(yy) * 2.**0.5
    guess_offset = np.mean(yy)
    guess = np.array([guess_amp, 2.*np.pi*guess_freq, guess_phase, guess_offset])
    
    def sinfunc(t, A, w, p, c):  return A * np.sin(w*t + p) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, c = popt
    f = w/(2.*np.pi)
    fitfunc = lambda t: A * np.sin(w*t + p) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "freq": f,
            "period": 1./f, "fitfunc": fitfunc, "maxcov": np.max(pcov),
            "rawres": (guess,popt,pcov)}


    

In [19]:

    

f1 = 2.5
phi1 = 0
A1 = 1
offs1 = 0

fs = 100
dt = 1/fs
N = 512

t = np.arange(0,N)*dt
y1 = A1*np.sin(2*np.pi*f1*t + phi1) + offs1

pFunc1 = hv.Curve((t,y1),('t','Time'),('y','amplitude'),label='Func1')

fit1 = fit_sin(t,y1)
pfit1 = hv.Curve((t,fit1['fitfunc'](t)),('t','Time'),('y','amplitude'),label='Fit1')

(pfit1*pFunc1).opts(legend_position='right',title='Fitting clean sine wave with no phase offset')


    

    
Out[19]:

In [20]:

    

f1 = 2.5
phi1 =  180*np.pi/180
A1 = 1
offs1 = 0

fs = 100
dt = 1/fs
N = 512

t = np.arange(0,N)*dt
y1 = A1*np.sin(2*np.pi*f1*t + phi1) + offs1

pFunc1 = hv.Curve((t,y1),('t','Time'),('y','amplitude'),label='Func1')

fit1 = fit_sin(t,y1)
pfit1 = hv.Curve((t,fit1['fitfunc'](t)),('t','Time'),('y','amplitude'),label='Fit1')

(pfit1*pFunc1).opts(legend_position='right',title='Fitting clean sine wave with phase offset')


    

    
Out[20]:

In [22]:

    

f1 = 2.5
phi1 = 35
A1 = 1
offs1 = 0.7

fs = 100
dt = 1/fs
N = 512

t = np.arange(0,N)*dt
y1 = A1*np.sin(2*np.pi*f1*t + phi1) + offs1

pFunc1 = hv.Curve((t,y1),('t','Time'),('y','amplitude'),label='Func1')

fit1 = fit_sin(t,y1)
pfit1 = hv.Curve((t,fit1['fitfunc'](t)),('t','Time'),('y','amplitude'),label='Fit1')

(pfit1*pFunc1).opts(legend_position='right',title='Fitting clean sine wave with phase and amplitude offset')


    

    
Out[22]:

In [23]:

    

f1 = 2.5
phi1 = 35
A1 = 1
offs1 = 0.7

fs = 100
dt = 1/fs
N = 512

Amp = A1 + np.random.rand(N)*A1*0.2

t = np.arange(0,N)*dt
y1 = Amp*np.sin(2*np.pi*f1*t + phi1) + offs1

pFunc1 = hv.Curve((t,y1),('t','Time'),('y','amplitude'),label='Func1')

fit1 = fit_sin(t,y1)
pfit1 = hv.Curve((t,fit1['fitfunc'](t)),('t','Time'),('y','amplitude'),label='Fit1')

(pfit1*pFunc1).opts(legend_position='right',title='Fitting sine wave with phase & amplitude offset, amplitude noise')


    

    
Out[23]:

In [24]:

    

f1 = 2.5
phi1 = 35
A1 = 1
offs1 = 0.7

fs = 100
dt = 1/fs
N = 512

Amp = A1 + np.random.rand(N)*A1*0.2
offset = offs1 + np.random.rand(N)*offs1*0.2

t = np.arange(0,N)*dt
y1 = Amp*np.sin(2*np.pi*f1*t + phi1) + offset

pFunc1 = hv.Curve((t,y1),('t','Time'),('y','amplitude'),label='Func1')

fit1 = fit_sin(t,y1)
pfit1 = hv.Curve((t,fit1['fitfunc'](t)),('t','Time'),('y','amplitude'),label='Fit1')

(pfit1*pFunc1).opts(legend_position='right',title='Fitting sine wave with phase & amplitude offset, amplitude & offset noise ')


    

    
Out[24]:

In [22]:

    

f1 = 2.5
phi1 = 180*np.pi/180
A1 = 1
offs1 = 0.7

fs = 100
dt = 1/fs
N = 512

Amp = A1 + np.random.rand(N)*A1*0.4
offset = offs1 + np.random.rand(N)*offs1*0.2
ph_offset = phi1 + np.random.rand(N)*phi1*0.3

t = np.arange(0,N)*dt
y1 = Amp*np.sin(2*np.pi*f1*t + ph_offset) + offset

pFunc1 = hv.Curve((t,y1),('t','Time'),('y','amplitude'),label='Func1')

fit1 = fit_sin(t,y1)
pfit1 = hv.Curve((t,fit1['fitfunc'](t)),('t','Time'),('y','amplitude'),label='Fit1')

(pfit1*pFunc1).opts(legend_position='right',title='Fitting sine wave with phase & amplitude offset, amplitude, offset & phase noise ')


    

    
Out[22]:

In [ ]: