C- 通过二次拟合进行峰值检测

Question

我有一个应用程序，我需要在给定数据集中找到峰值的位置。分辨率必须远高于数据点之间的间距（即，仅找到最高数据点是不够的，而是必须根据峰的形状估计“虚拟”峰位置）。一个峰值由大约 4 或 5 个数据点组成。每隔几毫秒获取一个数据集，并且必须实时执行峰值检测。

我比较了 LabVIEW 中的几种方法，发现最好的结果（在分辨率和速度方面）由 LabVIEW PeakDetector.vi 给出，它使用移动窗口（>= 3 点宽度）扫描数据集，并且对于每个position 执行二次拟合。得到的二次函数（抛物线）有一个局部最大值，然后与附近的点进行比较。

现在我想在C中实现相同的方法。多项式拟合实现如下（使用高斯矩阵）：

// Fits *y from x_start to (x_start + window) with a parabola and returns x_max and y_max
int polymax(uint16_t * y_data, int x_start, int window, double *x_max, double *y_max)
{
    float sum[10],mat[3][4],temp=0,temp1=0,a1,a2,a3;
    int i,j;

    float x[window];
    for(i = 0; i < window; i++)
        x[i] = (float)i;

    float y[window];
    for(i = 0; i < window; i++)
        y[i] = (float)(y_data[x_start + i] - y_data[x_start]);

    for(i = 0; i < window; i++)
    {
        temp=temp+x[i];
        temp1=temp1+y[i];
    }
    sum[0]=temp;
    sum[1]=temp1;
    sum[2]=sum[3]=sum[4]=sum[5]=sum[6]=0;

    for(i = 0;i < window;i++)
    {
        sum[2]=sum[2]+(x[i]*x[i]);
        sum[3]=sum[3]+(x[i]*x[i]*x[i]);
        sum[4]=sum[4]+(x[i]*x[i]*x[i]*x[i]);
        sum[5]=sum[5]+(x[i]*y[i]);
        sum[6]=sum[6]+(x[i]*x[i]*y[i]);
    }
    mat[0][0]=window;
    mat[0][1]=mat[1][0]=sum[0];
    mat[0][2]=mat[1][2]=mat[2][0]=sum[2];
    mat[1][2]=mat[2][3]=sum[3];
    mat[2][2]=sum[4];
    mat[0][3]=sum[1];
    mat[1][3]=sum[5];
    mat[2][3]=sum[6];

    temp=mat[1][0]/mat[0][0];
    temp1=mat[2][0]/mat[0][0];
    for(i = 0, j = 0; j < 3 + 1; j++)
    {
        mat[i+1][j]=mat[i+1][j]-(mat[i][j]*temp);
        mat[i+2][j]=mat[i+2][j]-(mat[i][j]*temp1);
    }

    temp=mat[2][4]/mat[1][5];
    temp1=mat[0][6]/mat[1][7];
    for(i = 1,j = 0; j < 3 + 1; j++)
    {
        mat[i+1][j]=mat[i+1][j]-(mat[i][j]*temp);
        mat[i-1][j]=mat[i-1][j]-(mat[i][j]*temp1);
    }

    temp=mat[0][2]/mat[2][2];
    temp1=mat[1][2]/mat[2][2];
    for(i = 0, j = 0; j < 3 + 1; j++)
    {
        mat[i][j]=mat[i][j]-(mat[i+2][j]*temp);
        mat[i+1][j]=mat[i+1][j]-(mat[i+2][j]*temp1);
    }

    a3 = mat[2][3]/mat[2][2];
    a2 = mat[1][3]/mat[1][8];
    a1 = mat[0][3]/mat[0][0];

    // zX^2 + yX + x

    if (a3 < 0)
    {
        temp = - a2 / (2*a3);
        *x_max = temp + x_start;
        *y_max = (a3*temp*temp + a2*temp + a1) + y_data[x_start];
        return 0;
    }
    else
       return -1;
}

扫描在一个外部函数中执行，它重复调用上述函数并选择最高的局部y_max。

上面的作品和山峰都找到了。只有噪声比 LabVIEW 对应物差得多（即，给定相同的输入数据集和相同的参数，我得到一个非常振荡的峰值位置）。由于算法有效，上述代码在概念上应该是正确的，所以我认为这可能是一个数值问题，因为我只是使用“浮点数”而没有进一步努力提高数值准确性。这是一个可能的答案吗？有没有人有提示，我应该在哪里寻找？

谢谢。

PS：我已经完成了搜索，发现这个非常 good overview 和 this question，类似于我的（不幸的是没有很多答案）。我会进一步研究这些。

编辑：我发现我的问题在其他地方。通过删除某些输出值来改进算法（一种只有在结果在移动窗口内时才接受结果的后验证）带来了问题的解决方案。现在我对结果感到满意，即它们与 LabVIEW 的结果相当。不过，非常感谢您的 cmets。

Answer 1

很抱歉迟到了，但如果您有 C/C++，使用 VS2013 Express（免费版）将其移植到 C# 代码并使用 .NET 工具集将其移植到 Labview 非常容易。