Geek asks and answers

Is a least squares finish approximation possible (using GSL)?

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

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

I think that it should be possible with non-linear GSL multi-fit, but I do not understand the mathematical background with all the Jacobian matrix material ...

+3
source share
2 answers

Yes,

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

You are required to provide two functions

target:

`

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]);
    }
    ...
    }


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

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

, - 3xN, N - , J (0,3) = sin (B * t_3 + C)

A, B, C x [0], x [1], x [2]

J (1,5) = A * t_5 * cos (B * t_5 + C) .

+2

, , !

, :

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;
}
0

More articles:


All Articles