信号処理におけるフィルタリング手法の基礎

フィルタリングとは

フィルタリングとは、観測された時系列データに含まれるノイズ成分を取り除き、本来の信号成分を抽出するための技術です。

状態空間モデル

フィルタリングでは、直接観測できない内部状態 \(x_t\) を、観測値 \(z_t\) から推定します。このとき、システムの振る舞いを以下の2つのモデルで表現します。

プロセスモデル（システムモデル）: 状態 \(x_t\) が時間と共にどう変化するかを記述します。

\[ x_t = f(x_{t-1}, u_{t-1}) + q_{t-1} \tag{1} \]

観測モデル: 状態 \(x_t\) から観測値 \(z_t\) がどのように生成されるかを記述します。

\[ z_t = h(x_t) + r_t \tag{2} \]

ここで、\(u_t\) は制御入力、\(q_t\) と \(r_t\) はそれぞれプロセスノイズと観測ノイズであり、一般的に平均0の正規分布（ガウス分布）に従うと仮定されます。

\[ q \sim \mathcal{N}(0, Q), \quad r \sim \mathcal{N}(0, R) \]

\(Q\) と \(R\) はノイズの共分散行列です。\(f\) と \(h\) がともに線形写像（\(f(x)=Ax\) , \(h(x)=Hx\) ）かつノイズがガウス分布のときに限り、後述のカルマンフィルタが厳密に最適な解を与えます。それ以外の非線形・非ガウスの場合には、EKF・UKF・粒子フィルタといった近似手法が必要になります。

カルマンフィルタ (Kalman Filter, KF)

概要

KFは、線形のシステムかつ、ノイズがガウス分布に従うという仮定のもとで、状態の平均と共分散を正確に推定する最適なフィルタです。「予測」と「更新（フィルタリング）」の2つのステップを繰り返すことで、状態を逐次的に推定します。

アルゴリズム

状態 \(x\) の確率分布を、平均 \(\mu\) と共分散 \(\Sigma\) で表現します。

1. 予測ステップ: 1ステップ前の状態推定値から、現在の状態を予測します。

事前状態推定値: 1ステップ前の状態 \(\mu_{t-1}\) から現在の状態 \(\hat{\mu}_t\) を予測。

\[ \hat{\mu}_t = A\mu_{t-1} + Bu_{t-1} \tag{3} \]

事前誤差共分散行列: 予測の不確かさ \(\hat{\Sigma}_t\) を計算。

\[ \hat{\Sigma}_t = A\Sigma_{t-1}A^T + Q \tag{4} \]

2. 更新ステップ: 予測結果を観測値 \(z_t\) で補正し、より確からしい現在の状態を推定します。

カルマンゲイン: 予測と観測のどちらをどの程度重視するかを決定する係数。

\[ K_t = \hat{\Sigma}_t H^T (H\hat{\Sigma}_t H^T + R)^{-1} \tag{5} \]

事後状態推定値: 予測値 \(\hat{\mu}_t\) を観測値 \(z_t\) で補正。

\[ \mu_t = \hat{\mu}_t + K_t(z_t - H\hat{\mu}_t) \tag{6} \]

事後誤差共分散行列: 更新後の不確かさ \(\Sigma_t\) を計算。

\[ \Sigma_t = (I - K_tH)\hat{\Sigma}_t \tag{7} \]

導出：なぜカルマンゲインの式が最適なのか

式(5)のカルマンゲインは天下り的に見えますが、実は「予測値と観測値という2つの不偏推定量を、事後分散が最小になるように線形結合する」という最適化問題を解いた結果です。詳細な多次元導出と実装はカルマンフィルタの理論とPython実装に譲り、ここではスカラー（1次元、\(H=1\) ）の場合で本質を示します。

状態の予測値 \(\hat{\mu}_t\) （誤差分散 \(\hat{\sigma}_t^2\) ）と観測値 \(z_t\) （観測ノイズ分散 \(R\) ）は、真の状態 \(x_t\) に対して独立な不偏推定量です。この2つを重み \(K\) で線形結合した推定量

\[ \mu_t = (1-K)\hat{\mu}_t + K z_t \]

を考えます。\(\mathbb{E}[\hat{\mu}_t]=\mathbb{E}[z_t]=x_t\) なので、\(K\) をどう選んでも \(\mu_t\) は不偏です。両者が独立であることから、分散は

\[ \mathrm{Var}(\mu_t) = (1-K)^2 \hat{\sigma}_t^2 + K^2 R \tag{8} \]

となります。この分散を最小にする \(K\) を、\(K\) について微分して0とおくことで求めます。

\[ \frac{d}{dK}\mathrm{Var}(\mu_t) = -2(1-K)\hat{\sigma}_t^2 + 2KR = 0 \]

これを \(K\) について解くと、

\[ K = \frac{\hat{\sigma}_t^2}{\hat{\sigma}_t^2 + R} \tag{9} \]

が得られます。これは式(5)で \(H=1\) とした場合と完全に一致します。つまりカルマンゲインは、予測誤差分散と観測ノイズ分散の比で決まる「精度（分散の逆数）による重み付け平均」であり、予測が不確か（\(\hat{\sigma}_t^2\) が大きい）なら観測を重視し、観測が不確か（\(R\) が大きい）なら予測を重視するように、分散最小の観点から自動的に調整されます。さらに式(9)を式(8)に代入すると、

\[ \mathrm{Var}(\mu_t)_{\min} = (1-K)\hat{\sigma}_t^2 = \frac{\hat{\sigma}_t^2 R}{\hat{\sigma}_t^2 + R} \]

となり、これは式(7)を \(H=1\) で評価した結果と一致します。多次元・行列の場合も同様に、事後共分散 \(\Sigma_t\) のトレース（各成分の分散の和）を \(K_t\) について最小化することで式(5)が導かれます（ガウス・マルコフの定理に基づく最良線形不偏推定量、BLUE）。ノイズがガウス分布に従う場合は、事後分布自体がガウス分布になるため、この線形推定量は線形推定量の中で最良であるだけでなく、非線形推定量も含めた全推定量の中で分散最小（最尤推定量かつ事後平均）になることが知られています。

数値実験：1次元ランダムウォークによる推定

導出したカルマンフィルタが実際にノイズを低減できることを、最も単純な1次元モデルで確認します。真の状態がランダムウォーク \(x_t = x_{t-1} + q_t\) （\(q_t \sim \mathcal{N}(0, Q)\) 、\(Q=0.25\) ）に従い、観測が \(z_t = x_t + r_t\) （\(r_t \sim \mathcal{N}(0, R)\) 、\(R=4.0\) ）で得られる状況を考えます（\(A=H=1\) のスカラー版カルマンフィルタ）。

import numpy as np

np.random.seed(0)

T = 50
Q = 0.25  # プロセスノイズ分散
R = 4.0  # 観測ノイズ分散

x_true = np.zeros(T + 1)
z = np.zeros(T + 1)
for t in range(1, T + 1):
    x_true[t] = x_true[t - 1] + np.random.normal(0, np.sqrt(Q))
z[1:] = x_true[1:] + np.random.normal(0, np.sqrt(R), size=T)

# --- スカラーカルマンフィルタ（A=H=1） ---
mu, sigma2 = 0.0, 1.0  # 初期状態推定値・初期分散
mu_hist, sigma2_hist, K_hist = [mu], [sigma2], []

for t in range(1, T + 1):
    mu_pred = mu  # 予測ステップ（式3、A=1）
    sigma2_pred = sigma2 + Q  # 予測ステップ（式4）

    K = sigma2_pred / (sigma2_pred + R)  # カルマンゲイン（式9）
    mu = mu_pred + K * (z[t] - mu_pred)  # 更新ステップ（式6）
    sigma2 = (1 - K) * sigma2_pred  # 更新ステップ（式7）

    mu_hist.append(mu)
    sigma2_hist.append(sigma2)
    K_hist.append(K)

mu_hist, sigma2_hist = np.array(mu_hist), np.array(sigma2_hist)
rmse_obs = np.sqrt(np.mean((z[1:] - x_true[1:]) ** 2))
rmse_kf = np.sqrt(np.mean((mu_hist[1:] - x_true[1:]) ** 2))

print(f"K_1 = {K_hist[0]:.4f}, K_50 = {K_hist[-1]:.4f}")
print(f"RMSE(観測) = {rmse_obs:.4f}, RMSE(KF推定) = {rmse_kf:.4f}")

実行すると以下の結果が得られます。

K_1 = 0.2381, K_50 = 0.2207
RMSE(観測) = 1.7346, RMSE(KF推定) = 1.4041

式(9)の通り、初期分散 \(\hat{\sigma}_1^2 = 1.0 + Q = 1.25\) から \(K_1 = 1.25/(1.25+4.0) = 0.2381\) が導かれ、実測値と一致しています。カルマンフィルタは観測ノイズの標準偏差が2.0（\(\sqrt{R}\) ）と大きいにもかかわらず、生の観測値のRMSE 1.735 に対して推定値のRMSEは1.404 と、約19%誤差を低減できています。下図は真値・観測値・推定値（±2σ区間つき）を示したものです。

1次元ランダムウォークの真の状態、観測値、カルマンフィルタによる推定値と±2σ信頼区間の推移

推定値（赤破線）が観測値（灰色点）より滑らかに真値（青線）へ追従し、±2σ帯もほぼ全ての真値を覆っていることが確認できます。1次元等速運動モデルによるより実践的な実装例はカルマンフィルタの理論とPython実装を参照してください。

拡張カルマンフィルタ (Extended Kalman Filter, EKF)

概要

KFを非線形システムに拡張した手法です。非線形な関数を、現在の状態推定値の周りで**線形近似（テーラー展開の1次項まで利用）**することで、KFの枠組みを適用します。

方法

プロセスモデル \(f\) と観測モデル \(h\) を、ヤコビ行列（偏微分）を用いて線形化します。

\[ F_t = \left.\frac{\partial f}{\partial x}\right|_{x=\mu_{t-1}}, \qquad H_t = \left.\frac{\partial h}{\partial x}\right|_{x=\hat{\mu}_t} \]

この \(F_t\) と \(H_t\) を、線形カルマンフィルタの \(A\) と \(H\) の代わりに用いて計算を行いますが、予測ステップの式は非線形のまま計算します。

課題: 非線形性が強いシステムでは、線形化による近似誤差が大きくなり、推定精度が劣化したり、発散したりすることがあります。

導出：1次線形近似はいつ妥当か

EKFがなぜ「ヤコビ行列による線形化」で近似できるのか、テイラー展開から確認します。観測モデル \(h\) を現在の推定値 \(\hat{\mu}_t\) の周りでテイラー展開すると、

\[ h(x) = h(\hat{\mu}_t) + H_t (x - \hat{\mu}_t) + \frac{1}{2}(x-\hat{\mu}_t)^T \nabla^2 h(\hat{\mu}_t) (x-\hat{\mu}_t) + O\!\left(\|x-\hat{\mu}_t\|^3\right) \tag{10} \]

となります。ここで \(\nabla^2 h\) は \(h\) のヘッセ行列（2階微分）です。EKFは式(10)の右辺第1項（定数項）と第2項（1次のヤコビ項）だけを残し、2次以降の曲率項を丸ごと無視します。この近似が妥当なのは、状態の誤差共分散 \(\hat{\Sigma}_t\) が十分小さく、\(h\) の曲率 \(\nabla^2 h\) が推定値近傍でほぼ一定とみなせる場合です。実際、\(x\) の分布に関する期待値をとると、

\[ \mathbb{E}[h(x)] \approx h(\hat{\mu}_t) + \frac{1}{2}\,\mathrm{tr}\!\left(\nabla^2 h(\hat{\mu}_t)\, \hat{\Sigma}_t\right) \]

となり、EKFが無視する2次項の大きさは曲率 \(\nabla^2 h\) と共分散 \(\hat{\Sigma}_t\) の積のトレースで決まります。したがって、共分散が大きい（不確かさが大きい）場合や、非線形性が強く曲率が大きい場合には、この無視した項が無視できない大きさになり、EKFの推定は系統的なバイアスを持ち、共分散を過小評価して発散に至ることがあります。この線形化誤差と発散の具体例、およびニューラルネットワークで曲率をデータから学習して補正する近年の拡張手法は拡張カルマンフィルタ（EKF）の理論とPython実装で詳しく扱っています。

Unscentedカルマンフィルタ (Unscented Kalman Filter, UKF)

概要

EKFと同様に非線形システムを扱いますが、関数を直接線形化するのではなく、Unscented変換という手法で状態の確率分布を扱います。

方法

現在の状態分布を表現する少数の代表点（シグマ点）をサンプリングし、それぞれの点を非線形関数に通します。変換後の点の分布から、重み付き平均と共分散を再計算することで、EKFよりも高精度な推定を実現します。

利点: ヤコビ行列の計算が不要で、EKFよりも非線形性の強いシステムに対して頑健です。

直感：なぜUKFは2次まで精度が上がるのか

EKFが1次のテイラー項までしか捉えられなかったのに対し、UKFが2次の項まで捉えられる理由を、スカラー（1変数）の場合で見てみます。状態 \(x \sim \mathcal{N}(\mu, \sigma^2)\) を非線形関数 \(g\) で変換した後の真の期待値は、\(g\) を \(\mu\) の周りでテイラー展開して期待値をとると、

\[ \mathbb{E}[g(x)] = g(\mu) + \frac{1}{2}g''(\mu)\sigma^2 + \frac{1}{6}g'''(\mu)\,\mathbb{E}[(x-\mu)^3] + O(\sigma^4) \]

となります。ガウス分布は左右対称なので奇数次の中心モーメント \(\mathbb{E}[(x-\mu)^3]=0\) が消え、

\[ \mathbb{E}[g(x)] = g(\mu) + \frac{1}{2}g''(\mu)\sigma^2 + O(\sigma^4) \tag{11} \]

が残ります。EKFの予測ステップは実質的に \(g(\mu)\) だけを使う（式(10)と同様に1次項までの線形化を平均に適用する）ため、式(11)の第2項 \(\frac{1}{2}g''(\mu)\sigma^2\) という曲率由来のバイアスをまるごと落としてしまいます。

これに対しUKFは、\(\mu\) と \(\mu \pm \sqrt{(n+\lambda)}\,\sigma\) （\(n\) は状態次元、\(\lambda\) はスケーリングパラメータ）に置いたシグマ点を、\(g\) を線形近似せずにそのまま通し、重み \(W_i\) で加重平均します。シグマ点と重みは、変換前の分布の平均・分散を厳密に再現するように設計されているため、加重平均 \(\sum_i W_i\, g(\chi_i)\) を各シグマ点周りでテイラー展開すると、式(11)の \(\frac12 g''(\mu)\sigma^2\) の項まで自動的に再現されます。さらに左右対称にシグマ点を配置することで3次のモーメントに対応する項も相殺されるため、UKFは非線形性を持つ任意の分布に対して2次まで、ガウス分布を仮定できる場合はおおむね3次まで、テイラー展開の意味で正確に平均・共分散を捉えられます（Julier & Uhlmann, 2004）。ヤコビ行列という微分情報を陽に計算する代わりに、有限個の点をサンプリングして通すことで暗黙的に曲率の効果を取り込む、という発想の転換がUKFの核心です。多変数の場合のシグマ点配置・重みの設計とPython実装の詳細は Unscented Transformation（アンセンテッド変換）のPython実装および Unscented Kalman Filter（UKF）の理論とPython実装を参照してください。

粒子フィルタ (Particle Filter, PF)

概要

非線形・非ガウスの状態空間モデルを対象とした、より汎用的なフィルタリング手法です。モンテカルロ法に基づいています。

方法

状態の確率分布を、パーティクルと呼ばれる多数のサンプル点の集合で近似します。各パーティクルが状態の「仮説」を表しており、それぞれの尤度に基づいて重みが付けられます。

PFは主に「予測」「更新」「リサンプリング」のステップで構成されます。

予測: 全てのパーティクルをプロセスモデルに従って時間発展させます。
更新: 観測値が得られると、各パーティクルの尤度（観測値らしさ）を計算し、重みを更新します。
リサンプリング: 重みに応じてパーティクルを再サンプリングします。これにより、尤度の低い（ありえない）パーティクルは消滅し、尤度の高いパーティクルが複製され、推定が効率的に行われます。

利点: ガウス分布以外の複雑な確率分布も表現できるため、非常に汎用性が高いです。
課題: 状態の次元が大きくなると、分布を適切に表現するために必要なパーティクル数が指数関数的に増大する「次元の呪い」という問題があります。リサンプリング手法（系統リサンプリング・層化リサンプリングなど）の比較とPython実装は粒子フィルタのPython実装：リサンプリング手法の比較で扱っています。

近年の研究動向

ここまで見てきたKF・EKF・UKF・粒子フィルタはいずれも「モデル（\(f, h, Q, R\) ）が既知である」ことを前提とした古典的なベイズフィルタです。2023年以降は、これらの逐次フィルタリングの枠組みとニューラルネットワークを組み合わせ、モデルの一部または全部をデータから学習する方向の研究が活発化しています。

Chen and Li (2023) は “An overview of differentiable particle filters for data-adaptive sequential Bayesian inference”（Foundations of Data Science 誌、 arXiv:2302.09639 ）で、この方向性を粒子フィルタについて体系的にまとめています。従来の粒子フィルタは提案分布や尤度関数を人手で設計する必要がありましたが、リサンプリング処理を勾配が流れるように微分可能な形で再定義し、提案分布・動特性モデル・観測モデルをニューラルネットワークでパラメータ化することで、フィルタ全体をデータから end-to-end に学習する「微分可能な粒子フィルタ（differentiable particle filters）」が近年の主要な研究テーマとなっていることをこのサーベイは示しています。

この流れは粒子フィルタに限らず、 EKFの理論とPython実装で紹介しているNeural EKFなど、EKFやUKFにも共通して見られる傾向であり、「線形化やシグマ点によるモデルベースの近似」と「ニューラルネットワークによるデータ駆動の補正」を組み合わせるハイブリッドなアプローチが、今後のフィルタリング研究の主流になりつつあります。

カルマンフィルタの理論とPython実装 - 線形カルマンフィルタのアルゴリズムを1次元等速運動モデルでPython実装しています。
拡張カルマンフィルタ（EKF）の理論とPython実装 - ヤコビ行列による線形化で非線形システムに対応するEKFをPython実装しています。
Unscented Kalman Filter（UKF）の理論とPython実装 - シグマ点による非線形変換でEKFより高精度な推定を実現するUKFをPython実装しています。
Unscented Transformation（アンセンテッド変換）のPython実装 - UKFの基盤となるシグマ点による非線形変換の推定手法をPython実装しています。
Cubature Kalman Filter（CKF）の理論とPython実装 - UKFの重みの問題を解決するキュバチャカルマンフィルタをPython実装しています。
粒子フィルタのPython実装：リサンプリング手法の比較 - 各種リサンプリング手法を比較しながら粒子フィルタをPython実装しています。
カルマンスムーザ（RTS Smoother）の理論とPython実装 - カルマンフィルタの結果を全観測データで平滑化するRTSスムーザを解説しています。
指数移動平均（EMA）フィルタの周波数特性 - EMAフィルタのZ変換による周波数特性解析をPythonで可視化しています。
移動平均フィルタの種類と比較 - 単純移動平均・加重移動平均・EMAの特性比較とPython実装を解説しています。
高速フーリエ変換（FFT）の仕組みとPython実装 - フィルタの周波数特性を理解する基盤となるDFT/FFTのアルゴリズムと周波数解析を解説しています。
適応フィルタの基礎理論とデジタル信号処理への応用 - FIR/IIRフィルタの基礎からLMSアルゴリズムまで解説しています。
ローパスフィルタの設計と比較 - 移動平均・バターワース・チェビシェフフィルタの周波数応答を比較しています。