In [1]:

    

from lmfit import Parameters
import matplotlib
import matplotlib.pyplot as plt
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.signal import wiener, filtfilt, butter, gaussian, freqz
from scipy.ndimage import filters
from core.util.units import compute_dft


    

In [2]:

    

def sine_testing():
    """ Sinus fit testing with a self defined estimator. """

    x_axis = np.linspace(0, 25, 75)
    x_axis1 = np.linspace(25, 50, 75)
    x_axis = np.append(x_axis, x_axis1)
    x_nice = np.linspace(x_axis[0],x_axis[-1], 1000)

    mod,params = fitlogic.make_sine_model()
    print('Parameters of the model',mod.param_names,
          ' with the independet variable',mod.independent_vars)

    print(1/(x_axis[1]-x_axis[0]))
    params['amplitude'].value=0.9 + np.random.normal(0,0.4)
    params['frequency'].value=0.1
    params['phase'].value=np.pi*1.0
    params['offset'].value=3.94+np.random.normal(0,0.4)

    data_noisy=(mod.eval(x=x_axis,params=params) + 0.3*np.random.normal(size=x_axis.shape))


    # set the offset as the average of the data
    offset = np.average(data_noisy)

    # level data
    data_level = data_noisy - offset

    # estimate amplitude
    ampl_val = max(np.abs(data_level.min()), np.abs(data_level.max()))

    # calculate dft with zeropadding to obtain nicer interpolation between the
    # appearing peaks.
    dft_x, dft_y = compute_dft(x_axis, data_level, zeropad_num=1)

    plt.figure()
    plt.plot(dft_x, dft_y, label='dft of data')
    plt.xlabel('Frequency')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
               ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    stepsize = x_axis[1]-x_axis[0]  # for frequency axis
    frequency_max = np.abs(dft_x[np.log(dft_y).argmax()])

    print(frequency_max)
    print('offset',offset)
#
    # find minimal distance to the next meas point in the corresponding time
    # value
    diff_array = np.ediff1d(x_axis)
    min_x_diff = diff_array.min()

    # if at least two identical values are in the array, then the difference is
    # of course zero, catch that case.
    for tries in range(len(diff_array)):
        if np.isclose(min_x_diff, 0.0):
            index = np.argmin(diff_array)
            diff_array = np.delete(diff_array, index)
            min_x_diff = diff_array.min()

        else:
            if len(diff_array) == 0:
                fitlogic.log.error('The passed x_axis for the sinus estimation contains the same values!'
                                   ' Cannot do the fit!')

                return -1, params
            else:
                min_x_diff = diff_array.min()
            break

    # How many points are used to sample the estimated frequency with min_x_diff:
    iter_steps = int(1/(frequency_max*min_x_diff))
    if iter_steps < 1:
        iter_steps = 1

    sum_res = np.zeros(iter_steps)

    # Procedure: Create sin waves with different phases and perform a summation.
    #            The sum shows how well the sine was fitting to the actual data.
    #            The best fitting sine should be a maximum of the summed time
    #            trace.

    for iter_s in range(iter_steps):
        func_val = ampl_val * np.sin(2*np.pi*frequency_max*x_axis + iter_s/iter_steps *2*np.pi)
        sum_res[iter_s] = np.abs(data_level - func_val).sum()

    plt.figure()
    plt.plot(list(range(iter_steps)), sum_res, label='sum of different phase iteration')
    plt.xlabel('iteration')
    plt.ylabel('summed value')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    # The minimum indicates where the sine function was fitting the worst,
    # therefore subtract pi. This will also ensure that the estimated phase will
    # be in the interval [-pi,pi].
    phase = sum_res.argmax()/iter_steps *2*np.pi - np.pi

    mod, params = fitlogic.make_sine_model()

    # values and bounds of initial parameters
    params['amplitude'].set(value=ampl_val)
    params['frequency'].set(value=frequency_max, min=0.0, max=1/(stepsize)*3)
    params['phase'].set(value=phase, min=-np.pi, max=np.pi)
    params['offset'].set(value=offset)

    # perform fit:
    result = mod.fit(data_noisy, x=x_axis, params=params)

    plt.plot(x_axis,data_noisy,'ob', label='noisy data')
    plt.plot(x_axis,result.init_fit,'-y', label='initial values')
    plt.plot(x_axis,result.best_fit,'-r',linewidth=2.0, label='actual fit')
    plt.xlabel('time')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)

    plt.show()


    

In [3]:

    

sine_testing()


    

    




Parameters of the model ['amplitude', 'frequency', 'phase', 'offset']  with the independet variable x
2.96
0.0986666666666
offset 3.47077663886






    













    













    

In [4]:

    

def sine_testing2():
    """ Sinus fit testing with the direct fit method. """

    x_axis = np.linspace(0, 250, 75)
    x_axis1 = np.linspace(250, 500, 75)
    x_axis = np.append(x_axis, x_axis1)
    x_nice = np.linspace(x_axis[0],x_axis[-1], 1000)

    mod, params = fitlogic.make_sine_model()

    params['phase'].value = np.pi/2 # np.random.uniform()*2*np.pi
    params['frequency'].value = 0.01
    params['amplitude'].value = 1.5
    params['offset'].value = 0.4

    data = mod.eval(x=x_axis, params=params)
    data_noisy = (mod.eval(x=x_axis, params=params) + 1.5 * np.random.normal(size=x_axis.shape))

    update_dict = {}
    update_dict['phase'] = {'vary': False, 'value': np.pi/2.}

    result = fitlogic.make_sine_fit(
        x_axis=x_axis,
        data=data_noisy,
        add_params=update_dict,
        estimator=fitlogic.estimate_sine)

    plt.figure()
    plt.plot(x_axis, data_noisy, label='noisy data')
    plt.plot(x_axis, result.init_fit, label='initial data')
    plt.plot(x_axis, result.best_fit, label='fit data')
    plt.xlabel('time')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()


    

In [5]:

    

sine_testing2()


    

In [6]:

    

def two_sine_offset_testing():
    """ Testing procedure for the estimator for a double sine with offset fit. """

    x_axis = np.linspace(5, 300 ,200)
    phase1 = np.random.uniform()*2*np.pi
    ampl1 = 1
    freq1 = 0.02

    phase2 = np.random.uniform()*2*np.pi
    ampl2 = 1
    freq2 = 0.01

    offset = 1

    data = ampl1 * np.sin(2*np.pi*freq1*x_axis +phase1) +ampl2 * np.sin(2*np.pi*freq2*x_axis +phase2) + offset

    noisy_data = data + data.mean() * np.random.normal(size=x_axis.shape)*2

    res = fitlogic.make_sine_fit(x_axis=x_axis, data=noisy_data, estimator=fitlogic.estimate_sine)

    data_sub = noisy_data - res.best_fit

    res2 = fitlogic.make_sine_fit(x_axis=x_axis, data=data_sub, estimator=fitlogic.estimate_sine)

    plt.figure()
    plt.plot(x_axis, res.best_fit+res2.best_fit,'-', label='fit')
    plt.plot(x_axis, noisy_data, 'o--', label='noisy_data')
    plt.plot(x_axis, data,'-', label='ideal data')
    plt.xlabel('Time micro-s')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    mod, params = fitlogic.make_sinedouble_model()

    params['s1_amplitude'].set(value=res.params['amplitude'].value)
    params['s1_frequency'].set(value=res.params['frequency'].value)
    params['s1_phase'].set(value=res.params['phase'].value)

    params['s2_amplitude'].set(value=res2.params['amplitude'].value)
    params['s2_frequency'].set(value=res2.params['frequency'].value)
    params['s2_phase'].set(value=res2.params['phase'].value)

    params['offset'].set(value=data.mean())

    result = mod.fit(noisy_data, x=x_axis, params=params)

    plt.figure()
    plt.plot(x_axis, result.best_fit,'-', label='fit')
    plt.plot(x_axis, noisy_data, 'o--', label='noisy_data')
    plt.plot(x_axis, data,'-', label='ideal data')
    plt.xlabel('Time micro-s')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    print("freq1:", result.params['s1_frequency'].value)
    print("freq2:", result.params['s2_frequency'].value)


    

In [7]:

    

two_sine_offset_testing()


    

    




freq1: 0.0106850987403
freq2: 0.0196366616392






    













    

In [8]:

    

def two_sine_offset_testing2():
    """ Testing procedure for the implemented two sine with exponential
        decay and offset fit. """

    x_axis = np.linspace(5, 300 ,200)
    phase1 = np.random.uniform()*2*np.pi
    ampl1 = 1
    freq1 = 0.02

    phase2 = np.random.uniform()*2*np.pi
    ampl2 = 1
    freq2 = 0.01

    offset = 1

    data = ampl1 * np.sin(2*np.pi*freq1*x_axis +phase1) +ampl2 * np.sin(2*np.pi*freq2*x_axis +phase2) + offset

    noisy_data = data + data.mean() * np.random.normal(size=x_axis.shape)*1.0

    result = fitlogic.make_sinedouble_fit(
        x_axis=x_axis,
        data=noisy_data,
        estimator=fitlogic.estimate_sinedouble)

    plt.figure()
    plt.plot(x_axis, result.best_fit,'-', label='fit')
    plt.plot(x_axis, noisy_data, 'o--', label='noisy_data')
    plt.plot(x_axis, data,'-', label='ideal data')
    plt.xlabel('Time micro-s')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    print(result.fit_report())


    

In [9]:

    

two_sine_offset_testing2()


    

    




[[Model]]
    (((Model(amplitude_function, prefix='s1_') * Model(bare_sine_function, prefix='s1_')) + (Model(amplitude_function, prefix='s2_') * Model(bare_sine_function, prefix='s2_'))) + Model(constant_function))
[[Fit Statistics]]
    # function evals   = 35
    # data points      = 200
    # variables        = 7
    chi-square         = 151.115
    reduced chi-square = 0.783
    Akaike info crit   = -42.056
    Bayesian info crit = -18.967
[[Variables]]
    s1_amplitude:   1.07466323 +/- 0.089091 (8.29%) (init= 1.085042)
    s1_frequency:   0.01011808 +/- 0.000159 (1.57%) (init= 0.01015681)
    s1_phase:      -0.79396713 +/- 0.167738 (21.13%) (init=-0.8541416)
    s2_amplitude:   0.97797253 +/- 0.088986 (9.10%) (init= 0.9765722)
    s2_frequency:   0.02006821 +/- 0.000174 (0.86%) (init= 0.02008157)
    s2_phase:       1.00041090 +/- 0.191526 (19.14%) (init= 0.9813603)
    offset:         1.01786233 +/- 0.063693 (6.26%) (init= 1.009824)
[[Correlations]] (unreported correlations are <  0.100)
    C(s2_frequency, s2_phase)    = -0.876 
    C(s1_frequency, s1_phase)    = -0.864 
    C(s1_phase, s2_phase)        =  0.183 
    C(s1_frequency, offset)      =  0.133 
    C(s2_frequency, offset)      = -0.132 
    C(s1_amplitude, s2_frequency)  =  0.104 
    C(s2_phase, offset)          =  0.103 
    C(s1_phase, s2_frequency)    = -0.102 
    C(s1_amplitude, s2_phase)    = -0.102 
    C(s1_phase, offset)          = -0.100 







    

In [10]:

    

def three_sine_offset_testing():
    """ Testing procedure for the estimator for a three sine with offset fit. """

    x_axis = np.linspace(5, 300 ,200)

    phase1 = np.random.uniform()*2*np.pi
    ampl1 = 3
    freq1 = 0.03

    phase2 = np.random.uniform()*2*np.pi
    ampl2 = 2
    freq2 = 0.01

    phase3 = np.random.uniform()*2*np.pi
    ampl3 = 1
    freq3 = 0.05

    offset = 1.1

    data = (  ampl1 * np.sin(2*np.pi*freq1*x_axis + phase1)
            + ampl2 * np.sin(2*np.pi*freq2*x_axis + phase2)
            + ampl3 * np.sin(2*np.pi*freq3*x_axis + phase3)
            + offset)

    noisy_data = data + data.mean() * np.random.normal(size=x_axis.shape)*3

    x_dft1, y_dft1 = compute_dft(x_val=x_axis, y_val=noisy_data, zeropad_num=1)

    plt.figure()
    plt.plot(x_axis, noisy_data, 'o--', label='noisy_data')
    plt.plot(x_axis, data,'-', label='ideal data')
    plt.xlabel('Time (micro-s)')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    res1 = fitlogic.make_sine_fit(x_axis=x_axis, data=noisy_data, estimator=fitlogic.estimate_sine)
    data_sub1 = noisy_data - res1.best_fit

    x_dft2, y_dft2 = compute_dft(x_val=x_axis, y_val=data_sub1, zeropad_num=1)

    res2 = fitlogic.make_sine_fit(x_axis=x_axis, data=data_sub1, estimator=fitlogic.estimate_sine)
    data_sub2 = data_sub1 - res2.best_fit

    res3 = fitlogic.make_sine_fit(x_axis=x_axis, data=data_sub2, estimator=fitlogic.estimate_sine)

    x_dft3, y_dft3 = compute_dft(x_val=x_axis, y_val=data_sub2, zeropad_num=1)

    plt.figure()
    plt.plot(x_dft1, y_dft1, '-', label='noisy_data (3 peaks)')
    plt.plot(x_dft2, y_dft2, '-', label='noisy_data (2 peaks)')
    plt.plot(x_dft3, y_dft3, '-', label='noisy_data (1 peak)')
    plt.xlabel('Frequency (MHz)')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()


    mod, params = fitlogic.make_sinetriple_model()

    params['s1_amplitude'].set(value=res1.params['amplitude'].value)
    params['s1_frequency'].set(value=res1.params['frequency'].value)
    params['s1_phase'].set(value=res1.params['phase'].value)

    params['s2_amplitude'].set(value=res2.params['amplitude'].value)
    params['s2_frequency'].set(value=res2.params['frequency'].value)
    params['s2_phase'].set(value=res2.params['phase'].value)

    params['s3_amplitude'].set(value=res3.params['amplitude'].value)
    params['s3_frequency'].set(value=res3.params['frequency'].value)
    params['s3_phase'].set(value=res3.params['phase'].value)

    params['offset'].set(value=data.mean())

    result = mod.fit(noisy_data, x=x_axis, params=params)

    plt.figure()
    plt.plot(x_axis, result.best_fit,'-', label='fit')
    plt.plot(x_axis, noisy_data, 'o--', label='noisy_data')
    plt.plot(x_axis, data,'-', label='ideal data')
    plt.xlabel('Time micro-s')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    print(result.fit_report())


    

In [11]:

    

three_sine_offset_testing()


    

    




[[Model]]
    ((((Model(amplitude_function, prefix='s1_') * Model(bare_sine_function, prefix='s1_')) + (Model(amplitude_function, prefix='s2_') * Model(bare_sine_function, prefix='s2_'))) + (Model(amplitude_function, prefix='s3_') * Model(bare_sine_function, prefix='s3_'))) + Model(constant_function))
[[Fit Statistics]]
    # function evals   = 36
    # data points      = 200
    # variables        = 10
    chi-square         = 2071.580
    reduced chi-square = 10.903
    Akaike info crit   = 487.550
    Bayesian info crit = 520.533
[[Variables]]
    s1_amplitude:   3.27318297 +/- 0.333511 (10.19%) (init= 3.276196)
    s1_frequency:   0.02977691 +/- 0.000188 (0.63%) (init= 0.02984499)
    s1_phase:       0.73485320 +/- 0.208776 (28.41%) (init= 0.7022821)
    s2_amplitude:   2.46123748 +/- 0.336709 (13.68%) (init= 2.458016)
    s2_frequency:   0.01048051 +/- 0.000239 (2.29%) (init= 0.01047821)
    s2_phase:      -0.63118707 +/- 0.265912 (42.13%) (init=-0.6117081)
    s3_amplitude:   1.38892165 +/- 0.333296 (24.00%) (init= 1.376176)
    s3_frequency:   0.04948674 +/- 0.000438 (0.89%) (init= 0.04949179)
    s3_phase:      -1.95673665 +/- 0.482066 (24.64%) (init=-1.958989)
    offset:         1.05161290 +/- 0.234413 (22.29%) (init= 1.084918)
[[Correlations]] (unreported correlations are <  0.100)
    C(s1_frequency, s1_phase)    = -0.875 
    C(s3_frequency, s3_phase)    = -0.871 
    C(s2_frequency, s2_phase)    = -0.866 
    C(s1_phase, s2_phase)        =  0.134 
    C(s1_phase, s2_frequency)    = -0.103 
    C(s1_frequency, s2_phase)    = -0.103 







    













    













    

In [12]:

    

def three_sine_offset_testing2():
    """ Testing procedure for the implemented three sine with offset fit. """

    x_axis = np.linspace(5, 300 ,200)

    phase1 = np.random.uniform()*2*np.pi
    ampl1 = 3
    freq1 = 0.03

    phase2 = np.random.uniform()*2*np.pi
    ampl2 = 2
    freq2 = 0.01

    phase3 = np.random.uniform()*2*np.pi
    ampl3 = 1
    freq3 = 0.05

    offset = 1.1

    data = (ampl1 * np.sin(2*np.pi*freq1*x_axis + phase1) 
            + ampl2 * np.sin(2*np.pi*freq2*x_axis + phase2)
            + ampl3 * np.sin(2*np.pi*freq3*x_axis + phase3)
            + offset)

    noisy_data = data + data.mean() * np.random.normal(size=x_axis.shape)*3

    result = fitlogic.make_sinetriple_fit(
            x_axis=x_axis,
            data=noisy_data,
            estimator=fitlogic.estimate_sinetriple)

    plt.figure()
    plt.plot(x_axis, result.best_fit,'-', label='fit')
    plt.plot(x_axis, noisy_data, 'o--', label='noisy_data')
    plt.plot(x_axis, data,'-', label='ideal data')
    plt.xlabel('Time micro-s')
    plt.ylabel('signal')
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3, ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

    print(result.fit_report())


    

In [13]:

    

three_sine_offset_testing2()


    

    




[[Model]]
    ((((Model(amplitude_function, prefix='s1_') * Model(bare_sine_function, prefix='s1_')) + (Model(amplitude_function, prefix='s2_') * Model(bare_sine_function, prefix='s2_'))) + (Model(amplitude_function, prefix='s3_') * Model(bare_sine_function, prefix='s3_'))) + Model(constant_function))
[[Fit Statistics]]
    # function evals   = 58
    # data points      = 200
    # variables        = 10
    chi-square         = 1839.963
    reduced chi-square = 9.684
    Akaike info crit   = 463.837
    Bayesian info crit = 496.820
[[Variables]]
    s1_amplitude:   2.67779394 +/- 0.311896 (11.65%) (init= 2.645571)
    s1_frequency:   0.02993159 +/- 0.000221 (0.74%) (init= 0.02986787)
    s1_phase:       2.84106400 +/- 0.242473 (8.53%) (init= 2.904878)
    s2_amplitude:   1.52303945 +/- 0.306025 (20.09%) (init= 1.517533)
    s2_frequency:   0.00962267 +/- 0.000409 (4.25%) (init= 0.009667001)
    s2_phase:       0.02184890 +/- 0.439243 (2010.37%) (init=-0.02278683)
    s3_amplitude:   1.29938426 +/- 0.311706 (23.99%) (init= 1.298443)
    s3_frequency:   0.22373179 +/- 0.000444 (0.20%) (init= 0.2237347)
    s3_phase:       2.84558743 +/- 0.486880 (17.11%) (init= 2.843349)
    offset:         1.50326286 +/- 0.220968 (14.70%) (init= 1.518072)
[[Correlations]] (unreported correlations are <  0.100)
    C(s2_frequency, s2_phase)    = -0.879 
    C(s1_frequency, s1_phase)    = -0.876 
    C(s3_frequency, s3_phase)    = -0.871 
    C(s1_amplitude, s2_frequency)  = -0.126 
    C(s1_amplitude, s2_phase)    =  0.112 







    

In [ ]: