SciPy で Savitzky-Golay フィルタ

SciPy での Savitzky-Golay フィルタについてメモ。
Savitzky-Golay フィルタは最小二乗法による多項式近似により信号を平準化します。
信号の高周波数成分を維持しつつ平準化したいときに効果的なフィルタです。

特定の周波数帯域のノイズを除去する場合は通常の FIR フィルタの方が優れています。

基本的には savgol_filter 関数に対象の信号、フィルタ窓の長さ(近似に使うデータの個数)、多項式の次数を指定すれば平準化された信号を取得できます。

また、オブション引数の mode で信号にパディングする値の拡張タイプを選択できます。

フィルタ係数を求めたいだけなら savgol_coeffs を使用します。

# coding: utf-8
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# 時系列のサンプルデータ作成
n = 512                         # データ数
dt = 0.01                       # サンプリング間隔
f = 1                           # 周波数
t = np.linspace(1, n, n)*dt-dt
y = np.sin(2*np.pi*f*t)+0.5*np.random.randn(t.size)

y1 = signal.savgol_filter(y, n/4+1, 5)
y2 = signal.savgol_filter(y, n/4+1, 5, mode="nearest")
y3 = signal.savgol_filter(y, n/4+1, 5, mode="mirror")
y4 = signal.savgol_filter(y, n/4+1, 5, mode="constant")
y5 = signal.savgol_filter(y, n/4+1, 5, mode="wrap")

plt.figure(figsize=(12,9))
plt.plot(t, y)
plt.plot(t, y1, "m", linewidth=2, label="interp")
plt.plot(t, y2, "r", linewidth=2, label="nearest")
plt.plot(t, y3, "c", linewidth=2, label="mirror")
plt.plot(t, y4, "y", linewidth=2, label="constant")
plt.plot(t, y5, "k", linewidth=2, label="wrap")
plt.axis("tight")
plt.legend(loc="upper right")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.show()

15071801.png

コメント

Comments powered by Disqus
書籍更新情報
2016-10-21
Pythonによる科学技術計算 基礎編
PDF版の販売を開始しました。
販売ページはこちら

2016-09-09
Pythonによる科学技術計算 基礎編
1.2版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
Kindle ストア、Leanpubで販売中です
Pythonによる科学技術計算 基礎編
PDF版の販売はこちら
同人誌
技術書典(2016/6/25)
Emacs/org-modeのPDF作成術
電子版をBOOTHで販売中です
Emacs/org-modeのPDF作成術
Share