正弦波最小二乘曲线拟合是否可行（使用GSL）？

Question

是否可能使用GSL或类似的库来拟合A*sin(B*t+C)函数？

我想获取4096个样本（8位）中正弦波的A和C参数，并且可以提供B的良好近似值。

我认为可以使用GSL的非线性多次拟合来实现，但是我不理解所有与雅可比矩阵相关的数学背景...

- Benjamin

---

对于这种特殊情况，您也可以使用FFT，并查找振幅最大绝对值A的频率B。然后从复幅度的相位中获得相位C。 - Lutz Lehmann

Answer 1

是的，

您可能已经阅读过这篇文章：http://www.gnu.org/software/gsl/manual/html_node/Overview-of-Nonlinear-Least_002dSquares-Fitting.html#Overview-of-Nonlinear-Least_002dSquares-Fitting 您需要提供两个函数：

目标函数：

`

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获取的。 因此，Jacobian矩阵只是一个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)， 这是关于B的导数。

- alexandre iolov

---

sigma[]-向量的目的是什么？我在示例中看到它被初始化为0.1？ - Benjamin

我认为这是每个数据点的权重 f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i我不知道为什么他们将其设置为0.1，因为如果所有权重都相同（如示例中），则它们不会影响最小化。只有当某些数据点比其他数据点更准确时，它才有意义，此时它们将具有较低的sigma_i。 - alexandre iolov

Answer 2

谢谢Alexandre，你帮了我很多！

这是我现在正在使用的代码：

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

- Benjamin