デジタルフィルタ
デジタルフィルタは、サンプリングされた時系列信号（例: \(x_0, x_1, \dots, x_N\) ）に対して処理を施すものです。アナログフィルタとは異なり、デジタル信号処理によって実現されます。代表的なデジタルフィルタに線形フィルタがあります。
線形フィルタ
線形フィルタは、入力信号の線形結合として出力信号を生成します。
FIR (Finite Impulse Response) フィルタ
FIRフィルタは、現在の入力と過去の有限個の入力サンプルのみを用いて出力を計算します。インパルス応答が有限の長さで終わるため、常に安定です。
\[ y_k = w_0 x_k + w_1 x_{k-1} + \dots + w_N x_{k-N} \]ここで \(w_i\) はフィルタ係数です。
IIR (Infinite Impulse Response) フィルタ
IIRフィルタは、現在の入力と過去の有限個の入力サンプルに加え、過去の有限個の出力サンプルも用いて出力を計算します。インパルス応答が無限に続く可能性があるため、設計によっては不安定になることがあります。
\[ y_k = w_0 x_k + w_1 x_{k-1} + \dots + w_N x_{k-N} - v_1 y_{k-1} - \dots - v_M y_{k-M} \]ここで \(w_i, v_i\) はフィルタ係数です。
適応フィルタ
適応フィルタは、入力信号の統計的性質が時間とともに変化する場合や、ノイズの特性が未知である場合など、環境の変化に合わせてフィルタ係数を自動的に調整（学習）するデジタルフィルタです。FIRフィルタとIIRフィルタのどちらの構造に対しても、同様の考え方でフィルタ係数を設計できます。本記事は 適応フィルタ（LMS/RLS）の理論とPython実装 など後続記事が「式(2)の導出過程は本記事で扱っています」と参照する、Wiener-Hopf方程式・MMSE理論の導出元（理論的な出発点）です。LMS/NLMS/RLSの具体的な実装・収束条件の証明・ノイズキャンセレーション応用は姉妹記事（ LMS／NLMSアルゴリズムの理論とPython実装 、 RLS（逐次最小二乗法）アルゴリズムの理論とPython実装 ）に譲り、本記事はMMSE最適化問題としてのWiener-Hopf方程式の完全な導出と、そこから生じる3つの根本的な論点（悪条件性・適応アルゴリズムとの関係・非定常性）に焦点を当てます。
適応フィルタの入出力と誤差
- 入力: \(x_k = x(kT)\) （\(k=0, 1, 2, \dots\) ）
- 出力: \(y_k = w_0 x_k + w_1 x_{k-1} + \dots + w_{N-1} x_{k-N+1}\) （ここではFIRフィルタの例を示していますが、IIRフィルタでも同様に考えられます）
- 目標値: \(d_k\) （理想的な出力信号）
- 誤差: \(e_k = d_k - y_k\)
フィルタ係数ベクトル \(W\) と入力信号ベクトル \(X_k\) を次のように定義します。
\[ W = [w_0, w_1, \dots, w_{N-1}]^T \] \[ X_k = [x_k, x_{k-1}, \dots, x_{k-N+1}]^T \]すると、フィルタ出力 \(y_k\) は \(y_k = W^T X_k = X_k^T W\) と書け、誤差 \(e_k\) は次のようになります。
\[ e_k = d_k - W^T X_k \]フィルタ係数の設計：MMSE問題としてのWiener-Hopf方程式の完全導出
適応フィルタの目的は、平均二乗誤差（MSE）を最小化する最適なフィルタ係数 \(W^o\) を見つけることです。これは次のMMSE（Minimum Mean Square Error）最適化問題として定式化できます。
\[ \min_{W} \; J(W) = \min_{W} \; \mathbb{E}\bigl[e_k^2\bigr] = \min_{W} \; \mathbb{E}\Bigl[\bigl(d_k - W^T X_k\bigr)^2\Bigr] \tag{1} \]コスト関数の展開
式(1)の2乗を展開します。
\[ J(W) = \mathbb{E}\bigl[d_k^2\bigr] - 2\, \mathbb{E}\bigl[d_k\, W^T X_k\bigr] + \mathbb{E}\bigl[(W^T X_k)^2\bigr] \]\(W\) は確率変数ではない（期待値の外に出せる）ことに注意して整理すると、
\[ J(W) = \mathbb{E}\bigl[d_k^2\bigr] - 2\, W^T \mathbb{E}\bigl[d_k X_k\bigr] + W^T \, \mathbb{E}\bigl[X_k X_k^T\bigr] \, W \]ここで、次の3つの統計量を定義します。
- \(\sigma_d^2 = \mathbb{E}[d_k^2]\) ：目標信号の分散（パワー）
- \(R = \mathbb{E}[X_k X_k^T]\) ：入力信号の自己相関行列（\(N \times N\) ）
- \(p = \mathbb{E}[d_k X_k]\) ：目標値と入力信号の相互相関ベクトル（\(N \times 1\) ）
これらを用いると、コスト関数は \(W\) についての2次形式として書けます。
\[ J(W) = \sigma_d^2 - 2 W^T p + W^T R W \tag{2} \]\(R\) は定義から常に半正定値（任意の \(a\) に対して \(a^T R a = \mathbb{E}[(a^T X_k)^2] \ge 0\) ）なので、\(J(W)\) は \(W\) の凸関数であり、勾配がゼロになる点は大域最小点です（\(R\) が正定値であれば、その最小点は一意に定まります）。
勾配をゼロと置く
式(2)を \(W\) で微分します。\(\nabla_W (W^T p) = p\) 、\(\nabla_W (W^T R W) = 2RW\) （\(R\) は対称行列）であることを使うと、
\[ \nabla_W J(W) = -2p + 2RW \]これをゼロと置くと、最適なフィルタ係数 \(W^o\) が満たすべき方程式が得られます。
\[ R W^o = p \tag{3} \]この関係式をWiener-Hopf方程式（正規方程式）と呼びます。\(R\) が正則（逆行列を持つ）であれば、\(W^o\) は解析的に求まります。
\[ W^o = R^{-1} p \tag{4} \]最小MSE \(J_{\min}\)
式(4)を式(2)に代入すると、Wiener-Hopf方程式を満たす最適解における最小MSEが得られます。\(W^{oT} R W^o = W^{oT} p\) （式(3)の両辺に \(W^{oT}\) を左からかけた関係）を使うと、
\[ J_{\min} = J(W^o) = \sigma_d^2 - 2 W^{oT} p + W^{oT} R W^o = \sigma_d^2 - 2 W^{oT} p + W^{oT} p = \sigma_d^2 - W^{oT} p \]\(W^o = R^{-1}p\) を代入すれば、
\[ J_{\min} = \sigma_d^2 - p^T R^{-1} p \tag{5} \]これは、目標信号のパワー \(\sigma_d^2\) から、入力信号によって「説明できる」分（\(p^T R^{-1} p\) ）を差し引いた残差が、達成可能な最小の誤差パワーであることを意味します。\(R\) が正定値であれば \(p^T R^{-1} p \ge 0\) かつ \(J_{\min} \ge 0\) が常に成り立ちます。
エッジケース1：\(R\) が特異・悪条件（ill-conditioned）な場合
式(4)の \(W^o = R^{-1}p\) は \(R\) が正則であることを前提としています。しかし実際の応用では、フィルタのタップ（入力）どうしが強く相関している、あるいは冗長な入力を含む場合に \(R\) が特異または悪条件になり、Wiener解が非一意になったり数値的に不安定になったりします。これを、8個の入力タップのうち複数が「ほぼ同じ信号（冗長なセンサ）」であるという典型的な状況で数値的に確認します。
import numpy as np
np.random.seed(0)
N = 3000
# 2つの「真に独立な」情報源
s1 = np.random.randn(N)
s2 = np.random.randn(N)
# 8タップの入力のうち、タップ2,3はタップ0(s1)のほぼ複製、
# タップ5,6はタップ1(s2)のほぼ複製という冗長な構成（例：オーバーサンプリングされたセンサ群）
eps = 1e-3
X_full = np.stack([
s1,
s2,
s1 + eps * np.random.randn(N),
s1 + eps * np.random.randn(N),
0.5 * s1 + 0.5 * s2,
s2 + eps * np.random.randn(N),
s2 + eps * np.random.randn(N),
0.3 * s1 + 0.7 * s2 + eps * np.random.randn(N),
], axis=1) # shape (N, 8)
# 真のシステムは独立な2つの情報源のみを使う
h_true = np.array([0.8, -0.5])
v = np.random.randn(N) * 0.05
d = h_true[0] * s1 + h_true[1] * s2 + v
M = X_full.shape[1]
R = (X_full.T @ X_full) / N
p = (X_full.T @ d) / N
print(f"cond(R) = {np.linalg.cond(R):.4e}")
print("固有値 =", np.array2string(np.linalg.eigvalsh(R), precision=6))
実行結果は次の通りでした。
cond(R) = 1.3337e+16
固有値 = [4.957129e-16 3.013686e-07 3.861934e-07 9.736620e-07 9.878075e-07
1.008506e-06 2.936495e+00 3.935144e+00]
条件数は \(1.33 \times 10^{16}\) と、倍精度浮動小数点の限界（約 \(10^{16}\) ）に達するほど悪条件です。固有値も8個中6個がほぼゼロ（\(10^{-7} \sim 10^{-16}\) オーダー）に潰れており、\(R\) の実効的な階数はほぼ2（真に独立な情報源の数）しかないことが確認できます。
try:
w_opt = np.linalg.solve(R, p)
print("solve() は成功:", np.array2string(w_opt, precision=4))
except np.linalg.LinAlgError as e:
print("solve() は失敗:", e)
w_lstsq, *_ = np.linalg.lstsq(R, p, rcond=None)
print("lstsq() の最小ノルム解:", np.array2string(w_lstsq, precision=4))
solve() は成功: [-3.9873 -5.4822 -0.4334 0.6601 9.2412 1.1976 -0.695 -0.2011]
lstsq() の最小ノルム解: [ 0.6714 -0.8236 -0.4334 0.6601 -0.0761 1.1976 -0.695 -0.2011]
ここが実務上もっとも危険な点です。np.linalg.solve はエラーを出さずに「成功」しますが、得られた係数（\(w_4 \approx 9.24\)
など）は真の関係（独立な2成分のみ）とかけ離れた無意味な値です。\(R\)
が完全な特異行列（ランク落ち）であれば LinAlgError: Singular matrix で明示的に失敗しますが、悪条件（特異に極めて近いが数値的には正則）な場合はエラーなく誤った解を返すため、条件数を確認しない限り異常に気づけません。
続いて、観測データにごく微小な摂動（標準偏差 \(10^{-6}\) ）を200回加えて解の感度を確認します。
rng = np.random.default_rng(1)
deltas = []
for _ in range(200):
Xp = X_full + rng.normal(0, 1e-6, size=X_full.shape)
Rp = (Xp.T @ Xp) / N
pp = (Xp.T @ d) / N
wp, *_ = np.linalg.lstsq(Rp, pp, rcond=None)
deltas.append(np.linalg.norm(wp - w_lstsq))
deltas = np.array(deltas)
print(f"摂動 1e-6 -> ||Δw|| の範囲: {deltas.min():.4e} - {deltas.max():.4e}（平均 {deltas.mean():.4e}）")
摂動 1e-6 -> ||Δw|| の範囲: 1.8197e+01 - 2.5618e+03（平均 6.9458e+02）
観測データへの \(10^{-6}\) スケールの微小な摂動が、解ベクトルを最大で \(2{,}562\) （入力摂動の実に \(10^9\) 倍以上）も動かしてしまいます。これは条件数が悪いほど解の相対誤差が \(\text{cond}(R)\) 倍まで増幅されうるという線形代数の一般的な性質の直接的な帰結です。
対処法は、\(R\) の対角に微小な正定数 \(\lambda\) を加えるリッジ正則化（\(R + \lambda I\) ）です。
for lam in [0, 1e-8, 1e-4, 1e-2, 1e-1]:
Rreg = R + lam * np.eye(M)
w_reg = np.linalg.solve(Rreg, p)
print(f"lambda={lam}: cond={np.linalg.cond(Rreg):.4e} ||w||={np.linalg.norm(w_reg):.4f}")
lambda=0: cond=1.3337e+16 ||w||=11.5730
lambda=1e-08: cond=3.9351e+08 ||w||=1.8992
lambda=0.0001: cond=3.9352e+04 ||w||=0.5412
lambda=0.01: cond=3.9451e+02 ||w||=0.5392
lambda=0.1: cond=4.0351e+01 ||w||=0.5232
\(\lambda=10^{-2}\) まで正則化すると条件数は \(394\) まで改善します。同じ手順で \(\lambda=10^{-2}\) での感度を再測定すると、\(\|\Delta w\|\) は \(6.3\times10^{-8} \sim 4.1\times10^{-7}\) まで低下し（正則化前と比べて実に10桁近い改善）、実務で使える安定性を回復できます。ただし正則化パラメータ \(\lambda\) を大きくしすぎると解自体が真値からバイアスを持つようになる（\(\|w\|\) が単調に縮小している点に注目）ため、条件数と推定バイアスのトレードオフを見ながら選ぶ必要があります。
エッジケース2：WienerフィルタとLMSの関係——適応フィルタは何に収束しているのか
式(4)の \(W^o = R^{-1}p\) は、\(R\) と \(p\) が既知であることを前提とした閉形式（バッチ処理）の解です。しかし実際の環境では信号の統計量は未知であるか、収集にコストがかかるため、この閉形式解を直接計算することは困難です。そこで登場するのが、勾配降下法に基づく逐次アルゴリズムです。
\[ W_{k+1} = W_k - \eta \, \nabla_W J(W_k) = W_k - \eta\,(-2p + 2RW_k) = W_k + 2\eta\,(p - RW_k) \tag{6} \]これは最急降下法と呼ばれ、\(R\) と\(p\) が既知であれば式(4)の \(W^o\) に指数関数的に収束します。しかし\(R, p\) が未知の場合、式(6)の勾配 \(\nabla_W J(W_k) = -2\mathbb{E}[e_k X_k]\) を計算する期待値を、瞬時的なサンプル \(-2 e_k X_k\) で置き換えた近似が LMS (Least Mean Squares) アルゴリズムです。
\[ W_{k+1} = W_k + \eta\, e_k X_k \tag{7} \]LMSは式(4)の \(W^o\) を直接計算せず、観測データ \(\{X_k, d_k\}\) だけから確率的勾配降下によって同じ解に収束することを期待するアルゴリズムです。この「同じ解に収束する」という主張を、定常な合成システム同定問題で数値的に検証します。
import numpy as np
np.random.seed(42)
N = 4000
M = 4
x = np.random.randn(N) # 白色雑音入力（定常）
h_true = np.array([1.0, -0.6, 0.3, -0.1])
v = np.random.randn(N) * 0.1 # 観測雑音
d = np.zeros(N)
for n in range(M - 1, N):
d[n] = h_true @ x[n - M + 1 : n + 1][::-1] + v[n]
# 閉形式のWiener解（式(4)：バッチ処理、全データから計算）
X = np.zeros((N - M + 1, M))
for i in range(M):
X[:, i] = x[M - 1 - i : N - i]
dd = d[M - 1 :]
R = (X.T @ X) / X.shape[0]
p = (X.T @ dd) / X.shape[0]
w_wiener = np.linalg.solve(R, p)
Jmin = np.mean(dd**2) - w_wiener @ p
print("Wiener解 w* =", np.array2string(w_wiener, precision=4))
print("真値 h_true =", h_true)
print(f"Jmin = {Jmin:.6f}（観測雑音分散 sigma_v^2 = {0.1**2}）")
# LMS（逐次処理：R, pを一切使わず観測データのみから更新）
mu = 0.01
w = np.zeros(M)
w_hist = np.zeros((N, M))
for n in range(M - 1, N):
xv = x[n - M + 1 : n + 1][::-1]
e = d[n] - w @ xv
w = w + mu * e * xv
w_hist[n] = w
print("LMS最終 w(N) =", np.array2string(w, precision=4))
print(f"||w_LMS(N) - w_wiener|| = {np.linalg.norm(w - w_wiener):.6e}")
実行結果は次の通りでした。
Wiener解 w* = [ 0.9989 -0.5984 0.3015 -0.1022]
真値 h_true = [ 1. -0.6 0.3 -0.1]
Jmin = 0.010337（観測雑音分散 sigma_v^2 = 0.010000000000000002）
LMS最終 w(N) = [ 1.0006 -0.5947 0.3051 -0.1077]
||w_LMS(N) - w_wiener|| = 7.692268e-03
まず、式(4)で計算した閉形式のWiener解 \(w^*=[0.9989, -0.5984, 0.3015, -0.1022]\) は真値 \(h_{\text{true}}=[1.0, -0.6, 0.3, -0.1]\) にほぼ一致し、\(J_{\min}=0.010337\) は観測雑音の分散 \(\sigma_v^2=0.01\) にほぼ等しくなっています（これは式(5)の理論通り、達成可能な最小誤差が観測雑音レベルに一致することを裏付けます）。そして、\(R\) や\(p\) を一度も明示的に計算していないLMSの最終係数 \([1.0006, -0.5947, 0.3051, -0.1077]\) も、Wiener解とのユークリッド距離 \(7.7\times10^{-3}\) という小さな差にまで収束しています。
反復回数ごとの重み誤差 \(\|W_k - W^o\|\) を追跡すると、収束の過程がより明確になります。
| 反復回数 \(k\) | \(\|W_k - W^o\|\) |
|---|---|
| 50 | 0.81228 |
| 100 | 0.51381 |
| 200 | 0.19565 |
| 500 | 0.01374 |
| 1000 | 0.01389 |
| 2000 | 0.01747 |
| 4000 | 0.00769 |
初期値 \(W_0=0\) から数百反復でWiener解の近傍まで急速に収束し、その後は完全にゼロへは収束せず \(10^{-2}\) オーダーの残留誤差（過剰MSE／ミスアジャストメント、瞬時勾配のばらつみに起因）の周りで揺らぎ続けています。この揺らぎ自体はLMSの原理的な性質であり、ステップサイズ \(\eta\) を小さくすれば揺らぎは減りますが収束速度も遅くなるというトレードオフがあります（この収束速度・定常誤差の定量的な理論はステップサイズの収束条件とあわせて LMS／NLMSアルゴリズムの理論とPython実装 で扱っています）。
エッジケース3：非定常環境でのWiener解の限界とLMSによる追従
式(4)のWiener解 \(W^o=R^{-1}p\) は、入力・目標信号の統計量（\(R, p\) ）が時間不変（定常）であることを暗黙に仮定しています。しかし実世界のシステム（音響エコー経路の変化、通信チャネルのフェージングなど）は時間とともに変化する（非定常）ことが多く、この場合「一度だけ計算した固定のWiener解」は、統計量が変化した後は最適ではなくなります。一方、LMSのような適応アルゴリズムは、逐次的に係数を更新し続けるため、統計量の変化を追従できる可能性があります。これを、未知システムが途中で急変するシステム同定問題で検証します。
import numpy as np
np.random.seed(7)
N = 6000
M = 4
change_point = N // 2 # n=3000でシステムが急変
x = np.random.randn(N)
h1 = np.array([1.0, -0.6, 0.3, -0.1])
h2 = np.array([-0.5, 0.8, -0.4, 0.2]) # 急変後の別システム
v = np.random.randn(N) * 0.1
d = np.zeros(N)
for n in range(M - 1, N):
xv = x[n - M + 1 : n + 1][::-1]
h = h1 if n < change_point else h2
d[n] = h @ xv + v[n]
def build_R_p(x, d, M, n0, n1):
X = np.zeros((n1 - n0, M))
dd = np.zeros(n1 - n0)
for k, n in enumerate(range(n0, n1)):
X[k] = x[n - M + 1 : n + 1][::-1]
dd[k] = d[n]
R = (X.T @ X) / X.shape[0]
p = (X.T @ dd) / X.shape[0]
return R, p
# 急変前のデータのみから計算した「固定Wiener解」
R1, p1 = build_R_p(x, d, M, M - 1, change_point)
w_fixed = np.linalg.solve(R1, p1)
print("固定Wiener解（前半データのみ） w_fixed =", np.array2string(w_fixed, precision=4))
print("h1 =", h1, " h2 =", h2)
# 固定解を全区間（前半・後半とも）にそのまま適用
e_fixed = np.zeros(N)
for n in range(M - 1, N):
xv = x[n - M + 1 : n + 1][::-1]
e_fixed[n] = d[n] - w_fixed @ xv
# 適応LMSを全区間にわたって実行
mu = 0.02
w = np.zeros(M)
e_lms = np.zeros(N)
for n in range(M - 1, N):
xv = x[n - M + 1 : n + 1][::-1]
e_lms[n] = d[n] - w @ xv
w = w + mu * e_lms[n] * xv
seg_pre = slice(M - 1, change_point)
seg_post_early = slice(change_point, change_point + 500)
seg_post_late = slice(N - 500, N)
for name, e in [("固定Wiener", e_fixed), ("LMS", e_lms)]:
mse_pre = np.mean(e[seg_pre] ** 2)
mse_early = np.mean(e[seg_post_early] ** 2)
mse_late = np.mean(e[seg_post_late] ** 2)
print(f"{name:10s}: 急変前MSE={mse_pre:.5f} 急変直後MSE(500点)={mse_early:.5f} "
f"十分後MSE(末尾500点)={mse_late:.5f}")
print(f"\nLMS最終 w(N) = {np.array2string(w, precision=4)}")
print(f"||w_LMS(N) - h2|| = {np.linalg.norm(w - h2):.5f}")
print(f"||w_fixed - h2|| = {np.linalg.norm(w_fixed - h2):.5f}")
実行結果は次の通りでした。
固定Wiener解（前半データのみ） w_fixed = [ 1.0011 -0.5988 0.3002 -0.0979]
h1 = [ 1. -0.6 0.3 -0.1] h2 = [-0.5 0.8 -0.4 0.2]
固定Wiener : 急変前MSE=0.00966 急変直後MSE(500点)=4.34709 十分後MSE(末尾500点)=5.09337
LMS : 急変前MSE=0.02295 急変直後MSE(500点)=0.25829 十分後MSE(末尾500点)=0.01016
LMS最終 w(N) = [-0.5012 0.8011 -0.3913 0.1942]
||w_LMS(N) - h2|| = 0.01060
||w_fixed - h2|| = 2.18833
急変前のデータから計算した固定Wiener解 \(w_{\text{fixed}}=[1.0011, -0.5988, 0.3002, -0.0979]\) は \(h_1\) にほぼ一致し、急変前区間ではMSE\(=0.00966\) （雑音分散 \(0.01\) とほぼ一致、理論上の最良値）を達成しています。しかしシステムが \(h_2\) に急変した後は、\(w_{\text{fixed}}\) は \(h_2\) から \(\|w_{\text{fixed}}-h_2\|=2.188\) もかけ離れたままとなり、MSEは急変直後で\(4.347\) 、十分時間が経った後も\(5.093\) と、固定解のままでは全く改善しないまま高止まりします。
一方、適応LMS（\(\mu=0.02\) ）は急変前区間ではMSE\(=0.02295\) と固定解よりやや劣ります（LMSの原理的な残留誤差のため）が、システム急変直後は一時的にMSEが\(0.258\) まで跳ね上がるものの、その後は継続的な係数更新によって\(h_2\) を再学習し、末尾500点でのMSEは\(0.01016\) （雑音床とほぼ同水準）まで回復します。最終的な重みも \(\|w_{\text{LMS}}(N)-h_2\|=0.0106\) と、\(h_2\) にほぼ完全に一致しています。
この結果は、非定常環境において「一度計算したWiener解」と「継続的に更新し続ける適応フィルタ」の本質的な違いを定量的に示しています。Wiener解はある瞬間の統計量に対する最適解にすぎず、統計量自体が変化する環境では自動的に追従しません。適応フィルタが実務で選ばれる最大の理由は、単にWiener解を「学習で近似できる」からではなく、Wiener解そのものが時間とともに動く場合にそれを追いかけられるからです。
最新研究動向
Wiener-Hopf方程式・MMSE線形フィルタの理論は半世紀以上前に確立されていますが、その悪条件性への対処や現代的な応用は現在も活発に研究されています。
- 正則化の自動決定: Zanco, Szczecinski, & Benesty (2023/2024, arXiv:2312.06560) は、本記事のエッジケース1で扱った「\(R\) の悪条件性に対するリッジ正則化」について、正則化パラメータ \(\lambda\) を観測信号から自動的に決定するベイズ統計的な枠組みを提案し、システム同定・ビームフォーミング応用で最適に近い性能を達成することを示しています。手動でグリッドサーチする代わりに \(\lambda\) をデータ駆動で決める点が本記事の実験に対する直接の発展です。
- Wiener解とKalman/LMS/NLMSの統一的理解: Szczecinski, Benesty, & Kuhn (2025, arXiv:2502.18325) は、ベイズ逐次推定の枠組みから適応フィルタを再導出し、ガウス雑音を仮定するとLMS・NLMS・Kalmanフィルタが、ラプラス分布などの非ガウス雑音を仮定すると符号誤差法などのロバストな亜種が自然に導かれることを示しました。本記事のエッジケース2（LMSはWiener解に何を目指して収束するのか）を、より広い統計的仮定のもとで捉え直す視点を提供します。
- 古典的Wienerフィルタの再評価: Bled & Pitié (2023, arXiv:2303.16640) は、画像ノイズ除去において実装を注意深く最適化したWienerフィルタが、DnCNNのような深層学習ベースの手法に匹敵する性能を達成できることを示しました。半世紀以上前の線形MMSE理論が、深層学習全盛の現在でも侮れない実用性を持つことを裏付ける結果です。
関連記事
- 信号処理におけるフィルタリング手法の基礎 - カルマンフィルタ、EKF、UKF、粒子フィルタの概要を解説しています。
- 指数移動平均（EMA）フィルタの周波数特性 - IIRフィルタの一種であるEMAの周波数特性をPythonで可視化しています。
- ローパスフィルタの設計と比較 - 移動平均・バターワース・チェビシェフフィルタの周波数応答を比較しています。
- FIRフィルタとIIRフィルタの比較 - 本記事で触れたデジタルフィルタの2大分類（FIR/IIR）の特性・設計思想を体系的に比較しています。
- チェビシェフフィルタの設計原理とPython実装 - 等リップル特性を持つ代表的なIIRフィルタの設計・実装を解説しています。
- ベッセルフィルタの理論とPython実装 - 群遅延が最大限平坦なIIRフィルタの設計と、位相歪み最小化の観点を解説しています。
- バンドパスフィルタの設計とPython実装 - 帯域通過フィルタの設計（ローパス＋ハイパスの合成、周波数変換）を実装例とともに解説しています。
- ウィーナーフィルタ：最適線形フィルタの理論とPython実装 - 本記事のWiener-Hopf方程式を周波数領域（SNRベースのゲイン設計）に拡張し、因果性の制約・悪条件な\(R\) への正則化・Kalmanフィルタとの数値的等価性まで掘り下げた発展編です。
- 適応フィルタ（LMS/RLS）の理論とPython実装 - 本記事のWiener-Hopf方程式・最急降下法を出発点に、LMS・NLMS・RLSの比較とノイズキャンセレーション応用を扱った実践編です。
- LMS／NLMSアルゴリズムの理論とPython実装 - 本記事のLMSアルゴリズムについて、収束条件の厳密な証明・ミスアジャストメントの定量化まで踏み込んだ発展編です。
- RLS（逐次最小二乗法）アルゴリズムの理論とPython実装 - 本記事のWiener-Hopf正規方程式を行列反転補題で逐次的に厳密解にするRLSアルゴリズムと、カルマンフィルタとの数値的等価性を扱っています。
参考文献
- 大松 繁, “信号処理ーデジタルフィルタ”
- Haykin, S. (2014). Adaptive Filter Theory (5th ed.). Pearson. Chapters 1-2.
- Widrow, B., & Stearns, S. D. (1985). Adaptive Signal Processing. Prentice-Hall.
- Zanco, D. G. P., Szczecinski, L., & Benesty, J. (2023/2024). Automatic Regularization for Linear MMSE Filters. arXiv:2312.06560. https://arxiv.org/abs/2312.06560
- Szczecinski, L., Benesty, J., & Kuhn, E. V. (2025). A Unified Bayesian Perspective for Conventional and Robust Adaptive Filters. arXiv:2502.18325. https://arxiv.org/abs/2502.18325
- Bled, C., & Pitié, F. (2023). Pushing The Limits of the Wiener Filter in Image Denoising. arXiv:2303.16640. https://arxiv.org/abs/2303.16640