Savitzky-Golay滤波器（2）

前几天写过一篇介绍 Savitzky-Golay滤波器的文章, 没想到最近做项目还真的用上了。
因此就顺便写了个 C 语言的自动计算生成 SG 滤波器系数的程序。利用这里的代码可以生成任意阶数的 SG 滤波器。对于一些需要动态调整 SG 滤波器系数的场合，这里的代码还是很实用的。 上次的文章中给出了计算 SG 滤波器系数的 matlab 代码。

    x=[zeros(1,M),1,zeros(1,M)]  
    a=polyfit([-M:M],x,N)  
    polyval(a,[-M:M])  

只有三行，非常的简单。但是现在项目需要用 C 语言来实现。要是 用 C 语言从头写起这个代码会非常的长，光一个 ployfit 函数实现起来就很麻烦。我的程序中直接调用了 GSL 的相关函数。使得整个程序简化了很多。
先说说 polyfit 函数。多项式函数拟合可以转换成超定线性代数方程组的最小二乘解的问题。标准解法是 SVD 分解。不过GSL库中提供了 gsl_multifit_linear 函数进行线性模型的拟合，可以直接使用。下面给出代码。很简单，就不多做解释了。这个例子也可以作为 gsl_multifit_linear 函数用法的一个小例子。

/**
 * 多项式拟合函数
 * 根据 x y 拟合出一个N次多项式。返回多项式的系数。
 */
gsl_vector* PolyFit(gsl_vector *x, gsl_vector *y, int N)
{
    int i, j;
    int M = GSL_MIN(x->size, y->size);
    double chisq, xi;

    gsl_matrix *XX = gsl_matrix_alloc(M, N + 1);
    gsl_vector *c = gsl_vector_alloc(N + 1);
    gsl_matrix *cov = gsl_matrix_alloc(N + 1, N + 1);

    for(i = 0; i < M; i++)
    {
        gsl_matrix_set(XX, i, 0, 1.0);
        xi = gsl_vector_get(x, i);
        for(j = 1; j <= N; j++)
        {
            gsl_matrix_set(XX, i, j, pow(xi, j));
        }
    }
    gsl_multifit_linear_workspace *workspace = gsl_multifit_linear_alloc(M, N + 1);
    gsl_multifit_linear(XX, y, c, cov, &chisq, workspace);

    gsl_matrix_free(XX);
    gsl_matrix_free(cov);
    gsl_multifit_linear_free(workspace);
    return c;
}

有了 PolyFit 函数就可以计算 SG 滤波器的系数了。
Matlab 中的 polyval 函数 可以用 gsl_poly_eval 来实现。

下面是具体的代码。

/**
 * 计算 Savitzky-Golay 滤波器系数
 * SG 滤波器的阶数为 2M+1，多项式的最高次数为 N
 */
gsl_vector* SG_FilterCreate(int M, int N /* Poly Order */)
{
    int i;
    gsl_vector *x = gsl_vector_alloc(2 * M + 1);
    gsl_vector *y = gsl_vector_alloc(2 * M + 1);

    gsl_vector_set_zero(y);
    gsl_vector_set(y, M, 1);

    for(i= -M; i <= M; i++)
    {
        gsl_vector_set(x, i + M, i);
    }
    gsl_vector *c = PolyFit(x, y, N);
    gsl_vector_free(x);
    gsl_vector_free(y);

    gsl_vector *fir = gsl_vector_alloc(2 * M + 1);
    for(i = -M; i <= M; i++)
    {
        gsl_vector_set(fir, i + M, gsl_poly_eval (c->data, N + 1, i));
    }
    gsl_vector_free(c);
    return fir;
}