c++ - Savitzky Golay滤波器的零填充不适用于C++数值配方

我的问题很简单:
我想使用Savitzgy Golay过滤器平滑一些数据。我使用C++。
该代码摘自1书，可分为两部分:

计算Savitzgy Golay系数并将其存储在 vector C中。

通过与C卷积来平滑信号数据S。

问题在于界限。由于信号S不是周期性的，因此必须考虑边界效应。这是通过所谓的0填充完成的，这意味着最后会将一些额外的0加到信号上。该过程在1的第13.1.1章中进行了详细描述。

但是，我找不到该过程的完整示例，尽管我绝对无法理解为什么，但我自己的实现似乎无法正常工作。下面是一个很好的例子。有人可以发现边界出了什么问题吗？

#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <cstdlib>
#include <string>
#include <fstream>
#include <vector>

#include "./numerical_recipes/other/nr.h"
#include "./numerical_recipes/recipes/savgol.cpp"
#include "./numerical_recipes/recipes/lubksb.cpp"
#include "./numerical_recipes/recipes/ludcmp.cpp"
#include "./numerical_recipes/recipes/convlv.cpp"
#include "./numerical_recipes/recipes/realft.cpp"
#include "./numerical_recipes/recipes/four1.cpp"

using namespace std;


int main()
{
    // set savgol parameters
    int nl = 6;  // left window length
    int nr = 6;  // right window length
    int m  = 3;  // order of interpolation polynomial

    // calculate savitzgy golay coefficients
    int np=nl+nr+1;         // number of coefficients
    Vec_DP coefs(np);       // vector that stores the coefficients
    NR::savgol(coefs,np,nl,nr,0,m); // calculate the coefficients

    // as example input data, generate sinh datapoints between -1 and 1
    int nvals = int(pow(2,7))-nl; // number of datapoints to analyze (equal to 2^7 including zero-padding)
    Vec_DP args(nvals); // stores arguments
    Vec_DP vals(nvals); // stores signal
    double next_arg; // help variable
    for(int i = 0; i < nvals; i++)
    {
        next_arg = i*2./(nvals-1)-1;    // next argument
        args[i] = next_arg;             // store argument point
        vals[i] = sinh(next_arg);        // evaluate next value
    }

    // for zero padding, we have to add nl datapoints to the right. The signal is then of length 2^7.
    // see also chapter 13.1.1 in [1]
    // [1] Press, William H., et al. "Numerical recipes: the art of scientific computing." (1987)
    Vec_DP input_signal(int(pow(2,7))); // create vector of length 2^7
    for(int i = 0; i < nvals; i++) input_signal[i] = vals[i]; // overwrite with actual signal
    for(int i = nvals; i < int(pow(2,7)); i++) input_signal[i] = 0.0; // add zeros for zero-patting

    // perfrom the convolution
    Vec_DP ans(int(pow(2,7)));  // stores the smoothed signal
    NR::convlv(input_signal,coefs,1,ans); // smoothen the data

    // write data to the output for visual inspection
    string filename = "test.csv"; // output filename
    string write_line;
    ofstream wto(filename,ios::app);
    for(int i = 0; i < nvals; i++) // write result to output, drop the values from 0-padding
    {
        write_line = to_string(args[i])+", "+to_string(vals[i])+= ", "+to_string(ans[i]);
        wto << write_line << endl;
    }
    wto.close();

    return 0;
}

这是输出的可视化。我们可以清楚地看到，尽管考虑了零填充，但拟合在边界处失败了。

最佳答案

在我的《数值接收》一书中，第13章是“傅立叶和频谱应用”。尽管零填充信号对于傅立叶变换来说是完美的，但对于Savitzky-Golay来说并不是一个好主意。

我看到了几种在信号边界应用Savitzky-Golay平滑的方法:

从卷积中排除信号的丢失位。将与丢失位相对应的系数设置为零，然后将剩余的系数重新归一化为1。

为具有不完整邻域的每个信号点计算一个特殊的Savitzky-Golay内核。实际上，这并不难。从概念上讲，使用Savitzky-Golay内核进行卷积等效于将多项式拟合到信号点的附近，然后从多项式中获取该信号点。没有什么可以阻止您拥有单侧或不对称邻域。为任意邻域构建Savitzky-Golay内核是将多项式拟合到信号的原点的问题，该信号的原点值为1，而其他各处为零。原点不必位于邻域的中心。那么，Savitzky-Golay核系数就是相应信号点处的拟合多项式函数的值。

关于c++ - Savitzky Golay滤波器的零填充不适用于C++数值配方，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/41980303/