[2.1]スペクトル解析の前処理 (②ノイズ低減と特徴強調)
応用編：Savitzky-Golayフィルタ

スペクトル解析における前処理は、最終的な解析精度を向上させるために重要です。
[2.1]スペクトル解析の前処理 (②ノイズ低減と特徴強調 基礎編)では、「スムージングによる、測定ノイズの低減」と、「微分処理によるピーク位置強調」について解説しました。
しかし、それぞれの手法には、「 スムージングは、ピークの頂点、位置がわずかにずれる」、「微分は、ノイズまで強調してしまう」という形で、両立が難しいという課題があります。

その「ノイズ低減」と「ピーク保持」を両立する手法として、分光解析の分野で広く利用されているのが Savitzky-Golayフィルタです。

本記事では、Savitzky-Golayフィルタの「概要」から、「Pythonを持ちいたフィルタ適用例」まで幅広く解説していきます。

1. Savitzky-Golayフィルタの概要

Savitzky-Golayフィルタは、スペクトルの形状をできるだけ保ちながらノイズを低減する平滑化手法です。
分光分析などピーク位置の精度が重要となる解析において標準的な前処理として使用されています。

1.1 Savitzky-Golayフィルタの原理

Savitzky-Golayフィルタは、「局所的なウィンドウ内で多項式回帰を行う」フィルタです。

まず、「局所的なウィンドウ」とは、スペクトルデータの各点を中心として、その前後の一定範囲のデータを取り出した領域のことを指します。

スペクトルは通常、数百点からなる離散データで構成されていますが、Savitzky-Golayフィルタでは、それぞれの点ごとに周囲のデータ（例えば前後に数点ずつ）を取り出し、その範囲の情報を用いてスムージング後の値を決定します。

局所的なウィンドウの例

以下は、スペクトルの6点目および7点目を対象に、それぞれの前後3点（計7点）を局所ウィンドウとして抽出した例です。

次に「多項式回帰」とは、取り出したウィンドウ内のデータに対して、低次の多項式（例えば2次や3次）を最小二乗法によってフィッティングする処理です。

最小二乗法とは、実際のデータと近似曲線とのずれ（誤差）の二乗和が最も小さくなるように、多項式の係数を決定する方法です。
つまり、ウィンドウ内のデータ全体に最もよく一致する滑らかな曲線を求める操作といえます。
そして、このようにして得られた多項式を用いて、ウィンドウの中心点における値を評価し、その値をスムージング後のデータとして採用します。このとき重要なのは、単純に平均値を取るのではなく、「周囲のデータに最も適合する曲線上の値」を採用している点です。

多項式回帰の例

以下は、スペクトルの6点目および7点目を対象に、それぞれの局所ウィンドウ内で２次多項式近似を実行した例です。

この処理をスペクトルの全データ点に対して順番に繰り返していくことで、ノイズが低減されつつもピーク形状や局所的な変化が保持された、滑らかなスペクトルが得られます。

Savitzky-Golayフィルタの例

以下は、各点に対して前後3点ずつ（計7点）の局所ウィンドウを取り出し、2次多項式による回帰を行った結果です。

1.2 Savitzky-Golayフィルタの特徴

①スペクトル波形のピークが崩れにくい

たとえば、一般的なスムージングである”移動平均”は、局所的なウィンドウ内で平均値を取る手法です。
移動平均の手順は、Savitzky-Golayフィルタと似ていますが、平均値を取ることで頂点が丸まり、ピークがなくなったり、ずれてしまうことがあります。

一方、Savitzky-Golayフィルタは、局所的なウィンドウ内のデータを曲線として近似します。
2次項以上の項は曲率を表現することができるため、「曲率」であるピークの丸みや対称性といった形状情報を保持することができます。
これにより、ピークの高さや位置を保ちながらノイズを抑えるという両立が可能になります。

②微分によるピーク検出が可能

Savitzky-Golayフィルタのもう一つの重要な特徴は、平滑化だけでなく微分処理にも適している点です。

2次多項式スムージングは、局所ウィンドウ内のデータを2次式で近似し、その中心値を評価する手法であり、”微分処理”を行うことまでを目的として設計されているものではありません。
そのため、”微分処理”を行う場合は、スムージング後のデータに対して改めて処理をすることになり、その段階で再びノイズが増幅されてしまう可能性があります。

一方、Savitzky-Golayフィルタは”微分処理”を含めて設計された手法です。
局所多項式で近似した後、その多項式を解析的に微分します。
この方法では、単純な差分による微分と比べてノイズの増幅を抑えることができ、より安定した微分スペクトルを得ることができます。
また、微分を行うことでスペクトルの変化が強調されるため、ピーク位置の検出や近接したピークの分離が容易になるという利点があります。

2. Savitzky-Golayフィルタにおける多項式近似と微分の関係

１章では、Savitzky-Golayフィルタが「局所ウィンドウ内で多項式回帰を行うことで、スペクトルの平滑化と微分を同時に行える手法」であることを説明しました。
２章では、その仕組みをもう少し具体的に理解するため、多項式近似と微分の関係を数式を用いて解説していきます。

2.1 多項式近似における"a₀"の意味

Savitzky-Golayフィルタでは、特定の点の前後を含む「局所ウィンドウ」を取り出して、その範囲のデータを、以下のような「多項式」で近似します。

$P\left(x\right)=a_0+a_{1}x+a_2{x}^2$

この数式における "x" は、ウィンドウの中心を "0" とした相対座標です。
例えば中心の左右に3点ずつ計7点ある場合、(x=-3,-2,-1,0,1,2,3)で考えることができます。

この式に、中心点 "x=0" を代入すると、

$P\left(0\right)=a_0$

となります。
よって、Savitzky-Golayフィルタ適用後の中心点の値は、a0 です。

このことから、Savitzky-Golayフィルタは「周囲のデータに最もよく合う曲線」を求め、その曲線上の中心値を使用することがわかります。

2.2 多項式近似の 1次微分と傾き（a₁の意味）

次に、この近似式を微分すると、

$\frac{dP}{dx}=a_1+2a_{2}x$

となります。
この数式に中心点 (x=0) を代入すると、

${\frac{dP}{dx}}_{x=0}=a_1$

となります。

つまり、多項式の1次係数である $a_1$ は、中心点の1次微分、すなわち傾きを表しています。
よって、Savitzky-Golayフィルタでは、多項式の$a_1$ を利用することで、中心点における1次微分値を求めることができます。

2.3 2次微分と曲率（a₂の意味）

さらに微分して2次微分を計算します。

$\frac{d^{2}P}{d{x}^2}=2a_2$

この式から、2次係数 $a_2$ は曲率に対応しており、 $2a_2$ が2次微分値になることがわかります。

スペクトルのピークでは曲線が上に凸の形になるため、2次微分値は負になりやすくなります。 そのため、2次微分スペクトルではピーク位置が極小値として現れ、ピーク検出や近接ピークの分離に利用できます。

2.4 Savitzky-Golayフィルタが安定して微分を行える理由

ここで重要なのは、Savitzky-Golayフィルタが、生データをそのまま微分しているわけではないという点です。

通常の微分は隣接点との傾きを取るため、ノイズも同時に強調されやすくなります。
一方、Savitzky-Golayフィルタでは、まず局所データを滑らかな多項式関数として近似し、その係数から微分値を計算しています。

つまり、「離散的な測定データ」→「 滑らかな関数」→「その関数の微分」という手順で計算をおこなうため、ノイズの影響を抑えながら、安定した微分スペクトルを得ることができます。

3. Savitzky-Golayフィルタのパラメータ設定

Savitzky-Golayフィルタの性能は、パラメータ設定に大きく依存するため、パラメータの意味を知っておくことが重要です。
ここからは、「ノイズ低減」と「ピーク保持」のバランスを最適化するための以下の３つのパラメータを紹介していきます。
(1)ウィンドウ幅(window_length)
(2)多項式次数(polyorder)
(3)微分次数(deriv)

3.1 ウィンドウ幅

局所多項式近似を行う際に使用するデータ点数です。ウィンドウの中心点を明確にするためには、必ず奇数を指定する必要があります。
たとえば、中心点の前後に3点ずつを含む領域をWindowにするには、window_length = 7を設定します。
window_lengthを大きくすると、より多くの点を使って近似するため、ノイズ低減効果は強くなりますが、ピークの高さが低下したり、細かい構造が失われます。逆に小さくすると、ピーク形状は保持されますが、ノイズ抑制効果は弱くなるというトレードオフ関係です。

3.2 多項式次数

局所近似に用いる多項式の次数で、2次または3次が一般的です。
2次項があれば曲率を表現できるため、スペクトルの局所構造を十分に捉えられます。
次数を高くすると、より複雑な形状を近似できるようになりますが、データ点間で不自然な振動が起こり、微分値が不安定になることがあります。
そのため、次数は大きくても３次元程度に抑えるのが、望ましいです。
また、次数はウィンドウ幅よりも大きくすることはできません。たとえば、ウィンドウ幅を 5 とした場合、多項式次数 は最大でも4以下にする必要があります。

3.3 微分次数

出力として求める微分次数を指定するパラメータです。
Savitzky-Golayフィルタの特徴は、微分次数を変更するだけで、平滑化と微分を切り替えることができます。

・平滑化として使う場合

微分次数が 0 のときは、局所多項式の定数項 $a_0$が出力されます。
これはノイズ低減を目的とした通常の平滑化です。

・1次微分として使う場合

微分次数を 1 にすると、局所多項式の1次係数 $a_1$ に対応する値が出力されます。
ピークの変化点や傾きを調べることができます。

・2次微分として使う場合

微分次数を 2にすると、局所多項式の2次係数$a_2$ に対応する値が出力されます。
ピーク位置では2次微分が極小値を取るため、近接ピークの分離やピーク検出が可能となります。
微分次数が上がるほどノイズの影響は強くなるため、window_lengthと合わせて調整が必要になります。

3.4 パラメータのまとめ

Savitzky-Golayフィルタの主なパラメータは、ウィンドウ幅、多項式次数、微分次数の3つであり、それぞれ以下のような役割を持っています。
ウィンドウ幅：スムージングの強度
多項式次数：スムージングの際の形状再現度
微分次数：ピーク強調

これらのパラメータのバランスを適切に取ることが、分光解析において安定したピーク検出や干渉周期解析の精度向上につながります。

4. PythonでのSavitzky-Golayフィルタの実装

Savitzky-Golayフィルタは、Pythonでは SciPyライブラリのsignal.savgol_filter を使用することで簡単に実装することが可能です。
本章では signal.savgol_filter の基本的な使い方を説明し、その後にスペクトルデータへ適用する実装例を紹介します。

4.1 signal.savgol_filter 関数とは

SciPyは、Pythonで科学技術計算を行うための代表的なライブラリで、数値積分や最適化、統計処理、線形代数、信号処理、画像処理など幅広い機能が実装されています。
中でもscipy.signalは信号処理に特化したモジュールで、ノイズ除去のためのフィルタリングやフーリエ変換、スペクトル解析、畳み込み処理などの機能が利用できます。
このモジュールには、Savitzky-Golayフィルタを実現する関数「scipy.signal.savgol_filter」が含まれています。
フィルター機能を一から実装するのは非常に手間がかかりますが、このようなライブラリを使用することで、簡単にフィルター処理を実現することが可能となります。

4.2 signal.savgol_filterの使い方

signal.savgol_filterコマンドは、以下のように実行します。

filtered = signal.savgol_filter(data, window_length=11, polyorder=3, deriv=2)

この記述により、入力データ data に対してSavitzky-Golayフィルタを適用し結果を返します。

この時、「3. Savitzky-Golayフィルタのパラメータ設定」で紹介した3つのパラメータを以下のように設定することが可能です。

例えば、
signal.savgol_filter(data, window_length=11, polyorder=3, deriv=2)
と設定すると、各点に対して、前後を含め11点のデータを用いて3次多項式で局所近似を行い、その2次微分の計算結果が戻り値になります。

4.3 スペクトルデータへの適用例

次に、実際のスペクトルデータにSavitzky-Golayフィルタを適用する例を紹介します。

ここではスペクトルデータを読み込み、Savitzky-Golayフィルタを用いて2次微分スペクトルを計算し、元スペクトルと比較表示するサンプルコードを紹介します。

サンプルコード

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from spectral import envi

# --- (1) スペクトルイメージ ENVI (.hdr) の読み込み 
img = envi.open("your_file.hdr") 
cube = img.load() # shape: (H, W, B)

# --- (2)スペクトルイメージから任意ピクセルのスペクトル取得 
H, W, B = cube.shape
y = H // 2
x = W // 2

spectrum = cube[y, x, :] # shape: (B,)

# hdrから波長を取得
wavelength = np.array(img.metadata['wavelength']).astype(float)

# --- (3) Savitzky-Golayフィルタの適用 
second_derivative = signal.savgol_filter(
spectrum,
window_length=11,
polyorder=3,
deriv=2
)

# --- (4) Savitzky-Golayフィルタ適用後のグラフ表示 ---
fig, ax1 = plt.subplots(figsize=(8,6))

# 元スペクトル（左軸）
ax1.plot(wavelength, spectrum, color="blue", label="Original Spectrum")
ax1.set_xlabel("Wavelength")
ax1.set_ylabel("Intensity", color="blue")

# 2次微分スペクトル（右軸）
ax2 = ax1.twinx()
ax2.plot(wavelength, second_derivative, color="red", label="2nd Derivative")
ax2.set_ylabel("2nd Derivative", color="red")

plt.tight_layout()
plt.show()

(1) スペクトルイメージ ENVI (.hdr) の読み込み

ENVI形式のハイパースペクトルデータを読み込む処理です。
.hdr ファイルにはデータのサイズや波長情報などのメタデータが記録されており、Spectral Python (SPy) を使うことで対応するデータファイルを読み込むことが可能です。
img.load() によってデータ本体が読み込まれ、(H, W, B) の形状を持つ3次元配列として取得します。

関連記事

Pythonでハイパースペクトルデータを解析するためには、まずスペクトルデータを正しく読み込む必要があります。
詳しい読み込み方法については、以下の記事をご参照ください。

→ 【Python×スペクトル解析】ハイパースペクトルデータの読み込み方法

(2)スペクトルイメージから任意ピクセルのスペクトル取得

H, W, B = cube.shape で、読み込んだデータのサイズ（高さ H、幅 W、バンド数 B ）を取得し、中央付近の画素を選択するために H // 2 と W // 2 で高さと幅の真ん中の値を抽出しました。 最後に、spectrum = cube[y, x, :]で、指定した座標における各バンドの輝度値を１次元配列で抽出します。

関連記事

読み込んだスペクトルデータを画像表示して、クリックした点のスペクトルを取得する方法に関しては、以下の記事をご参照ください。

→ 【Python×スペクトル解析】スペクトルデータのグラフ表示

(3) Savitzky-Golayフィルタの適用

前節で説明した SciPyのsignal.savgol_filter を用いて、取得したスペクトルデータ spectrum にSavitzky-Golayフィルタを適用しました。
今回のフィルターでは、周囲11点のスペクトルデータを滑らかな3次元の曲線で近似し、その曲線の曲率を求めることで2次微分スペクトルを算出しています。

(4) Savitzky-Golayフィルタ適用後のグラフ表示

元スペクトルとSavitzky-Golayフィルタによる2次微分スペクトルを、1つのグラフ上に重ねて表示します。
左軸(ax1)に元スペクトルを描画し、右軸(ax2)にフィルタ後のスペクトルを重ねます。これにより、ピーク位置と微分結果の関係を同時に確認できます。
それぞれの軸ラベルを設定し、tight_layout() でレイアウトを整えた後、show() によって図を表示します。

delta = wavelength[1] - wavelength[0]

# Savitzky-Golayフィルタ（2次微分）
    second_derivative = signal.savgol_filter(
        spectrum,
        window_length=11,
        polyorder=3,
        deriv=2,
        delta=delta
    )

5. まとめ

本記事では、スペクトル解析において非常によく使用される前処理手法として、Savitzky-Golayフィルタについて解説しました。

Savitzky-Golayフィルタは、指定したサイズの局所ウィンドウ内で多項式近似を行うことで、スペクトルの形状を保ちながらノイズを低減することができます。
また、スムージングと微分を順に行った際に、微分でノイズが増幅されやすいという課題がありますが、Savitzky-Golayフィルタはそれらを同時に行うので、スムージングとピーク強調をバランスよく実現することができます。

Pythonでは、SciPyの signal.savgol_filter を用いることで簡単に実装できるため、実務においても非常に扱いやすい手法です。
パラメータである「ウィンドウ幅」「多項式次数」「微分次数」を適切に設定することで、ノイズ低減とピーク保持のバランスを調整できるため、目的に応じた柔軟な前処理が可能となります。

次の記事

前処理によって整えられたスペクトルデータは、その後の統計解析や機械学習によって「分類」や「特徴抽出」に活用されます。
次の記事では、主成分分析（PCA）を用いてスペクトルデータを可視化・分類する方法について解説します。

→ 【Python×スペクトル解析】主成分分析（PCA）による分類入門

このカテゴリのコンテンツ

2.1 前処理 ③ノイズ低減/特徴強調(2)

Pythonで始めるスペクトル解析

1.1 スペクトルデータの読み込み

1.2 スペクトルデータ画像表示

1.3 スペクトルグラフ表示

1.4 スペクトルビューア作成

関連コンテンツ

スペクトルカメラのインライン導入

論文・研究報告（LemnaTec）

低価格マルチスペクトルカメラの実力

UAV×ハイパー実機撮影デモ

ハイパースペクトルカメラの知識や撮影事例に関する情報がたくさん！

KLV大学
ハイパースペクトルカメラコース