Fitting a variable Sinc function in python

Question

I would like to fit a sinc function to a bunch of datalines. Using a gauss the fit itself does work but the data does not seem to be sufficiently gaussian, so I figured I could just switch to sinc..

I just tried to put together a short piece of self running code but realized, that I probably do not fully understand, how arrays are handled if handed over to a function, which could be part of the reason, why I get error messages calling my program

So my code currently looks as follows:

from numpy import exp
from scipy.optimize import curve_fit
from math import sin, pi

def gauss(x,*p):
    print(p)
    A, mu, sigma = p
    return A*exp(-1*(x[:]-mu)*(x[:]-mu)/sigma/sigma)

def sincSquare_mod(x,*p):
    A, mu, sigma = p
    return A * (sin(pi*(x[:]-mu)*sigma) / (pi*(x[:]-mu)*sigma))**2


p0 = [1., 30., 5.]
xpos = range(100)
fitdata = gauss(xpos,p0)
p1, var_matrix = curve_fit(sincSquare_mod, xpos, fitdata, p0)

What I get is:

Traceback (most recent call last):
File "orthogonal_fit_test.py", line 18, in <module>
fitdata = gauss(xpos,p0)
File "orthogonal_fit_test.py", line 7, in gauss
A, mu, sigma = p
ValueError: need more than 1 value to unpack

From my understanding p is not handed over correctly, which is odd, because it is in my actual code. I then get a similar message from the sincSquare function, when fitted, which could probably be the same type of error. I am fairly new to the star operator, so there might be a glitch hidden...

Anybody some ideas? :)

Thanks!

Answer 1

You need to make three changes,

def gauss(x, A, mu, sigma):
    return A*exp(-1*(x[:]-mu)*(x[:]-mu)/sigma/sigma)

def sincSquare_mod(x, A, mu, sigma):
    x=np.array(x)
    return A * (np.sin(pi*(x[:]-mu)*sigma) / (pi*(x[:]-mu)*sigma))**2

fitdata = gauss(xpos,*p0)

1, See Documentation

2, replace sin by the numpy version for array broadcasting

3, straight forward right? :P

Note, i think you are looking for p1, var_matrix = curve_fit(gauss,... rather than the one in the OP, which appears do not have a solution.

Answer 2

Also worth noting is that you will get rounding errors as x*Pi gets close to zero that might get magnified. You can approximate as demonstrated below for better results (VB.NET, sorry):

Private Function sinc(x As Double) As Double
    x = (x * Math.PI)
    'The Taylor Series expansion of Sin(x)/x is used to limit rounding errors for small values of x
    If x < 0.01 And x > -0.01 Then
        Return 1.0 - x ^ 2 / 6.0 + x ^ 4 / 120.0
    End If
    Return Math.Sin(x) / x
End Function

http://www.wolframalpha.com/input/?i=taylor+series+sin+%28x%29+%2F+x&dataset=&equal=Submit