fit sine wave to data

I have some astronomical data which I know has a sinusoidal component to it of the form

y = Asin(Ωt+Φ). The period of the sinusoid is equal to a sidereal day. So I know the frequency Ω. So, I just need a way to find out optimal values for A and Φ.

Now, since the period of the sinusoid is so huge and we can track a source only to an extent of about 6 hours, I never get data that would have a full cycle of a sine wave in it. most of the cases I might have only 25 percent of it, at the maximum. In this situation where I don't even have a full cycle, what is the best way to estimate the value of the amplitude and phase of the sinusoid?

-

migrated from astronomy.stackexchange.comAug 18 '14 at 17:24

This question came from our site for astronomers and astrophysicists.

You said by yourself: the fit. – Py-ser Aug 12 '14 at 5:50
If your data covers 1/4 of the period, you should be able to determine the phase fairly easily. Unless you collect data within 20° of the peak, it would be difficult to calculate the amplitude. – LDC3 Aug 12 '14 at 5:56
@LDC3 can you tell me the procedure how to go about it so that its failsafe at all sections of the sinusoid? – Vivek V K Aug 12 '14 at 6:00
Maybe this is a nice link to folow: math.unl.edu/~sdunbar1/sinefit.html – agtoever Aug 12 '14 at 10:03
This really should go to stats.SE, not math.SE. Don't know why it was forwarded here. – heropup Aug 18 '14 at 17:36

Write $y(t) = A\sin(\Omega t)\cos(\phi) + A\cos(\Omega t)\sin(\phi)$. Let $A_1 = A\cos(\phi)$, $A_2 = A\sin(\phi)$. Define the variables $w(t) = \sin(\Omega t)$ and $z(t) = \cos(\Omega t)$. Now you can write $$y(t) = \begin{bmatrix}w(t) & z(t)\end{bmatrix}\begin{bmatrix}A_1 \\ A_2\end{bmatrix}$$

The rest is linear regression. To be more clear you can write your datapoints as $$\begin{bmatrix}y_1 \\ \vdots \\ y_n\end{bmatrix} = \begin{bmatrix}w_1 & z_1\\ \vdots & \vdots\\ w_n & z_n \end{bmatrix}\begin{bmatrix}A_1 \\ A_2\end{bmatrix}$$ Using the matrix notation $$Y = XB$$ You want to solve the following minimization problem $$\min_{B} \lVert Y - XB \rVert = \min_{B} (Y-XB)^\top (Y-XB)$$

If you differentiate the cost function with respect to $B$ and set it to $0$, you will see that $$B = (X^\top X)^{-1}X^\top Y$$.

Having $B$ you can solve for $A$ and $\phi$.

-

If $y(t) = A\sin(\Omega t + \phi)$, then $y'(t) = \Omega A\cos(\Omega t + \phi)$ and hence

$$A^2 = y(t)^2 + \frac{1}{\Omega^2} y'(t)^2$$

You just need to have a rough estimate of $y(0)$ and $y'(0)$ (or same information at some other $t$), plug it into above formula to get $A$. $\phi$ is then equal to $\sin^{-1}( y(0) / A )$. If you want to improve accuracy, you can peform a least square fit of the from

$$y = \tilde{A}\sin(\tilde{\Omega} t + \tilde{\phi})$$ with $( \tilde{A}, \tilde{\Omega}, \tilde{\phi} )$ near the numbers $(A, \Omega, \phi)$ you just estimated.

-

What I have done before is expand the equation to $y(t) = S \sin(\Omega t) + C \cos(\Omega t)$ and then calculate the following integrals:

\begin{align} \int_0^{\frac{2\pi}{\Omega}} \sin(\Omega t) y(t)\,{\rm d}t & = \frac{\pi S}{\Omega} \\ \int_0^{\frac{2\pi}{\Omega}} \cos(\Omega t) y(t)\,{\rm d}t & = \frac{\pi C}{\Omega} \end{align}

So if you need average of $n$ cycles you have

\begin{align} S & = \frac{\Omega}{n \pi} \int_0^{\frac{2\pi\,n}{\Omega}} \sin(\Omega t) y(t)\,{\rm d}t \\ C & = \frac{\Omega}{n \pi} \int_0^{\frac{2\pi\,n}{\Omega}} \cos(\Omega t) y(t)\,{\rm d}t \end{align}

I have deployed a simple trapezoidal numerical integrator with very good results. Of course this depends on how many data points you have for each cycle. The more points then the more noise you have and the less you might get skewed results due to aliasing.

In the end your $A=\sqrt{S^2+C^2}$ and $\Phi = \arctan\left( \frac{C}{S} \right)$

For the fun of it here some VBA code I used in Excel to do exactly this:

Const PI As Double = 3.14159265358979

'---------------------------------------------------------------------------------------
' Procedure : IntegrateSin
' Author    : ja72
' Date      : 7/22/2014
' Purpose   : Perform a integral of SIN(x)*f(x) between x_1 and x_2
'             To be used as a formula like: "=1/180*IntegrateSin(X$55:X$83,0,360)"
'---------------------------------------------------------------------------------------
Public Function IntegrateSin(ByRef r_f As Range, x_1 As Double, x_2 As Double) As Double
Dim i As Integer, N As Integer, sum As Double, h As Double, f() As Variant, x As Double
N = r_f.Rows.Count - 1
h = (x_2 - x_1) / N
' Convert range of cells into array f
f = r_f.Value: sum = 0#
' Transform f values to sin(x)*f
For i = 1 To N + 1
x = x_1 + h * (i - 1)
f(i, 1) = Sin(PI / 180# * x) * f(i, 1)
Next i
' Trapezoidal integration
sum = sum + h * f(1, 1) / 2
For i = 2 To N
sum = sum + h * f(i, 1)
Next i
sum = sum + h * f(N + 1, 1) / 2
IntegrateSin = sum
End Function

'---------------------------------------------------------------------------------------
' Procedure : IntegrateCos
' Author    : ja72
' Date      : 7/22/2014
' Purpose   : Perform a integral of COS(x)*f(x) between x_1 and x_2
'             To be used as a formula like: "=1/180*IntegrateCos(X$55:X$83,0,360)"    '---------------------------------------------------------------------------------------
Public Function IntegrateCos(ByRef r_f As Range, x_1 As Double, x_2 As Double) As Double
Dim i As Integer, N As Integer, sum As Double, h As Double, f() As Variant, x As Double
N = r_f.Rows.Count - 1
h = (x_2 - x_1) / N
' Convert range of cells into array f
f = r_f.Value: sum = 0#
' Transform f values to cos(x)*f
For i = 1 To N + 1
x = x_1 + h * (i - 1)
f(i, 1) = Cos(PI / 180# * x) * f(i, 1)
Next i
' Trapezoidal integration
sum = sum + h * f(1, 1) / 2
For i = 2 To N
sum = sum + h * f(i, 1)
Next i
sum = sum + h * f(N + 1, 1) / 2
IntegrateCos = sum
End Function
-