fit sine wave to data

Question

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?

Answer 1

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$.

Answer 2

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.

Answer 3

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

Answer 4

Two samples with $90^\circ$ difference method for $y=A\sin(\omega t+\pi)$ and $T$ as the period

first sample : $y_1=A\sin(\omega t+\pi)$
second sample : $y_2=A\sin\left(\omega\left[t+\frac T4\right]+\pi\right)=A\cos(\omega t+\pi)$
so $A^2=y_1^2+y_2^2$

B.R s.mosaviasl