Sine wave least squares curve fitting possible (using GSL)?

Question

Is it possible to fit an A*sin(B*t+C) function with GSL or a similar library?

i want to get the A and C parameter of a sine wave present in 4096 samples (8bit) and can provide an good approximation of B.

A think that should be possible with GSLs non-linear multifit but I don’t understand the mathematical background with all that Jacobian matrix stuff...

Answer 1

Yes,

You have probably read this: http://www.gnu.org/software/gsl/manual/html_node/Overview-of-Nonlinear-Least_002dSquares-Fitting.html#Overview-of-Nonlinear-Least_002dSquares-Fitting

What is required from you is to provide two functions

the objective:

`

int sine_f (const gsl_vector * x, void *data, 
        gsl_vector * f){
    ...
    for(...){
    ...
      double Yi = A * sin(B*t +C);
      gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
    }
    ...
    }

and then the derivative of the objective with respect to the parameters

int
sine_df (const gsl_vector * x, void *data, 
         gsl_matrix * J)
//the derivatives of Asin(Bt +C) wrt A,B,C for each t

This is straight from http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Nonlinear-Least_002dSquares-Fitting.html#Example-programs-for-Nonlinear-Least_002dSquares-Fitting

So the Jacobian is just a 3xN matrix, where N is the number of data points For example J(0,3) = sin(B*t_3 + C)

if A,B,C correspond to x[0],x[1],x[2]

And J(1,5) = A*t_5*cos(B*t_5 + C) This is the derivative wrt. B

Answer 2

Thank you alexandre, you helped me a lot!

Here is the code I'am using now:

typedef struct{
  uint32    u32_n;
  float64*  pf64_y;
  float64*  pf64_sigma;
}ST_DATA;

int expb_f (const gsl_vector* x, void* p_data, gsl_vector* f)
{
    ST_DATA* pst_data = (ST_DATA*)p_data;
    uint32   u32_n      = pst_data->u32_n;
    float64* pf64_y     = pst_data->pf64_y;
    float64* pf64_sigma = pst_data->pf64_sigma;

    float64  A  = /* x[0] */ gsl_vector_get (x, 0);
    float64  B  = /* x[1] */ gsl_vector_get (x, 1);
    float64  C  = /* x[2] */ gsl_vector_get (x, 2);
    float64  Yi,Fi;
    uint32 i;
    for (i=0; i<u32_n; i++)
    {
        Yi = A * sin(B*i + C);
        Fi = (Yi - pf64_y[i])/pf64_sigma[i];
        /* f[i] = Fi; */    gsl_vector_set(f,i,Fi);
    }

    return GSL_SUCCESS;
}

int expb_df (const gsl_vector* x, void* p_data, gsl_matrix* J)
{
    ST_DATA* pst_data = (ST_DATA*)p_data;
    uint32   u32_n      = pst_data->u32_n;
    float64* pf64_y     = pst_data->pf64_y;
    float64* pf64_sigma = pst_data->pf64_sigma;

    float64  A  = /* x[0] */ gsl_vector_get (x, 0);
    float64  B  = /* x[1] */ gsl_vector_get (x, 1);
    float64  C  = /* x[2] */ gsl_vector_get (x, 2);
    float64  Yi;
    uint32 i;
    for (i=0; i<u32_n; i++)
    {
        /* J[i][0] =    sin(B*i+C);   */gsl_matrix_set (J, i, 0,   sin(B*i+C)  );
        /* J[i][1] =  A*cos(B*i+C)*i; */gsl_matrix_set (J, i, 1, A*cos(B*i+C)*i);
        /* J[i][2] =  A*cos(B*i+C);   */gsl_matrix_set (J, i, 2, A*cos(B*i+C)  );
    }
  return GSL_SUCCESS;
}

int expb_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J)
{
  expb_f(x, data, f);
  expb_df(x, data, J);

  return GSL_SUCCESS;
}