Sine Wave Curve Fit Question

[Iain Mackay]

Python Folks

I'm a newbie to Python and am looking for a library / function that can help
me fit a 1D data vector to a sine wave. I know the frequency of the wave,
so its really only phase and amplitude information I need.

I can't find anything in the most widely known libraries (they seem to be
strong on polynomial fitting, but not, apparently, on trig functions) and I
wondered if any one here had recommendations?

Something that implemented IEEE 1057 , or similar, would be perfect.


TIA
Iain

 

[marek.rocki]

Iain Mackay napisal(a):

Python Folks

I'm a newbie to Python and am looking for a library / function that can help
me fit a 1D data vector to a sine wave. I know the frequency of the wave,
so its really only phase and amplitude information I need.

I can't find anything in the most widely known libraries (they seem to be
strong on polynomial fitting, but not, apparently, on trig functions) and I
wondered if any one here had recommendations?

Something that implemented IEEE 1057 , or similar, would be perfect.


TIA
Iain

I'm not aware of any specialized library, but how about using numpy
and running FFT on your data? If you already know the frequency, you
can easily extract phase and scaled amplitude from the result.

Regards,
Marek

 

[Helmut Jarausch]

Python Folks

I'm a newbie to Python and am looking for a library / function that can help
me fit a 1D data vector to a sine wave. I know the frequency of the wave,
so its really only phase and amplitude information I need.

I can't find anything in the most widely known libraries (they seem to be
strong on polynomial fitting, but not, apparently, on trig functions) and I
wondered if any one here had recommendations?

Something that implemented IEEE 1057 , or similar, would be perfect.

Let's do a bit math first.

Your model is A*sin(omega*t+alpha) where A and alpha are sought.
Let T=(t_1,...,t_N)' and Y=(y_1,..,y_N)' your measurements (t_i,y_i)
( ' denotes transposition )

First, A*sin(omega*t+alpha) =
A*cos(alpha)*sin(omega*t) + A*sin(alpha)*cos(omega*t) =
B*sin(omega*t) + D*cos(omega*t)

by setting B=A*cos(alpha) and D=A*sin(alpha)

Once, you have B and D, tan(alpha)= D/B A=sqrt(B^2+D^2)
Then in vector notation S=sin(omega*T) C=cos(omega*T)
you get the 2x2 system for B and D :

(S'*S) * B + (S'*C) * D = S'*Y
(S'*C) * B + (C'*C) * D = C'*Y

where S'*C is the scalar product of the vectors S and C and similarly.

Now, for Python, to handle vectors and scalar products efficiently, have a look
at numpy.



--
Helmut Jarausch

Lehrstuhl fuer Numerische Mathematik
RWTH - Aachen University
D 52056 Aachen, Germany

 

[Iain Mackay]

Thanks folks - I'll have a think about both of these options.

 

[Mikael Olofsson]

Your model is A*sin(omega*t+alpha) where A and alpha are sought.
Let T=(t_1,...,t_N)' and Y=(y_1,..,y_N)' your measurements (t_i,y_i)
( ' denotes transposition )

First, A*sin(omega*t+alpha) =
A*cos(alpha)*sin(omega*t) + A*sin(alpha)*cos(omega*t) =
B*sin(omega*t) + D*cos(omega*t)

by setting B=A*cos(alpha) and D=A*sin(alpha)

Once, you have B and D, tan(alpha)= D/B A=sqrt(B^2+D^2)

This is all very true, but the equation tan(alpha)=D/B may fool you.
This may lead you to believe that alpha=arctan(D/B) is a solution, which
is not always the case. The point (B,D) may be in any of the four
quadrants of the plane. Assuming B!=0, the solutions to this equation
fall into the two classes

alpha = arctan(D/B) + 2*k*pi

and

alpha = arctan(D/B) + (2*k+1)*pi,

where k is an integer. The sign of B tells you which class gives you the
solution. If B is positive, the solutions are those in the first class.
If B is negative, the solutions are instead those in the second class.
Whithin the correct class, you may of course choose any alternative.

Then we have the case B=0. Then the sign of D determines alpha. If D is
positive, we have alpha=pi/2, and if D is negative, we have alpha=-pi/2.

Last if both B and D are zero, any alpha will do.

/MiO

 

[Helmut Jarausch]

This is all very true, but the equation tan(alpha)=D/B may fool you.

You're right: standard Python's math library missing the function arctan2.
But you can get it from numpy or numarray.
In case of doubt just compute the residuum with different possibilities .

This may lead you to believe that alpha=arctan(D/B) is a solution, which
is not always the case. The point (B,D) may be in any of the four
quadrants of the plane. Assuming B!=0, the solutions to this equation
fall into the two classes

alpha = arctan(D/B) + 2*k*pi

and

alpha = arctan(D/B) + (2*k+1)*pi,

where k is an integer. The sign of B tells you which class gives you the
solution. If B is positive, the solutions are those in the first class.
If B is negative, the solutions are instead those in the second class.
Whithin the correct class, you may of course choose any alternative.

Then we have the case B=0. Then the sign of D determines alpha. If D is
positive, we have alpha=pi/2, and if D is negative, we have alpha=-pi/2.

Last if both B and D are zero, any alpha will do.

/MiO

--
Helmut Jarausch

Lehrstuhl fuer Numerische Mathematik
RWTH - Aachen University
D 52056 Aachen, Germany

 

[Paul Rubin]

You're right: standard Python's math library missing the function arctan2.

It's math.atan2 .