@Seine_A_Shintani(新谷正嶺)Team Coffeein

中部大学新谷研究室 / Shintani Lab, Chubu University

フィルタをかけると波形はきれいになる。でも真実に近づいたとは限らない：Pythonで見る平滑化の落とし穴

Posted at 2026-07-01

はじめに

時系列データを扱っていると、こんな場面によく出会います。

波形がノイズでギザギザしている。
フィルタをかけると、なめらかできれいになった。
では、この波形は「真の波形」に近づいたと言ってよいのだろうか？

直感的には、ノイズが減って見えると「良くなった」と感じます。

しかし、ここに落とし穴があります。

フィルタは、波形をきれいにする処理ではありますが、真の波形を復元する魔法ではありません。

特に、ピーク値、ピーク時刻、立ち上がり速度、面積などを定量したい場合、フィルタによって結果そのものが変わることがあります。

この記事では、Pythonでダミーデータを作り、

ノイズを含む観測波形
移動平均
Savitzky-Golayフィルタ
ローパスフィルタ
強すぎるローパスフィルタ

を比較します。

そして、

見た目はきれいになったのに、解析値は真値から遠ざかることがある

ということを、図と数値で確認します。

この記事で伝えたいこと

この記事の結論は次の通りです。

フィルタは、波形に含まれる成分を選別する処理である
なめらかな波形が、必ずしも正しい波形とは限らない
フィルタの強さによって、ピーク値、ピーク時刻、傾きが変わる
特に微分や速度を計算する場合、フィルタの影響は大きい
フィルタ条件は、解析目的に合わせて決める必要がある

対象読者

この記事は、次のような方を想定しています。

Pythonで時系列データを扱い始めた人
実験データ、センサーデータ、波形データを解析したい人
「フィルタ」「平滑化」「ローパス」という言葉を聞いたことはあるが、直感的に理解したい人
データをきれいにする前処理の落とし穴を知りたい人

専門的な信号処理の理論には深く入りません。
まずは、図で理解することを重視します。

用語の整理

波形

ここでは、時間とともに値が変化するデータを波形と呼びます。

たとえば、

センサー値
電圧
明るさ
位置
速度
生体信号
実験装置の出力

などは、時間方向に並べると波形として扱えます。

真の信号と観測信号

実際の測定では、本当に知りたい信号にノイズが混ざります。

この記事では、次のように考えます。

$$
x_{\mathrm{obs}}(t)=
x_{\mathrm{true}}(t)
+
\eta(t)
$$

ここで、

$x_{\mathrm{true}}(t)$：本当は知りたい真の信号
$x_{\mathrm{obs}}(t)$：実際に観測された信号
$\eta(t)$：ノイズやアーティファクト

です。

実データでは、通常 $x_{\mathrm{true}}(t)$ は分かりません。

しかし、今回は教育用のダミーデータを作るので、あえて「真の信号」をこちらで定義します。
これにより、フィルタ後の波形が真値に近いのか遠いのかを比較できます。

ノイズ

ノイズとは、知りたい信号に混ざる不要な変動のことです。

ただし、実際には何がノイズで何が信号かは、解析目的によって変わります。

たとえば、高周波のギザギザが不要なノイズに見える場合もあります。
一方で、速い立ち上がりや短いイベントを見たい場合、その高周波成分は重要な信号かもしれません。

この点が、フィルタを難しくしている理由の一つです。

フィルタ

フィルタとは、波形に含まれる成分を選別する処理です。

簡単に言えば、

残したい成分を残し、弱めたい成分を弱める処理

です。

たとえば、ローパスフィルタは低い周波数の成分を残し、高い周波数の成分を弱めます。

ただし、高周波成分を弱めるということは、ノイズだけでなく、波形の急な変化も弱める可能性がある、ということです。

平滑化

平滑化は、波形をなめらかにする処理です。

代表例が移動平均です。

中心化した移動平均は、次のように書けます。

\tilde{x}_i=
\frac{1}{2m+1}
\sum_{k=-m}^{m}
x_{i+k}

これは、ある点の前後の値を平均して、その点の値にする処理です。

分かりやすい一方で、急なピークや立ち上がりは丸められやすくなります。

微分とノイズ

時系列データから速度のような量を求めたい場合、差分を使って微分を近似します。

$$
\frac{dx}{dt}
\approx
\frac{x_{i+1}-x_i}{\Delta t}
$$

ここで重要なのは、隣り合う点の差を取るということです。

ノイズがあると、隣り合う点の差が大きくなりやすいため、微分するとノイズが強調されます。

つまり、

位置や強度の波形では小さく見えるノイズが、速度や傾きでは大きく暴れる

ことがあります。

今回作るダミーデータ

今回は、急に立ち上がってゆっくり戻る波形を作ります。

これは、特定の実験系や生体信号を厳密に再現するものではありません。
あくまで、フィルタの性質を学ぶための教育用ダミーデータです。

観測信号には、次の3種類のノイズ・アーティファクトを加えます。

ランダムな白色ノイズ
高周波の周期的な揺れ
ゆっくりしたベースラインドリフト

そのうえで、複数のフィルタを比較します。

Google Colabで実行するコード

以下のコードは、Google Colabでそのまま実行できます。

# ============================================
# Filtering demo:
# A smooth-looking waveform is not always closer to truth
# ============================================

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from scipy.signal import savgol_filter, butter, filtfilt
from scipy.integrate import trapezoid

# -----------------------------
# 1. Settings
# -----------------------------
SEED = 202601
rng = np.random.default_rng(SEED)

fs = 200       # sampling frequency [Hz]
dt = 1 / fs
T = 12         # duration [s]
t = np.arange(0, T, dt)

os.makedirs("fig", exist_ok=True)

# -----------------------------
# 2. Create dummy "true" waveform
# -----------------------------
def transient(t, t0, amp=1.0, tau_rise=0.04, tau_decay=0.35):
    """
    Asymmetric transient:
    - rapid rise
    - slow decay

    This is a toy waveform, not a specific biological/physical model.
    """
    x = t - t0
    y = np.zeros_like(t)
    mask = x >= 0

    raw = (1 - np.exp(-x[mask] / tau_rise)) * np.exp(-x[mask] / tau_decay)

    # Normalize the peak to 1.0
    if raw.size > 0:
        peak = np.max(raw)
        if peak > 0:
            raw = raw / peak

    y[mask] = amp * raw
    return y

event_times = np.array([1.0, 3.1, 5.4, 7.9, 10.2])
event_amps = np.array([1.0, 0.85, 1.15, 0.95, 1.05])

true_signal = np.zeros_like(t)
for t0, amp in zip(event_times, event_amps):
    true_signal += transient(t, t0, amp=amp)

# -----------------------------
# 3. Add measurement noise/artifacts
# -----------------------------
white_noise = 0.07 * rng.normal(size=t.size)
high_freq_artifact = 0.06 * np.sin(2 * np.pi * 25 * t)
baseline_drift = 0.04 * np.sin(2 * np.pi * 0.20 * t + 0.5)

observed_signal = true_signal + white_noise + high_freq_artifact + baseline_drift

# -----------------------------
# 4. Define filters
# -----------------------------
def moving_average_centered(x, window):
    """
    Centered moving average with edge padding.
    window must be an odd integer.
    """
    if window % 2 == 0:
        raise ValueError("window must be odd")

    kernel = np.ones(window) / window
    pad = window // 2
    x_pad = np.pad(x, (pad, pad), mode="edge")

    return np.convolve(x_pad, kernel, mode="valid")

def butter_lowpass_filtfilt(x, cutoff_hz, fs, order=4):
    """
    Zero-phase Butterworth low-pass filter using filtfilt.
    """
    nyquist = fs / 2
    b, a = butter(order, cutoff_hz / nyquist, btype="low")

    return filtfilt(b, a, x)

moving_avg_105ms = moving_average_centered(observed_signal, window=21)
savgol_105ms = savgol_filter(observed_signal, window_length=21, polyorder=3)
lowpass_8hz = butter_lowpass_filtfilt(observed_signal, cutoff_hz=8, fs=fs, order=4)
lowpass_2hz = butter_lowpass_filtfilt(observed_signal, cutoff_hz=2, fs=fs, order=4)

signals = {
    "true": true_signal,
    "observed": observed_signal,
    "moving_avg_105ms": moving_avg_105ms,
    "savgol_105ms": savgol_105ms,
    "lowpass_8Hz": lowpass_8hz,
    "lowpass_2Hz": lowpass_2hz,
}

# -----------------------------
# 5. Metrics around one target event
# -----------------------------
target_t0 = event_times[2]

analysis_window = (t >= target_t0 - 0.2) & (t <= target_t0 + 1.2)
baseline_window = (t >= target_t0 - 0.2) & (t < target_t0)

def compute_metrics(y):
    tw = t[analysis_window]
    yw = y[analysis_window]

    # Simple baseline correction
    baseline = np.mean(y[baseline_window])
    ywc = yw - baseline

    idx_peak = np.argmax(ywc)
    peak_amp = ywc[idx_peak]
    peak_time = tw[idx_peak]

    dy_dt = np.gradient(y - baseline, dt)
    max_rising_slope = np.max(dy_dt[analysis_window])

    auc = trapezoid(ywc, tw)

    true_baseline = np.mean(true_signal[baseline_window])
    true_ywc = true_signal[analysis_window] - true_baseline
    rmse = np.sqrt(np.mean((ywc - true_ywc) ** 2))

    return peak_amp, peak_time, max_rising_slope, auc, rmse

rows = []
for name, y in signals.items():
    rows.append((name,) + compute_metrics(y))

metrics = pd.DataFrame(
    rows,
    columns=[
        "signal",
        "peak_amp",
        "peak_time_s",
        "max_rising_slope",
        "auc",
        "rmse_vs_true",
    ],
)

true_peak_amp = metrics.loc[metrics["signal"] == "true", "peak_amp"].iloc[0]
true_peak_time = metrics.loc[metrics["signal"] == "true", "peak_time_s"].iloc[0]
true_max_slope = metrics.loc[metrics["signal"] == "true", "max_rising_slope"].iloc[0]

metrics["peak_amp_error_%"] = (metrics["peak_amp"] / true_peak_amp - 1) * 100
metrics["peak_time_delay_ms"] = (metrics["peak_time_s"] - true_peak_time) * 1000
metrics["max_slope_error_%"] = (metrics["max_rising_slope"] / true_max_slope - 1) * 100

metrics_rounded = metrics.round({
    "peak_amp": 3,
    "peak_time_s": 3,
    "max_rising_slope": 2,
    "auc": 3,
    "rmse_vs_true": 3,
    "peak_amp_error_%": 1,
    "peak_time_delay_ms": 1,
    "max_slope_error_%": 1,
})

metrics_rounded.to_csv("fig/metrics_table.csv", index=False)

try:
    display(metrics_rounded)
except NameError:
    print(metrics_rounded.to_string(index=False))

# -----------------------------
# 6. Figure 1: true and observed signals
# -----------------------------
plt.figure(figsize=(11, 4))
plt.plot(t, true_signal, label="true signal")
plt.plot(t, observed_signal, label="observed signal", alpha=0.7)
plt.title("Fig. 1: Dummy waveform with noise and artifacts")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude [a.u.]")
plt.legend()
plt.tight_layout()
plt.savefig("fig/fig1_true_observed.png", dpi=160)
plt.show()

# -----------------------------
# 7. Figure 2: filtered signals
# -----------------------------
plt.figure(figsize=(11, 4))
plt.plot(t, true_signal, label="true signal", linewidth=2)
plt.plot(t, observed_signal, label="observed", alpha=0.35)
plt.plot(t, moving_avg_105ms, label="moving average 105 ms")
plt.plot(t, savgol_105ms, label="Savitzky-Golay 105 ms")
plt.plot(t, lowpass_8hz, label="low-pass 8 Hz")
plt.plot(t, lowpass_2hz, label="low-pass 2 Hz")
plt.title("Fig. 2: Filtering makes the waveform look smoother")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude [a.u.]")
plt.legend(ncol=2)
plt.tight_layout()
plt.savefig("fig/fig2_filter_overview.png", dpi=160)
plt.show()

# -----------------------------
# 8. Figure 3: zoom around one event
# -----------------------------
zoom = (t >= target_t0 - 0.2) & (t <= target_t0 + 1.2)

plt.figure(figsize=(10, 5))
for name, y in signals.items():
    if name == "observed":
        plt.plot(t[zoom], y[zoom], label=name, alpha=0.35)
    else:
        plt.plot(t[zoom], y[zoom], label=name)

# Peak markers
for name, y in signals.items():
    peak_amp, peak_time, *_ = compute_metrics(y)
    baseline = np.mean(y[baseline_window])
    peak_y = peak_amp + baseline
    plt.scatter(peak_time, peak_y, s=35)

plt.title("Fig. 3: Zoom view around one transient")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude [a.u.]")
plt.legend(ncol=2)
plt.tight_layout()
plt.savefig("fig/fig3_zoom_peak.png", dpi=160)
plt.show()

# -----------------------------
# 9. Figure 4: derivative comparison
# -----------------------------
plt.figure(figsize=(10, 5))
for name in [
    "true",
    "observed",
    "moving_avg_105ms",
    "savgol_105ms",
    "lowpass_8Hz",
    "lowpass_2Hz",
]:
    y = signals[name]
    dy_dt = np.gradient(y, dt)

    if name == "observed":
        plt.plot(t[zoom], dy_dt[zoom], label=name, alpha=0.25)
    else:
        plt.plot(t[zoom], dy_dt[zoom], label=name)

plt.title("Fig. 4: Differentiation amplifies noise and reveals smoothing effects")
plt.xlabel("Time [s]")
plt.ylabel("Derivative [a.u./s]")
plt.legend(ncol=2)
plt.tight_layout()
plt.savefig("fig/fig4_derivative.png", dpi=160)
plt.show()

# -----------------------------
# 10. Figure 5: error summary
# -----------------------------
plot_df = metrics[metrics["signal"] != "true"].copy()

fig, axes = plt.subplots(1, 3, figsize=(13, 4))

axes[0].bar(plot_df["signal"], plot_df["peak_amp_error_%"])
axes[0].axhline(0, linestyle="--", linewidth=1)
axes[0].set_title("Peak amplitude error")
axes[0].set_ylabel("Error [%]")
axes[0].tick_params(axis="x", rotation=45)

axes[1].bar(plot_df["signal"], plot_df["peak_time_delay_ms"])
axes[1].axhline(0, linestyle="--", linewidth=1)
axes[1].set_title("Peak time delay")
axes[1].set_ylabel("Delay [ms]")
axes[1].tick_params(axis="x", rotation=45)

axes[2].bar(plot_df["signal"], plot_df["max_slope_error_%"])
axes[2].axhline(0, linestyle="--", linewidth=1)
axes[2].set_title("Max rising slope error")
axes[2].set_ylabel("Error [%]")
axes[2].tick_params(axis="x", rotation=45)

plt.tight_layout()
plt.savefig("fig/fig5_metric_errors.png", dpi=160)
plt.show()

print("Saved figures to: fig/")
print("Saved metrics table to: fig/metrics_table.csv")

図1：真の波形と観測波形

図1では、青線が今回こちらで定義した「真の波形」、オレンジ色の線がノイズやアーティファクトを加えた「観測波形」です。

観測波形には、細かいギザギザが含まれています。

このような波形を見ると、自然に「フィルタをかけてきれいにしたい」と感じます。

ただし、ここで大切なのは、観測波形が汚いからといって、フィルタ後の波形が必ず真の波形に近づくとは限らないという点です。

図2：フィルタをかけると見た目はなめらかになる

図2では、観測波形に複数のフィルタをかけた結果を比較しています。

今回比較したのは、次の4種類です。

名前	内容
`moving_avg_105ms`	105 ms幅の移動平均
`savgol_105ms`	105 ms幅のSavitzky-Golayフィルタ
`lowpass_8Hz`	8 Hzのローパスフィルタ
`lowpass_2Hz`	2 Hzの強いローパスフィルタ

どの方法でも、観測波形よりはなめらかに見えます。

特に lowpass_2Hz はかなりなめらかです。

しかし、なめらかに見えることと、真値に近いことは別問題です。

図3：ピーク付近を拡大して見る

図3では、1つのイベント周辺を拡大しています。

ここで注目したいのは、ピーク付近です。

lowpass_2Hz は非常になめらかな波形になっています。
しかし、真の波形と比べると、ピークが低くなり、ピーク時刻も遅れています。

つまり、

強いフィルタによって、ノイズだけでなく、本来の速い変化も削られている

ということです。

これは非常に重要です。

波形の見た目だけで判断すると、lowpass_2Hz は「きれい」に見えるかもしれません。
しかし、ピーク値やピーク時刻を解析したい場合には、むしろ歪んだ波形になっています。

図4：微分すると違いがさらに目立つ

図4では、波形の微分を比較しています。

微分は、ざっくり言えば「どれくらい速く変化しているか」を表します。

差分で近似すると、次のようになります。

$$
\frac{dx}{dt}
\approx
\frac{x_{i+1}-x_i}{\Delta t}
$$

観測波形をそのまま微分すると、ノイズが大きく暴れます。

一方、フィルタをかけると微分は落ち着きます。
しかし、ここでも注意が必要です。

強いフィルタをかけると、立ち上がりの速さまで小さく見積もられます。

つまり、

微分前のフィルタは必要なことが多いが、フィルタによって速度や傾きが変わる

ということです。

「速度」「立ち上がり速度」「最大変化率」などを解析する場合、フィルタ条件は結果に強く影響します。

数値で比較する

今回のコードでは、3番目のイベントについて次の量を計算しました。

指標	意味
`peak_amp`	ベースライン補正後のピーク振幅
`peak_time_s`	ピーク時刻
`max_rising_slope`	最大立ち上がり速度
`auc`	波形の面積
`rmse_vs_true`	真の波形からのずれ
`peak_amp_error_%`	ピーク振幅の誤差
`peak_time_delay_ms`	ピーク時刻のずれ
`max_slope_error_%`	最大立ち上がり速度の誤差

RMSEは、次のような指標です。

RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(\tilde{x}_i-x_{\mathrm{true},i}\right)^2}

ここで、

$\tilde{x}_i$：フィルタ後の値
$x_{\mathrm{true},i}$：真の値
$N$：比較した点の数

です。

RMSEが小さいほど、今回のダミーデータにおける真の波形に近いことを意味します。

実行結果の例

乱数seedを固定しているので、基本的には同じ結果になるはずです。
私の環境では、次のような結果になりました。

signal	peak_amp	peak_time_s	max_rising_slope	auc	rmse_vs_true	peak_amp_error_%	peak_time_delay_ms	max_slope_error_%
true	1.149	5.490	35.73	0.501	0.000	0.0	0.0	0.0
observed	1.339	5.490	50.07	0.501	0.076	16.6	0.0	40.1
moving_avg_105ms	1.092	5.495	12.90	0.463	0.055	-5.0	5.0	-63.9
savgol_105ms	1.162	5.480	22.12	0.502	0.025	1.1	-10.0	-38.1
lowpass_8Hz	1.186	5.485	16.38	0.502	0.029	3.2	-5.0	-54.2
lowpass_2Hz	0.835	5.580	4.29	0.318	0.180	-27.4	90.0	-88.0

図5：誤差を可視化する

図5では、ピーク振幅、ピーク時刻、最大立ち上がり速度の誤差をまとめています。

ここで特に重要なのは、lowpass_2Hz です。

lowpass_2Hz は、見た目にはかなりなめらかです。
しかし、数値を見ると、

ピーク振幅は約27%小さい
ピーク時刻は約90 ms遅れる
最大立ち上がり速度は約88%小さい
RMSEは観測波形よりも大きい

という結果になっています。

つまり、今回の例では、

最もなめらかに見える波形が、真の波形に最も近いとは限らない

ことが確認できます。

なぜこうなるのか

ローパスフィルタは、高い周波数成分を弱めます。

これはノイズ除去には有効です。

しかし、波形の急な立ち上がりや鋭いピークも、高い周波数成分を含みます。

そのため、フィルタを強くしすぎると、

ピークが低くなる
立ち上がりが遅くなる
ピーク時刻がずれる
傾きが小さくなる
面積が変わる

といった変化が起こります。

つまり、フィルタは単なる「ノイズ除去」ではなく、波形そのものを変換する処理です。

フィルタ条件は解析目的で変わる

フィルタの良し悪しは、目的によって変わります。

たとえば、目的が「大まかなトレンドを見たい」だけなら、強めのフィルタでもよいかもしれません。

一方で、次のような量を測りたい場合は注意が必要です。

ピーク値
ピーク時刻
立ち上がり速度
減衰速度
周期
イベントの開始時刻
微分値
面積

これらは、フィルタ条件によって変わりやすい指標です。

特に、微分値はフィルタの影響を強く受けます。

実データ解析での注意点

実データでは、今回のように真の波形は分かりません。

だからこそ、次のような確認が重要です。

1. 生データを必ず残す

フィルタ後のデータだけを保存するのは危険です。

生データ、フィルタ後データ、フィルタ条件をセットで保存しておくべきです。

2. フィルタ条件を明記する

たとえば、論文、レポート、解析メモでは、

フィルタの種類
窓幅
カットオフ周波数
フィルタ次数
前後方向フィルタかどうか
解析に使ったライブラリ

を明記したほうがよいです。

「平滑化した」とだけ書いても、再現できません。

3. 可視化用と解析用を分けて考える

見やすい図を作るためのフィルタと、定量解析に使うフィルタは同じでよいとは限りません。

可視化用に強めに平滑化した波形を、そのままピーク解析に使うと、値が歪むことがあります。

4. 指標ごとに影響を確認する

フィルタの影響は、指標によって異なります。

ピーク振幅にはあまり影響しなくても、最大立ち上がり速度には大きく影響することがあります。

今回の例でも、savgol_105ms はピーク振幅を比較的よく保っています。
しかし、最大立ち上がり速度は真値よりかなり小さくなっています。

5. 合成データで検証する

今回のように、真の波形が分かっている合成データを使うと、フィルタの性質を調べやすくなります。

実データにいきなり適用する前に、

どの程度のノイズを消せるか
ピーク値がどれくらい変わるか
タイミングがずれるか
微分値がどう変わるか

を確認できます。

おまけ：パラメータを変えて遊んでみる

以下の値を変えると、結果が大きく変わります。

ノイズを強くする

white_noise = 0.15 * rng.normal(size=t.size)

ノイズが強くなると、フィルタの必要性は高くなります。
一方で、強いフィルタを使いたくなるため、波形の歪みにも注意が必要になります。

移動平均の窓幅を変える

moving_avg_105ms = moving_average_centered(observed_signal, window=41)

窓幅を大きくすると、よりなめらかになります。
しかし、ピークや立ち上がりは丸まりやすくなります。

ローパスフィルタのカットオフ周波数を変える

lowpass_4hz = butter_lowpass_filtfilt(observed_signal, cutoff_hz=4, fs=fs, order=4)

カットオフ周波数を低くすると、より強い平滑化になります。

ただし、急な変化も失われます。

まとめ

この記事では、ダミーデータを使って、フィルタが波形解析に与える影響を見ました。

ポイントは次の通りです。

フィルタをかけると波形はなめらかになる
しかし、なめらかさは正しさを保証しない
強いフィルタは、ノイズだけでなく本来の速い変化も削る
ピーク値、ピーク時刻、最大立ち上がり速度はフィルタで変わる
微分値は特にフィルタ条件の影響を受けやすい
フィルタは、目的に合わせて選び、条件を明記する必要がある

最後に、この記事の一番大事なメッセージをもう一度書きます。

フィルタは「真実に近づける魔法」ではなく、「仮定に基づいて波形を変換する処理」である。

波形がきれいになったときほど、解析値がどう変わったかを確認することが大切です。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up