Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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?

share|improve this question

migrated from astronomy.stackexchange.com Aug 18 at 17:24

This question came from our site for astronomers and astrophysicists.

    
You said by yourself: the fit. –  Py-ser Aug 12 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 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 at 6:00
    
Maybe this is a nice link to folow: math.unl.edu/~sdunbar1/sinefit.html –  agtoever Aug 12 at 10:03
    
This really should go to stats.SE, not math.SE. Don't know why it was forwarded here. –  heropup Aug 18 at 17:36

3 Answers 3

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

share|improve this answer

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.

share|improve this answer

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
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.