サビツキー・ゴレイフィルタ

Savitzky -Golay フィルタは、デジタル データポイントのセットに適用してデータを平滑化できるデジタル フィルタです。つまり、信号の傾向を歪めることなくデータの精度を高めることができます。これは、畳み込みと呼ばれるプロセスで、線形最小二乗法によって、隣接するデータ ポイントの連続サブセットを低次多項式に当てはめることによって実現されます。データ ポイントが等間隔に配置されている場合は、最小二乗方程式の解析解が、すべてのデータ サブセットに適用できる単一の「畳み込み係数」セットの形で見つかります。これにより、各サブセットの中心点で平滑化された信号 (または平滑化された信号の導関数) の推定値が得られます。この方法は、確立された数学的手順に基づいており、^{[1] [2]} 1964年にさまざまな多項式とサブセットのサイズに対する畳み込み係数の表を発表したアブラハム・サヴィツキーとマルセル・J・E・ゴレイによって普及しました。 ^{[3] [4]}表のいくつかのエラーは修正されています。^[5]この方法は、2次元および3次元データの処理に拡張されています。

サヴィツキーとゴレイの論文は、分析化学^誌[6]で最も引用されている論文の一つであり、同誌では「コンピューター制御の分析機器の黎明期はこの論文に遡ると言える」と評され、「10の重要な論文」の一つに分類されている^{[7] 。}

データは点の集合、で構成され、 は独立変数、 は観測値である。それらは、次の式に従って 畳み込み係数の集合、で処理される。${\displaystyle \{x_{j}}$${\displaystyle y_{j}\};j=1,...,n}$${\displaystyle x_{j}}$${\displaystyle y_{j}}$${\displaystyle m}$${\displaystyle C_{i}}$

${\displaystyle Y_{j}=\sum _{i={\tfrac {1-m}{2}}}^{\tfrac {m-1}{2}}C_{i}\,y_{j+i},\qquad {\frac {m+1}{2}}\leq j\leq n-{\frac {m-1}{2}}}$

選択された畳み込み係数は、以下の表に示されています。たとえば、5点の2次多項式による平滑化の場合、平滑化されたデータポイントは次のように表されます。 ${\displaystyle m=5,i=-2,-1,0,1,2}$${\displaystyle j^{th}}$${\displaystyle Y_{j}}$

${\displaystyle Y_{j}={\frac {1}{35}}(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+ 2})}$、

ここで、、などです。平滑化には、アルゴリズム チェーンを通じてノイズが伝播するのを防ぐなど、さまざまな用途があります。また、単にデータのノイズを実際よりも少なく見せるだけの用途もあります。 ${\displaystyle C_{-2}=-3/35,C_{-1}=12/35}$

以下はデータの数値微分法の応用例である。^[8] 注n次の導関数を計算する場合、絶対値を得るために、計算されたすべてのデータポイントに の追加のスケーリング係数を適用することができる（詳細については、以下の の式を参照）。 ${\displaystyle {\frac {n!}{h^{n}}}}$${\displaystyle {\frac {d^{n}Y}{dx^{n}}}}$

1. 実験データ曲線における最大値と最小値の位置。これがサビツキーを最初に刺激した応用でした。^[4]関数の一次導関数は、最大値または最小値でゼロになります。この図は、合成ローレンツ曲線に属するデータ ポイントを示しており、ノイズ (青いダイヤモンド) が追加されています。データは、ゼロでのピーク最大値を基準として、半分の幅のスケールでプロットされています。平滑化された曲線 (赤い線) と一次導関数 (緑) は、7 ポイントの 3 次サビツキー ゴレイ フィルターを使用して計算されました。ゼロ交差の両側の位置での一次導関数値の線形補間により、ピーク最大値の位置が得られます。この目的で 3 次導関数を使用することもできます。
2. 滴定曲線における終点の位置。終点は関数の2次導関数がゼロになる変曲点です。 ^[9]マロン酸の滴定曲線はこの方法の威力を示しています。4mlの最初の終点はほとんど見えませんが、2次導関数により、ゼロ交差を見つける線形補間によってその値を簡単に決定できます。
3. ベースラインの平坦化。分析化学では、湾曲したベースラインに対して吸収帯の高さを測定する必要がある場合があります。 ^[10]ベースラインの曲率は吸収帯の曲率よりもはるかに小さいため、2 次導関数はベースラインを効果的に平坦化します。吸収帯の高さに比例する導関数の高さの 3 つの尺度は、「ピークから谷まで」の距離 h1 と h2、およびベースラインからの高さ h3 です。^[11]
4. 分光法における分解能の向上。分光曲線の2次導関数のバンドはスペクトルのバンドよりも狭く、半値幅が狭くなっています。これにより、部分的に重なり合ったバンドを別々の（負の）ピークに「分解」することができます。^{[12]この図は、ピークから谷までの距離の測定を使用して、これを}化学分析にも使用できることを示しています。この場合、谷はローレンツ曲線の2次導関数の特性です。（x軸の位置は、半分の幅と半分の高さのスケール上のピーク最大値の位置に相対的です）。
5. 4 次導関数 (正のピーク) による解像度の向上。最小値は、ローレンツ関数の 4 次導関数の特性です。

「移動平均フィルター」は、短期的な変動を平滑化し、長期的な傾向やサイクルを強調するために時系列データでよく使用される Savitzky-Golay フィルターの簡単な例です。データ セットの各サブセットは、高次多項式ではなく、直線の水平線でフィッティングされます。加重なしの移動平均フィルターは、最も単純な畳み込みフィルターです。

移動平均は、株価、収益、取引量などの金融データの迅速なテクニカル分析によく使用されます。また、経済学では、国内総生産、雇用、その他のマクロ経済の時系列を調査するためにも使用されます。

すべての係数値が で同一であるため、Savitzky-Golay 畳み込み係数の一部の表には含まれませんでした。 ${\displaystyle {\frac {1}{m}}}$

データポイントが等間隔に配置されている場合、最小二乗方程式の解析解を求めることができる。 ^{[2]この解は、数値平滑化と微分化の}畳み込み法の基礎となる。データがn個の 点 ( x _j、y _j ) ( j = 1、...、n ) の集合で構成されているとする。ここで、x _jは独立変数、y _{j は}データ値である。多項式は、間隔hで区切られたm (奇数) 個の隣接するデータポイント の集合に線形最小二乗法で当てはめられる。まず、変数の変更が行われる。

${\displaystyle z={{x-{\bar {x}}} \over h}}$

ここで、は中心点の値である。zは値を取る（例：m = 5 → z = −2, −1, 0, 1, 2）。^{[注 1]} k次の多項式は次のように定義される。 ${\displaystyle {\bar {x}}}$${\displaystyle {\tfrac {1-m}{2}},\cdots ,0,\cdots ,{\tfrac {m-1}{2}}}$

${\displaystyle Y=a_{0}+a_{1}z+a_{2}z^{2}\cdots +a_{k}z^{k}.}$^{[注 2]}

係数a ₀、a _{1などは、}正規方程式を解くことによって得られます(太字のaはベクトル、太字のJ は行列を表します)。

${\displaystyle {\mathbf {a} }=\left({{\mathbf {J} }^{\mathbf {T} }{\mathbf {J} }}\right)^{-{\mathbf {1} }}{\mathbf {J} }^{\mathbf {T} }{\mathbf {y} },}$

ここで はヴァンデルモンド行列であり、 の- 行目には値があります。 ${\displaystyle \mathbf {J} }$${\displaystyle i}$${\displaystyle \mathbf {J} }$${\displaystyle 1,z_{i},z_{i}^{2},\dots }$

たとえば、5 つの点に近似された 3 次多項式z = −2、−1、0、1、2 の場合、正規方程式は次のように解かれます。

${\displaystyle \mathbf {J} ={\begin{pmatrix}1&-2&4&-8\\1&-1&1&-1\\1&0&0&0\\1&1&1&1\\1&2&4&8\end{pmatrix}}}$

${\displaystyle \mathbf {J^{T}J} ={\begin{pmatrix}m&\sum z&\sum z^{2}&\sum z^{3}\\\sum z&\sum z^{2}&\sum z^{3}&\sum z^{4}\\\sum z^{2}&\sum z^{3}&\sum z^{4}&\sum z^{5}\\\sum z^{3}&\sum z^{4}&\sum z^{5}&\sum z^{6}\\\end{pmatrix}}={\begin{pmatrix}m&0&\sum z^{2}&0\\0&\sum z^{2}&0&\sum z^{4}\\\sum z^{2}&0&\sum z^{4}&0\\0&\sum z^{4}&0&\sum z^{6}\\\end{pmatrix}}={\begin{pmatrix}5&0&10&0\\0&10&0&34\\10&0&34&0\\0&34&0&130\\\end{pmatrix}}}$

さて、正規方程式は、行と列を並べ替えることで、2つの別々の方程式に分解することができます。

${\displaystyle \mathbf {J^{T}J} _{\text{even}}={\begin{pmatrix}5&10\\10&34\\\end{pmatrix}}\quad \mathrm {and} \quad \mathbf {J^{T}J} _{\text{odd}}={\begin{pmatrix}10&34\\34&130\\\end{pmatrix}}}$

これらの各行列の逆行列の表現はクラメールの定理を使って得られる。

${\displaystyle (\mathbf {J^{T}J} )_{\text{even}}^{-1}={1 \over 70}{\begin{pmatrix}34&-10\\-10&5\\\end{pmatrix}}\quad \mathrm {and} \quad (\mathbf {J^{T}J} )_{\text{odd}}^{-1}={1 \over 144}{\begin{pmatrix}130&-34\\-34&10\\\end{pmatrix}}}$

正規方程式は

${\displaystyle {\begin{pmatrix}{a_{0}}\\{a_{2}}\\\end{pmatrix}}_{j}={1 \over 70}{\begin{pmatrix}34&-10\\-10&5\end{pmatrix}}{\begin{pmatrix}1&1&1&1&1\\4&1&0&1&4\\\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\end{pmatrix}}}$

そして

${\displaystyle {\begin{pmatrix}a_{1}\\a_{3}\\\end{pmatrix}}_{j}={1 \over 144}{\begin{pmatrix}130&-34\\-34&10\\\end{pmatrix}}{\begin{pmatrix}-2&-1&0&1&2\\-8&-1&0&1&8\\\end{pmatrix}}{\begin{pmatrix}y_{j-2}\\y_{j-1}\\y_{j}\\y_{j+1}\\y_{j+2}\\\end{pmatrix}}}$

共通の因数を掛け算して取り除くと、

${\displaystyle {\begin{aligned}a_{0,j}&={1 \over 35}(-3y_{j-2}+12y_{j-1}+17y_{j}+12y_{j+1}-3y_{j+2})\\a_{1,j}&={1 \over 12}(y_{j-2}-8y_{j-1}+8y_{j+1}-y_{j+2})\\a_{2,j}&={1 \over 14}(2y_{j-2}-y_{j-1}-2y_{j}-y_{j+1}+2y_{j+2})\\a_{3,j}&={1 \over 12}(-y_{j-2}+2y_{j-1}-2y_{j+1}+y_{j+2})\end{aligned}}}$

これらの式におけるyの係数は畳み込み係数と呼ばれ、行列の要素である。

${\displaystyle \mathbf {C=(J^{T}J)^{-1}J^{T}} }$

一般的に、

${\displaystyle (C\otimes y)_{j}\ =Y_{j}=\sum _{i=-{\tfrac {m-1}{2}}}^{\tfrac {m-1}{2}}C_{i}\,y_{j+i},\qquad {\frac {m+1}{2}}\leq j\leq n-{\frac {m-1}{2}}}$

この例は行列表記では次のように表される。

${\displaystyle {\begin{pmatrix}Y_{3}\\Y_{4}\\Y_{5}\\\vdots \end{pmatrix}}={1 \over 35}{\begin{pmatrix}-3&12&17&12&-3&0&0&\cdots \\0&-3&12&17&12&-3&0&\cdots \\0&0&-3&12&17&12&-3&\cdots \\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\vdots &\ddots \end{pmatrix}}{\begin{pmatrix}y_{1}\\y_{2}\\y_{3}\\y_{4}\\y_{5}\\y_{6}\\y_{7}\\\vdots \end{pmatrix}}}$

1964年に、Savitzky–Golay平滑化フィルタの畳み込み係数の表が、 mが25まで同じ方法で計算されて公開されました。 ^{[3] [5]}中心点の値z = 0は、係数の単一のセットから得られます。0は_平滑化、 1次導関数は1などです。数値導関数は、 Yを_微分することによって得られます。つまり、導関数は平滑化されたデータ曲線に対して計算されます。3次多項式の場合

${\displaystyle {\begin{aligned}Y&=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}=a_{0}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {dY}{dx}}&={\frac {1}{h}}\left({a_{1}+2a_{2}z+3a_{3}z^{2}}\right)={\frac {1}{h}}a_{1}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {d^{2}Y}{dx^{2}}}&={\frac {1}{h^{2}}}\left({2a_{2}+6a_{3}z}\right)={\frac {2}{h^{2}}}a_{2}&{\text{ at }}z=0,x={\bar {x}}\\{\frac {d^{3}Y}{dx^{3}}}&={\frac {6}{h^{3}}}a_{3}\end{aligned}}}$

一般に、次数(0と1)、^[注3] (2と3)、(4と5)などの多項式は、平滑化と偶数導関数に対して同じ係数を与えます。次数(1と2)、(3と4)などの多項式は、奇数導関数に対して同じ係数を与えます。

必ずしもサビツキー・ゴレイ表を使用する必要はありません。行列J ^T Jの合計は閉じた形で評価できます。

${\displaystyle {\begin{aligned}\sum _{z=-{\frac {m-1}{2}}}^{\frac {m-1}{2}}z^{2}&={m(m^{2}-1) \over 12}\\\sum z^{4}&={m(m^{2}-1)(3m^{2}-7) \over 240}\\\sum z^{6}&={m(m^{2}-1)(3m^{4}-18m^{2}+31) \over 1344}\end{aligned}}}$

畳み込み係数の代数式を導出することができる。^{[13] [注 4]}変曲点を持つ曲線に適した関数は以下のとおりである。

平滑化、多項式次数2,3: ( iの値の範囲は以下の式にも適用されます)${\displaystyle C_{0i}={\frac {\left({3m^{2}-7-20i^{2}}\right)/4}{m\left({m^{2}-4}\right)/3}};\quad {\frac {1-m}{2}}\leq i\leq {\frac {m-1}{2}}}$

1次導関数: 多項式次数3,4${\displaystyle C_{1i}={\frac {5\left({3m^{4}-18m^{2}+31}\right)i-28\left({3m^{2}-7}\right)i^{3}}{m\left({m^{2}-1}\right)\left({3m^{4}-39m^{2}+108}\right)/15}}}$

2次導関数: 多項式次数2,3${\displaystyle C_{2i}={\frac {12mi^{2}-m\left({m^{2}-1}\right)}{m^{2}\left({m^{2}-1}\right)\left({m^{2}-4}\right)/30}}}$

3次導関数: 多項式次数3,4${\displaystyle C_{3i}={\frac {-\left({3m^{2}-7}\right)i+20i^{3}}{m\left({m^{2}-1}\right)\left({3m^{4}-39m^{2}+108}\right)/420}}}$

変曲点を持たない曲線で使用できるより単純な式は次のとおりです。

平滑化、多項式次数 0,1 (移動平均):${\displaystyle C_{0i}={\frac {1}{m}}}$

1次導関数、多項式次数1,2:${\displaystyle C_{1i}={\frac {i}{m(m^{2}-1)/12}}}$

より高次の導関数を得ることもできる。例えば、2次導関数を2回実行することで4次導関数を得ることができる。^[14]

補助変数zの単純な多項式を使用してm 個のデータ ポイントを近似する代わりに、直交多項式を使用することもできます。

${\displaystyle Y=b_{0}P^{0}(z)+b_{1}P^{1}(z)\cdots +b_{k}P^{k}(z).}$

ここで、P ⁰、...、  P ^{k は}、次数 0、...、  kの相互に直交する多項式の集合である。直交多項式の表現方法と係数bとaの関係の詳細は、ゲストによって与えられている。^[2] 畳み込み係数の表現は、正規方程式行列J ^T Jが対角行列であるため簡単に得られる。これは、任意の 2 つの直交多項式の積が相互直交性によりゼロになるからである。したがって、その逆行列の非ゼロ要素はそれぞれ、正規方程式行列の対応する要素の逆数にすぎない。計算は、再帰を使用して直交グラム多項式を構築することでさらに簡略化される。計算全体は、再帰を含む計算に適したコンピュータ言語であるPASCALの数行でコード化できる。 ^[15]

Savitzky-Golay フィルタは、畳み込み係数の単一セットを使用して、中心点z = 0 での平滑化値または微分値を取得するために最も一般的に使用されます。このプロセスでは、系列の最初と最後の ( m  − 1)/2 点を計算することはできません。この不便さを回避するために、さまざまな戦略を採用できます。

• 逆の順序で、最初の ( m  − 1)/2 点のコピーを先頭に追加し、最後の ( m  − 1)/2 点のコピーを末尾に追加することで、データを人工的に拡張できます。たとえば、m = 5 の場合、データy ₁、...、y _n の先頭と末尾に 2 つの点が追加されます。

y ₃、y ₂、 y ₁、...、 y _n、 y _{n −1}、y _{n −2}。

• フィッティング多項式をもう一度見てみると、単一の多項式 a ₀ .. a _{kのすべての畳み込み係数セットを使用して、} zのすべての値に対してデータを計算できることは明らかです。

3次多項式の場合

${\displaystyle {\begin{aligned}Y&=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}\\{\frac {dY}{dx}}&={\frac {1}{h}}(a_{1}+2a_{2}z+3a_{3}z^{2})\\{\frac {d^{2}Y}{dx^{2}}}&={\frac {1}{h^{2}}}(2a_{2}+6a_{3}z)\\{\frac {d^{3}Y}{dx^{3}}}&={\frac {6}{h^{3}}}a_{3}\end{aligned}}}$

• 欠落している最初の点と最後の点の畳み込み係数も簡単に得ることができます。^[15]これは、最初の( m  +1)/2点を同じ多項式でフィッティングすることと同等であり、最後の点についても同様です。

上記の処理では、データポイントはすべて等しい重み付けがされていると暗黙的に示されています。技術的には、目的関数

${\displaystyle U=\sum _{i}w_{i}(Y_{i}-y_{i})^{2}}$

最小二乗法で最小化されるものは単位重みw _i  = 1を持つ。重みがすべて同じでない場合、正規方程式は

${\displaystyle \mathbf {a} =\left(\mathbf {J^{T}W} \mathbf {J} \right)^{-1}\mathbf {J^{T}W} \mathbf {y} \qquad W_{i,i}\neq 1}$、

すべてのデータサブセットに同じ対角重みのセットが使用される場合、正規方程式の解析解を書き出すことができます。たとえば、2次多項式では、 ${\displaystyle W={\text{diag}}(w_{1},w_{2},...,w_{m})}$

${\displaystyle \mathbf {J^{T}WJ} ={\begin{pmatrix}\sum w_{i}~~~&\sum w_{i}z_{i}&\sum w_{i}z_{i}^{2}\\\sum w_{i}z_{i}&\sum w_{i}z_{i}^{2}&\sum w_{i}z_{i}^{3}\\\sum w_{i}z_{i}^{2}&\sum w_{i}z_{i}^{3}&\sum w_{i}z_{i}^{4}\end{pmatrix}}}$

この行列の逆行列の明示的な表現はクラメールの規則を使って得ることができる。畳み込み係数の集合は次のように導出される。

${\displaystyle \mathbf {C} =\left(\mathbf {J^{T}W} \mathbf {J} \right)^{-1}\mathbf {J^{T}W} .}$

あるいは、係数Cは、スプレッドシートで計算され、組み込みの行列反転ルーチンを使用して正規方程式行列の逆行列を取得することもできます。この係数セットは、一度計算して保存すると、同じ重み付けスキームが適用されるすべての計算で使用できます。異なる重み付けスキームごとに異なる係数セットが必要です。

サビツキー・ゴレイフィルタは、フィッティング区間の両端で減少する重みを導入することで改善できることが示された。^[16]

2 次元の平滑化と微分化は、ピクセルの長方形のグリッドで構成される写真画像の強度値などのデータ値のテーブルにも適用できます。^{[17] [18]}このようなグリッドはカーネルと呼ばれ、カーネルを構成するデータ ポイントはノードと呼ばれます。秘訣は、ノードのインデックスを単純に順序付けることで、長方形のカーネルを 1 行に変換することです。1 次元のフィルター係数は、補助変数zの多項式をm 個のデータ ポイントのセットに当てはめることによって求められますが、2 次元の係数は、補助変数vとw の多項式をm × nカーネル ノードの値のセットに当てはめることによって求められます。次の例は、合計次数 3、 m = 7、n = 5の 2 変量多項式で、上記の 1 次元の場合のプロセスと並行するプロセスを示しています。^[19]

${\displaystyle v={\frac {x-{\bar {x}}}{h(x)}};w={\frac {y-{\bar {y}}}{h(y)}}}$

${\displaystyle Y=a_{00}+a_{10}v+a_{01}w+a_{20}v^{2}+a_{11}vw+a_{02}w^{2}+a_{30}v^{3}+a_{21}v^{2}w+a_{12}vw^{2}+a_{03}w^{3}}$

35個のデータ値の長方形カーネル、d ₁  −  d ₃₅

行を次々に配置するとベクトルになります。

d = ( d ₁ ... d ₃₅ ) ^T

ヤコビアンには、パラメータa ₀₀  −  a ₀₃のそれぞれに対応する10列と、vとwの値のペアに対応する35行があります。各行は次の形式です。

${\displaystyle J_{\text{row}}=1\quad v\quad w\quad v^{2}\quad vw\quad w^{2}\quad v^{3}\quad v^{2}w\quad vw^{2}\quad w^{3}}$

畳み込み係数は次のように計算される。

${\displaystyle \mathbf {C} =\left(\mathbf {J} ^{T}\mathbf {J} \right)^{-1}\mathbf {J} ^{T}}$

Cの最初の行には35 個の畳み込み係数が含まれており、それぞれ 35 個のデータ値と乗算して多項式係数 を得ることができます。これは、カーネルの中央ノード (つまり、上記の表の 18 番目のノード) での平滑化された値です。同様に、Cの他の行に 35 個の値を乗算して他の多項式係数を取得できます。この係数を使用して、異なるノードでの平滑化された値と異なる平滑化された偏導関数を取得できます。 ${\displaystyle a_{00}}$

ニキタスとパパ・ルイジは、使用される多項式の形式に応じて、平滑化の品質が大幅に変わる可能性があることを示した。^[20]彼らは、次の形式の多項式を使用することを推奨している 。

${\displaystyle Y=\sum _{i=0}^{p}\sum _{j=0}^{q}a_{ij}v^{i}w^{j}}$

このような多項式はカーネルの中心領域と境界付近の両方で良好な平滑化を実現できるため、サンプリングされた領域の内部データ ポイントと境界付近のデータ ポイントの両方で平滑化を行う際に自信を持って使用できます。最小二乗問題を解くときに悪条件を回避するには、p < mおよびq < nです。2 次元係数を計算するソフトウェアと、このようなCのデータベースについては、以下の多次元畳み込み係数のセクションを参照してください。

2次元畳み込み係数の考え方は、カーネルノードの多次元分布を1行に配置することで、直接的に高次元空間にも拡張できます^{[17] [21]}。2次元の場合における前述のNikitasとPappa-Louisiによる発見^[20]に従って、多次元の場合は次の形式の多項式の使用が推奨されます。

${\displaystyle Y=\sum _{i_{1}=0}^{p_{1}}\sum _{i_{2}=0}^{p_{2}}\cdots \sum _{i_{D}=0}^{p_{D}}(a_{i_{1}i_{2}\cdots i_{D}}\times u_{1}^{i_{1}}u_{2}^{i_{2}}\cdots u_{D}^{i_{D}})}$

ここで、Dは空間の次元、は多項式の係数、uは異なる空間方向の座標です。混合であろうとなかろうと、任意の次数の偏微分に対する代数式は、上記の式から簡単に導くことができます。^[21]ヤコビアンを作成するときに、C はカーネルノードが一列に配置される方法と、上記の多項式の展開形式のさまざまな項が配置される方法に依存することに注意してください。 ${\displaystyle a}$

コンピュータプログラミング言語で利用できる標準的な浮動小数点数の精度がもはや十分ではなくなったため、多次元の場合のCの正確な計算は困難になっています。精度が不十分なため、浮動小数点の切り捨て誤差がC の一部の要素の大きさに匹敵するようになり、その結果、精度が著しく低下して役に立たなくなります。Chandra Shekhar は、これらの精度の問題を適切に処理する 2 つのオープンソースソフトウェア、Advanced Convolution Coefficient Calculator (ACCC) と Precise Convolution Coefficient Calculator (PCCC)を発表しました。ACCCは、浮動小数点数を使用して反復的に計算を実行します。^{[22]浮動小数点数の精度は、} GNU MPFR を使用して、各反復で徐々に増加します。2 つの連続する反復で取得されたCが、事前に指定された距離まで同じ有効桁数を持ち始めると、収束に達したと見なされます。距離が十分に大きい場合、計算により非常に正確なC が生成されます。 PCCCは、 GNU Multiple Precision Arithmetic Libraryを使用して有理数計算を行い、有理数形式で完全に正確なCを生成します。^[23]最終的に、これらの有理数は、事前に指定された有効桁数になるまで浮動小数点数に変換されます。

1、2、3、4次元空間の単位間隔カーネルノード上の対称カーネルお​​よび対称多項式と非対称多項式の両方についてACCCを使用して計算されたCのデータベースが公開されています。^{[24] Chandra Shekharはまた、単位間隔カーネルノードで計算された}Cを使用して、不均一間隔カーネルノードでフィルタリングと偏微分（さまざまな順序）を実行することを記述する数学的フレームワークを提示し、 ^{[21]前述のデータベースで提供される}Cの使用を可能にしています。この方法では近似結果しか得られませんが、カーネルノードの非均一性が弱い限り、ほとんどのエンジニアリングアプリケーションで受け入れられます。

1. 平滑化のための畳み込み係数の合計は1に等しい。奇数導関数の係数の合計は0である。^[25]
2. 平滑化のための畳み込み係数の二乗和は中心係数の値に等しい。^[26]
3. 関数を平滑化しても関数の下の領域は変化しない。^[25]
4. 対称関数を偶数階微分係数で畳み込むと対称中心が保存される。^[25]
5. 微分フィルタの特性。^[27]

畳み込み処理で信号が歪むのは避けられません。上記の性質3から、ピークを持つデータを平滑化すると、ピークの高さは下がり、半値幅は広がります。歪みの程度とS/N（信号対雑音比）はともに向上します。

• 多項式の次数が増加するにつれて減少する
• 畳み込み関数の幅mが増加するにつれて増加する

例えば、すべてのデータポイントのノイズが無相関で、一定の標準偏差σを持つ場合、ノイズの標準偏差はmポイントの平滑化関数との畳み込みによって減少します^{[26] [注5]}

多項式次数0または1: (移動平均)${\displaystyle {\sqrt {1 \over m}}\sigma }$

多項式の次数2または3: .${\displaystyle {\sqrt {\frac {3(3m^{2}-7)}{4m(m^{2}-4)}}}\sigma }$

これらの関数は、右側のグラフに示されています。たとえば、9 ポイントの線形関数 (移動平均) ではノイズの 3 分の 2 が除去されますが、9 ポイントの 2 次/3 次平滑化関数ではノイズの約半分しか除去されません。残っているノイズのほとんどは低周波ノイズです (以下の「畳み込みフィルタの周波数特性」を参照)。

移動平均関数はノイズ低減に優れていますが、m点以上の曲率を持つデータを平滑化するには適していません。2次フィルタ関数は、2次多項式には変曲点がないため、変曲点を持つデータ曲線の導関数を得るのには適していません。多項式の次数と畳み込み係数の数の最適な選択は、ノイズ低減と歪みの間の妥協点になります。^[28]

歪みを軽減し、ノイズ除去を改善する1つの方法は、より狭い幅のフィルタを使用し、それに対して複数回の畳み込みを実行することです。同じフィルタを2回通過させると、元のフィルタとそれ自身との畳み込みによって得られたフィルタを1回通過させるのと同等になります。^[29]たとえば、係数(1/3、1/3、1/3)を持つフィルタを2回通過させることは、係数(1/9、2/9、3/9、2/9、1/9)を持つフィルタを1回通過させることと同等です。

マルチパスの欠点は、-ポイント関数のパスの等価フィルタ幅が広いため、マルチパスではより大きなエンドエフェクトの影響を受けることです。それでも、マルチパスは大きな利点があります。たとえば、信号対雑音比がわずか5のデータに対して40～80回のパスを行うと、有用な結果が得られました。^[30] 上記のノイズ低減式は適用されません。計算されたデータポイント間の相関はパスごとに増加するためです。 ${\displaystyle n}$${\displaystyle m}$${\displaystyle n(m-1)+1}$

畳み込みはフーリエ 変換領域における乗算にマッピングされます。畳み込みフィルタの離散フーリエ変換は、次のように表される 実数値関数です。

${\displaystyle FT(\theta )=\sum _{j={\tfrac {1-m}{2}}}^{\tfrac {m-1}{2}}C_{j}\cos(j\theta )}$

θ は 0 度から 180度まで実行され、その後は関数が単に繰り返されます。9 ポイントの 2 次/3 次平滑化関数のプロットは典型的です。非常に低い角度では、プロットはほぼ平坦です。つまり、データの低周波成分は平滑化操作によって実質的に変化しません。角度が増加すると、値が減少するため、高周波成分はますます減衰します。これは、畳み込みフィルターがローパスフィルターとして記述できることを示しています。つまり、除去されるノイズは主に高周波ノイズであり、低周波ノイズはフィルターを通過します。^[31]大きな角度でのフーリエ変換のうねりによって示されるように、一部の高周波ノイズ成分は他の成分よりも減衰します。これにより、平滑化されたデータに小さな振動が発生し^[32]、位相反転、つまりデータ内の高周波振動が Savitzky–Golay フィルタリングによって反転される可能性があります。^[33]

畳み込みは、データ内のエラー間の相関に影響します。畳み込みの効果は線形変換として表現できます。

${\displaystyle Y_{j}\ =\sum _{i=-{\frac {m-1}{2}}}^{\frac {m-1}{2}}C_{i}\,y_{j+i}}$

誤差伝播の法則により、データAの分散共分散行列は次のように B に変換される。

${\displaystyle \mathbf {B} =\mathbf {C} \mathbf {A} \mathbf {C} ^{T}}$

これが実際にどのように適用されるかを確認するには、データ ポイントの分散が等しく、データ ポイント間に相関関係がないと仮定して、計算された最初の 3 つのポイントY ₂  −  Y ₄に対する 3 ポイント移動平均の効果を考えてみましょう。Aは、各ポイントの分散である定数σ ^{2を乗じた}単位行列になります。

${\displaystyle \mathbf {B} ={\sigma ^{2} \over 9}{\begin{pmatrix}1&1&1&0&0\\0&1&1&1&0\\0&0&1&1&1\end{pmatrix}}{\begin{pmatrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1&0&0\\0&0&0&1&0\\0&0&0&0&1\end{pmatrix}}{\begin{pmatrix}1&0&0\\1&1&0\\1&1&1\\0&1&1\\0&0&1\end{pmatrix}}={\sigma ^{2} \over 9}{\begin{pmatrix}3&2&1\\2&3&2\\1&2&3\end{pmatrix}}}$

この場合、相関係数は、

${\displaystyle \rho _{ij}={\frac {B_{ij}}{\sqrt {B_{ii}B_{jj}}}},\ (i\neq j)}$

計算された点iとjの間の値は

${\displaystyle \rho _{i,i+1}={2 \over 3}=0.66,\,\rho _{i,i+2}={1 \over 3}=0.33}$

一般的に、観測値が相関していなくても計算値は相関している。相関は一度にm  −1個の計算点に及ぶ。 ^[34]

マルチパスがデータセットのノイズと相関に与える影響を説明するために、3点移動平均フィルタの2回目のパスの影響を考えてみましょう。2回目のパス^[注6]

${\displaystyle \mathbf {CBC^{T}} ={\sigma ^{2} \over 81}{\begin{pmatrix}1&1&1&0&0\\0&1&1&1&0\\0&0&1&1&1\end{pmatrix}}{\begin{pmatrix}3&2&1&0&0\\2&3&2&0&0\\1&2&3&0&0\\0&0&0&0&0\\0&0&0&0&0\end{pmatrix}}{\begin{pmatrix}1&0&0\\1&1&0\\1&1&1\\0&1&1\\0&0&1\end{pmatrix}}={\sigma ^{2} \over 81}{\begin{pmatrix}19&16&10&4&1\\16&19&16&10&4\\10&16&19&16&10\\4&10&16&19&16\\1&4&10&16&19\end{pmatrix}}}$

2 回のパスの後、中心点の標準偏差は に減少し、1 回のパスの 0.58 σと比較して減少しています。ノイズの減少は、同じ条件下で 5 ポイント移動平均の 1 回のパスで得られるノイズの減少よりもわずかに少なく、平滑化されたポイントの標準偏差は 0.45 σと小さくなります。 ${\displaystyle {\sqrt {\tfrac {19}{81}}}\sigma =0.48\sigma }$

相関関係は、相関係数を持つ4つの連続したポイントの範囲にまで広がります。

${\displaystyle \rho _{i,i+1}={16 \over 19}=0.84,\rho _{i,i+2}={10 \over 19}=0.53,\rho _{i,i+3}={4 \over 19}=0.21,\rho _{i,i+4}={1 \over 19}=0.05}$

より狭いスムージング関数を使用して 2 回のパスを実行することによって得られる利点は、計算されたデータに生じる歪みが少なくなることです。

他の平滑化フィルタ、例えばガウス分布との畳み込みやマルチパス移動平均フィルタリングと比較すると、Savitzky-Golay フィルタは、特に近似多項式の高次数において、周波数領域での初期応答がより平坦でカットオフがより急峻です (周波数特性を参照)。信号帯域幅が制限されたデータの場合、これは Savitzky-Golay フィルタリングが他の多くのフィルタよりも優れた信号対雑音比を提供できることを意味します。例えば、スペクトルのピーク高さは、同様のノイズ抑制を持つ他のフィルタよりもよく保持されます。Savitzky-Golay フィルタの欠点は、一部の高周波数の抑制が比較的貧弱であること (阻止帯域抑制が貧弱) と、最初と最後のポイントに多項式近似を使用した場合にアーティファクトが発生することです。^[16]

サビツキー・ゴレイフィルタの利点を共有し、少なくともいくつかの欠点を軽減する代替の平滑化手法としては、適切に選択された代替フィッティングウェイトを使用したサビツキー・ゴレイフィルタ、ウィテカー・ヘンダーソン平滑化とホドリック・プレスコットフィルタ（平滑化スプラインに密接に関連する同等の手法）、およびウィンドウ付き sinc関数を使用した畳み込みがある。^[16]

• flattenモジュール Lightkurve から。Lightkurve は、Kepler および TESS 望遠鏡のデータを分析するための公式ライブラリです。
• scipy.signal.savgol_filterモジュールSciPyから。SciPy は、学術コミュニティで科学計算に広く使用されている堅牢なライブラリです。

• カーネルスムーザー–統計で使用される多くの同じプロセスに対する異なる用語
• 局所回帰- LOESS法とLOWESS法
• 数値微分–関数の微分への応用
• スムージングスプライン
• ステンシル（数値解析） - 微分方程式の解法への応用
• ホドリック・プレスコットフィルター
• カルマンフィルター

データポイントのセット⁠ ⁠を考えます。Savitzky–Golay 表は、ステップ${\displaystyle (x_{j},y_{j})_{1\leq j\leq n}}$⁠ ⁠${\displaystyle x_{j}-x_{j-1}}$が一定、h の場合を指します。3 次多項式と 5 ポイントのウィンドウサイズmを使用した、いわゆる畳み込み係数の使用例は次のとおりです。

スムージング

${\displaystyle Y_{j}={\frac {1}{35}}(-3\times y_{j-2}+12\times y_{j-1}+17\times y_{j}+12\times y_{j+1}-3\times y_{j+2})}$ ;

1次導関数

${\displaystyle Y'_{j}={\frac {1}{12h}}(1\times y_{j-2}-8\times y_{j-1}+0\times y_{j}+8\times y_{j+1}-1\times y_{j+2})}$ ;

2次導関数

${\displaystyle Y''_{j}={\frac {1}{7h^{2}}}(2\times y_{j-2}-1\times y_{j-1}-2\times y_{j}-1\times y_{j+1}+2\times y_{j+2})}$。

1、2、3、4、5次の多項式の畳み込み係数の選択された値が次の表に示されています^[注7]。値はGorryで提供されているPASCALコードを使用して計算されました。^[15]

• ガンズ、ピーター（1992）。化学科学におけるデータフィッティング：最小二乗法による。Wiley。ISBN 9780471934127。

• 多次元最小二乗フィルタ用の高度な畳み込み係数計算機 (ACCC)
• 統計学の基礎におけるサビツキー・ゴレイフィルタ
• データセットのサイズ、適合度、中心点からのオフセットの範囲に応じた幅広い係数