Нулевая прокладка для фильтра Savitzgy Golay, не работающая для числового рецепта C++

Question

Моя проблема прямолинейна: Я хочу сгладить некоторые данные, используя фильтр Savitzgy Golay. Я использую C++. Код взят из книги 1 и может быть разделен на две части:

1. вычисления коэффициентов Savitzgy Голея и хранить их в векторном С.
2. Гладких S данных сигнала свертки его с C.

Проблема заключается в границах. Поскольку сигнал S не является периодическим, необходимо учитывать граничные эффекты. Это делается с так называемым 0-дополнением, что означает, что в конец сигнала добавляются дополнительные 0. Процедура описана точно в главе 13.1.1 от 1.

Однако я не могу найти полный пример этой процедуры, и моя собственная реализация, похоже, не работает, хотя я абсолютно не понимаю, почему. Ниже приведен хороший комментарий. Может ли кто-нибудь определить, что происходит на границах?

1 William H., et al. «Численные рецепты: научное искусство вычислений». (1987)

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

Вот визуализация выходного сигнала. Мы можем ясно видеть, что сбой заканчивается на границах, хотя учитывается нулевое дополнение.

Answer 1

Проблема заключается в том, что границы. Поскольку сигнал S не является периодическим, необходимо учитывать граничные эффекты. Это делается с так называемым 0-дополнением, что означает, что в конец сигнала добавляются дополнительные 0. Процедура описана точно в главе 13.1.1 из 1.

В моем издании Numerical Recipies глава 13 представляет собой «Фурье и спектральные приложения». В то время как нулевое заполнение сигнала идеально подходит для преобразования Фурье, это не очень хорошо для Савицки-Голая.

Я вижу несколько способов применить Савицкий-Гол сглаживание на границах сигнала:

1. Исключить недостающие биты сигнала из свертка. Установите коэффициенты, соответствующие отсутствующим битам, на ноль и измените нормализацию остальной части их на сумму до 0.
2. Вычислить специальное ядро ​​Savitzky-Golay для каждой точки сигнала с неполной окрестностью. На самом деле это не сложно. Концептуально свертка с ядром Савицки-Голея эквивалентна подгонке полинома к окрестности точки сигнала, а затем взятию этой точки сигнала из многочлена. Ничто не мешает вам иметь одностороннюю или асимметричную окрестность. Построение ядра Савицкого-Голея для произвольной окрестности - это вопрос об подгонке полинома к сигналу, где значение в начале координат равно 1 и нулю всюду. Происхождение не должно быть в центре окрестности. Коэффициенты ядра Савицкого-Голея являются тогда значениями установленной полиномиальной функции в соответствующих точках сигнала.

Answer 2

Я решил проблему сейчас, аналогичную предложению @ Tulon's 2.. В частности, я позабочусь о левой и правой границах, установив там дополнительный многочлен с обеих сторон. Это мотивируется внедрением фильтра Savitzgy Golay в библиотеке Scipy.signal для Python.