可以进行正弦波最小二乘曲线拟合（使用 GSL）？答案

【问题标题】：Sine wave least squares curve fitting possible (using GSL)?可以进行正弦波最小二乘曲线拟合（使用 GSL）？
【发布时间】：2014-03-12 05:26:45
【问题描述】：

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

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

认为 GSL 非线性多重拟合应该是可能的，但我不了解所有雅可比矩阵的数学背景...

【问题讨论】：

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

标签： math curve-fitting gsl least-squares

【解决方案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

所以雅可比矩阵只是一个 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) 这是衍生品。乙

【讨论】：

sigma[]-Vector 的用途是什么？我在示例中看到它是用 0.1 初始化的？
我认为这是每个数据点的权重 f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i 我不知道他们为什么将其设置为.1，因为如果所有权重都相同（如示例中所示），则它们不会影响最小化。仅当某些数据点的准确度高于其他数据点时才有意义，在这种情况下，它们的 sigma_i 会更低

【解决方案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;
}

【讨论】：