Savitzgy GolayフィルタのゼロパディングがC++の数値処方では機能しない

Question

問題はまっすぐです： Savitzgy Golayフィルタを使用してデータを滑らかにしたいですか？私はC + +を使用します。 コードブック1から取られ、2つの部分に分割することができる。

1. Cとを畳み込むことにより、信号データSスムースSavitzgyゴレイ係数を計算し、ベクトルC
2. に保存します

問題は境界です。信号Sは周期的ではないので、境界効果を考慮する必要がある。これは、いわゆる0パッディングで行われます。つまり、末尾の信号に余分な0が追加されます。この手順は、1の第13.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が追加されます。手順は1の第13.1.1章に正確に記載されています。

第13章数値レシピでは、第13章は「フーリエおよびスペクトルアプリケーション」です。フーリエ変換のために信号をゼロパディングするのはまったく問題はありませんが、Savitzky-Golayにとっては良い考えではありません。

は、私はいくつかの方法は、信号境界でサビツキー-Golay平滑化を適用するために参照してください。

1. は畳み込みからの信号の欠落ビットを除外します。欠けているビットに対応する係数をゼロに設定し、それらの残りの部分を1に合計するために再正規化します。
2. 不完全な近傍の各信号点に対して特別なSavitzky-Golayカーネルを計算します。それは実際には難しくありません。概念的には、Savitzky-Golayカーネルと畳み込むことは、多項式を信号点の近傍に当てはめ、その信号点を多項式からとることと等価です。片面または非対称の近所を避けることはできません。任意の近傍に対してSavitzky-Golayカーネルを構築するには、多項式を原点の値が1でどこでもゼロである信号に当てはめるという問題があります。起源は近所の中心にある必要はありません。 Savitzky-Golayカーネル係数は、対応する信号点における近似された多項式関数の値です。

Answer 2

@ Tulonの2つの提案に似た問題を解決しました。具体的には、両辺に余分な多項式を当てはめることによって、左右の境界を処理します。これは、Python用のScipy.signalライブラリのSavitzgy Golayフィルタの実装によって動機づけられます。