[lnkForumImage]
TotalShareware - Download Free Software
Confronta i prezzi di migliaia di prodotti.
Asp Forum
 Home | Login | Register | Search 


 

Forums >

comp.lang.python

Sine Wave Curve Fit Question

Ian Mackay

1/30/2008 2:26:00 PM

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
6 Answers

marek.rocki

1/30/2008 6:06:00 PM

0

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

1/30/2008 6:58:00 PM

0

Iain Mackay wrote:
> 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

Ian Mackay

1/31/2008 6:55:00 AM

0

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

Mikael Olofsson

1/31/2008 8:47:00 AM

0

Helmut Jarausch wrote:
> 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

1/31/2008 3:52:00 PM

0

Mikael Olofsson wrote:
> Helmut Jarausch wrote:
>> 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.

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

1/31/2008 3:59:00 PM

0

Helmut Jarausch <jarausch@skynet.be> writes:
> You're right: standard Python's math library missing the function arctan2.

It's math.atan2 .