サビツキー・ゴーレイフィルタの理論とPython実装

import numpy as np
from math import factorial
from scipy.signal import savgol_coeffs

def sg_coeffs_manual(window_length, polyorder, deriv=0, delta=1.0):
    m = (window_length - 1) // 2
    n = np.arange(-m, m + 1)
    # Vandermonde行列: A[i, k] = n_i^k  (式3と同じ定義)
    A = np.vander(n, N=polyorder + 1, increasing=True)
    # 正規方程式の解行列 (A^T A)^{-1} A^T  (式4)
    AtA_inv_At = np.linalg.inv(A.T @ A) @ A.T
    # d次微分係数は d! 倍した d 行目 (式6)
    row = AtA_inv_At[deriv, :]
    coeffs = factorial(deriv) * row / (delta ** deriv)
    return coeffs

# window=5, p=2 の平滑化係数（手計算した値と比較）
c_5_2 = sg_coeffs_manual(5, 2, deriv=0)
print("window=5, p=2:", np.round(c_5_2, 6))
# scipyのデフォルト(use='conv')ではなく use='dot' で比較する
# （理由は次の実行結果の後で説明）
c_scipy_dot = savgol_coeffs(5, 2, deriv=0, use='dot')
print("scipy (use=dot):", np.round(c_scipy_dot, 6))

# window/order/derivを変えた網羅的な比較
cases = [(5,2,0), (7,2,0), (9,3,0), (11,3,0), (11,4,0),
         (7,2,1), (9,3,1), (11,4,1),
         (9,4,2), (11,5,2)]
max_diff = 0.0
for w, p, d in cases:
    manual = sg_coeffs_manual(w, p, deriv=d)
    scipy_dot = savgol_coeffs(w, p, deriv=d, use='dot')
    max_diff = max(max_diff, np.max(np.abs(manual - scipy_dot)))
print("全ケースでの最大差 (use='dot'):", max_diff)

# デフォルト設定(use='conv')との比較も試す
max_diff_conv = 0.0
for w, p, d in cases:
    manual = sg_coeffs_manual(w, p, deriv=d)
    scipy_conv = savgol_coeffs(w, p, deriv=d)  # デフォルト use='conv'
    max_diff_conv = max(max_diff_conv, np.max(np.abs(manual - scipy_conv)))
print("全ケースでの最大差 (デフォルトuse='conv'):", max_diff_conv)

import numpy as np
from scipy.signal import savgol_filter
import matplotlib.pyplot as plt

# テストデータ生成
np.random.seed(42)
t = np.linspace(0, 1, 500)
# ピークを含む信号 + ノイズ
signal_clean = np.sin(2 * np.pi * 3 * t) + 0.5 * np.exp(-((t - 0.3) / 0.02)**2)
noise = 0.2 * np.random.randn(len(t))
signal_noisy = signal_clean + noise

# サビツキー・ゴーレイフィルタ適用
# window_length: 窓幅（奇数）, polyorder: 多項式次数
y_sg = savgol_filter(signal_noisy, window_length=21, polyorder=3)

# 移動平均との比較
kernel = np.ones(21) / 21
y_ma = np.convolve(signal_noisy, kernel, mode='same')

# プロット
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

axes[0].plot(t, signal_clean, 'g', label='Original', linewidth=2)
axes[0].plot(t, signal_noisy, 'gray', alpha=0.5, label='Noisy')
axes[0].plot(t, y_sg, 'b', label='Savitzky-Golay (w=21, p=3)')
axes[0].plot(t, y_ma, 'r--', label='Moving Average (w=21)')
axes[0].set_title('Smoothing Comparison')
axes[0].legend()
axes[0].grid(True)

# ピーク付近の拡大
mask = (t > 0.2) & (t < 0.4)
axes[1].plot(t[mask], signal_clean[mask], 'g', label='Original', linewidth=2)
axes[1].plot(t[mask], y_sg[mask], 'b', label='Savitzky-Golay')
axes[1].plot(t[mask], y_ma[mask], 'r--', label='Moving Average')
axes[1].set_title('Peak Region (Zoomed)')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()

import numpy as np
from scipy.signal import savgol_filter

t = np.linspace(0, 1, 500)
# 非周期・非平坦な信号（境界でも傾き/曲率がゼロでない）
signal_clean = t**2 + 0.3 * np.sin(2 * np.pi * 5 * t)

window_length = 21
polyorder = 3
m = window_length // 2  # 10

modes = ['mirror', 'nearest', 'constant', 'wrap', 'interp']
n_trials = 200
noise_sigma = 0.15

boundary_rmse = {mode: [] for mode in modes}
interior_rmse = {mode: [] for mode in modes}

rng = np.random.default_rng(7)
for _ in range(n_trials):
    noise = noise_sigma * rng.standard_normal(len(t))
    signal_noisy = signal_clean + noise
    for mode in modes:
        y = savgol_filter(signal_noisy, window_length, polyorder, mode=mode)
        err_boundary = np.concatenate([y[:m] - signal_clean[:m], y[-m:] - signal_clean[-m:]])
        err_interior = y[m:-m] - signal_clean[m:-m]
        boundary_rmse[mode].append(np.sqrt(np.mean(err_boundary**2)))
        interior_rmse[mode].append(np.sqrt(np.mean(err_interior**2)))

print(f"{'mode':>10} {'boundary RMSE':>15} {'interior RMSE':>15}")
for mode in modes:
    print(f"{mode:>10} {np.mean(boundary_rmse[mode]):15.5f} {np.mean(interior_rmse[mode]):15.5f}")

# 境界の最初のサンプル(t=0, 真値=0)におけるバイアスと分散を分解
first_vals = {mode: [] for mode in modes}
rng = np.random.default_rng(7)
for _ in range(n_trials):
    noise = noise_sigma * rng.standard_normal(len(t))
    signal_noisy = signal_clean + noise
    for mode in modes:
        y = savgol_filter(signal_noisy, window_length, polyorder, mode=mode)
        first_vals[mode].append(y[0])
print(f"\ntrue value at t=0: {signal_clean[0]:.5f}")
for mode in modes:
    vals = np.array(first_vals[mode])
    print(f"{mode:>10}: bias={np.mean(vals) - signal_clean[0]:+.5f}  std={np.std(vals):.5f}")

import numpy as np
from scipy.signal import savgol_filter
import matplotlib.pyplot as plt

np.random.seed(42)
t = np.linspace(0, 1, 500)
signal_clean = np.sin(2 * np.pi * 3 * t) + 0.5 * np.exp(-((t - 0.3) / 0.02) ** 2)
signal_noisy = signal_clean + 0.2 * np.random.randn(len(t))

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 窓幅の影響（多項式次数固定: p=3）
for w in [5, 11, 21, 41]:
    y = savgol_filter(signal_noisy, window_length=w, polyorder=3)
    axes[0].plot(t, y, label=f'w={w}')
axes[0].plot(t, signal_clean, 'k--', alpha=0.5, label='Original')
axes[0].set_title('Effect of Window Length (polyorder=3)')
axes[0].legend()
axes[0].grid(True)

# 多項式次数の影響（窓幅固定: w=21）
for p in [0, 1, 2, 3, 5]:
    y = savgol_filter(signal_noisy, window_length=21, polyorder=p)
    axes[1].plot(t, y, label=f'p={p}')
axes[1].plot(t, signal_clean, 'k--', alpha=0.5, label='Original')
axes[1].set_title('Effect of Polynomial Order (window=21)')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.show()

import numpy as np
from scipy.signal import savgol_filter

t = np.linspace(0, 1, 500)
signal_clean = np.sin(2 * np.pi * 3 * t) + 0.5 * np.exp(-((t - 0.3) / 0.02) ** 2)

windows = [5, 9, 11, 15, 21, 31, 41, 61]
orders = [0, 2, 3, 5]
sigma = 0.2
n_trials = 80
m_exclude = 30

results = {p: [] for p in orders}
for p in orders:
    for w in windows:
        if p >= w:
            results[p].append(None)
            continue
        rng = np.random.default_rng(hash((w, p)) % (2**31))
        errs = []
        for _ in range(n_trials):
            noise = sigma * rng.standard_normal(len(t))
            y = signal_clean + noise
            ys = savgol_filter(y, w, p)
            e = ys[m_exclude:-m_exclude] - signal_clean[m_exclude:-m_exclude]
            errs.append(np.sqrt(np.mean(e**2)))
        results[p].append(np.mean(errs))

for p in orders:
    row = ", ".join(f"w={w}:{v:.4f}" for w, v in zip(windows, results[p]) if v is not None)
    print(f"p={p}: {row}")

import numpy as np
from scipy.signal import savgol_filter

t = np.linspace(0, 1, 500)
peak_center, peak_width, peak_height = 0.5, 0.02, 1.0
signal_clean = peak_height * np.exp(-((t - peak_center) / peak_width) ** 2)
true_peak_idx = np.argmax(signal_clean)
true_peak_val = signal_clean[true_peak_idx]

def moving_average(x, w):
    kernel = np.ones(w) / w
    return np.convolve(x, kernel, mode='same')

def ema(x, alpha):
    y = np.zeros_like(x)
    y[0] = x[0]
    for i in range(1, len(x)):
        y[i] = alpha * x[i] + (1 - alpha) * y[i - 1]
    return y

n_trials = 200
sigma = 0.1
w = 21
alpha = 2 / (w + 1)  # windowと同等のスパンに合わせたEMA係数

methods = {
    "Moving Average": lambda y: moving_average(y, w),
    "EMA": lambda y: ema(y, alpha),
    "SG (p=2)": lambda y: savgol_filter(y, w, 2),
    "SG (p=3)": lambda y: savgol_filter(y, w, 3),
    "SG (p=5)": lambda y: savgol_filter(y, w, 5),
}

peak_heights = {name: [] for name in methods}
rng = np.random.default_rng(5)
for _ in range(n_trials):
    noise = sigma * rng.standard_normal(len(t))
    y = signal_clean + noise
    search = slice(max(0, true_peak_idx - 30), true_peak_idx + 31)
    for name, fn in methods.items():
        yf = fn(y)
        peak_heights[name].append(np.max(yf[search]))

print(f"true peak height = {true_peak_val:.3f}")
for name, vals in peak_heights.items():
    vals = np.array(vals)
    mean = np.mean(vals)
    bias_pct = (mean / true_peak_val - 1) * 100
    print(f"{name:>16}: mean peak={mean:.4f} (bias {bias_pct:+.1f}%), std={np.std(vals):.4f}")

h = 1e-6
f = lambda tt: peak_height * np.exp(-((tt - peak_center) / peak_width) ** 2)
f2 = (f(peak_center + h) - 2 * f(peak_center) + f(peak_center - h)) / h**2
dt = t[1] - t[0]
f2_sample = f2 * dt**2  # サンプル単位の2階差分に変換
m = w // 2
predicted_bias_pct = (f2_sample / 2) * (m * (m + 1) / 3) / true_peak_val * 100
print(f"式(8)による予測バイアス: {predicted_bias_pct:.1f}%  (実測: -26.9%)")

import numpy as np
from scipy.signal import savgol_coeffs

rng = np.random.default_rng(0)
sigma = 1.0
N = 200000
noise = sigma * rng.standard_normal(N)

cases = [(11, 2, 0), (21, 2, 0), (41, 2, 0),
         (11, 3, 1), (21, 3, 1), (41, 3, 1),
         (11, 4, 2), (21, 4, 2), (41, 4, 2)]

for w, p, d in cases:
    c = savgol_coeffs(w, p, deriv=d, use='dot')
    theory = sigma**2 * np.sum(c**2)
    filtered = np.convolve(noise, c[::-1], mode='valid')  # dot規約でのフィルタ適用
    emp = np.var(filtered)
    print(f"w={w:3d} p={p} d={d}: theory={theory:.6f}  empirical={emp:.6f}")

import numpy as np
from scipy.signal import savgol_filter

t = np.linspace(0, 1, 500)
dt = t[1] - t[0]
signal_clean = np.sin(2 * np.pi * 3 * t) + 0.5 * np.exp(-((t - 0.3) / 0.02) ** 2)
dy_true = 2 * np.pi * 3 * np.cos(2 * np.pi * 3 * t) + \
          0.5 * (-2 * (t - 0.3) / 0.02**2) * np.exp(-((t - 0.3) / 0.02) ** 2)

windows_b = [5, 7, 9, 11, 15, 21, 31, 41, 61, 81]
noise_levels = [0.02, 0.1, 0.3]
polyorder = 3
m_exclude = 30

for sigma in noise_levels:
    rmses = []
    for w in windows_b:
        rng = np.random.default_rng(int(sigma * 1000) + w)
        errs = []
        for _ in range(60):
            noise = sigma * rng.standard_normal(len(t))
            y = signal_clean + noise
            dy = savgol_filter(y, w, polyorder, deriv=1, delta=dt)
            e = dy[m_exclude:-m_exclude] - dy_true[m_exclude:-m_exclude]
            errs.append(np.sqrt(np.mean(e**2)))
        rmses.append(np.mean(errs))
    best_idx = int(np.argmin(rmses))
    print(f"sigma={sigma}: best window={windows_b[best_idx]}, RMSE={rmses[best_idx]:.4f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import savgol_coeffs, freqz

# フィルタ係数の取得
for w, p in [(11, 2), (21, 3), (41, 3)]:
    coeffs = savgol_coeffs(w, p)
    freq, response = freqz(coeffs, worN=2048)
    plt.plot(freq / np.pi, 20 * np.log10(np.abs(response)),
             label=f'w={w}, p={p}')

plt.xlabel('Normalized Frequency (×π rad/sample)')
plt.ylabel('Gain (dB)')
plt.title('Savitzky-Golay Filter Frequency Response')
plt.legend()
plt.grid(True)
plt.ylim(-60, 5)
plt.tight_layout()
plt.show()