はじめに
ウィーナーフィルタは、平均二乗誤差（MSE）を最小化する最適線形フィルタです。ノイズを含む観測信号から元の信号を推定する問題において、定常確率過程の枠組みで理論的に最適な解を与えます。
問題設定
ノイズを含む観測信号:
\[y(t) = s(t) + n(t) \tag{1}\]ここで \(s(t)\) は所望の信号、\(n(t)\) はノイズです。フィルタ \(h(t)\) を通した出力 \(\hat{s}(t)\) で \(s(t)\) を推定します。
\[\hat{s}(t) = \int_{-\infty}^{\infty} h(\tau) y(t - \tau) d\tau \tag{2}\]MSE最小化とウィーナー・ホップ方程式の完全導出
推定誤差の二乗平均を最小化します。
\[\min_h E\left[|s(t) - \hat{s}(t)|^2\right] \tag{3}\]式(2)は無限次元の畳み込み積分ですが、実装上は有限タップのFIRフィルタで近似します。ここでは離散・有限次元のベクトル形式で、なぜ最適解が正規方程式（ウィーナー・ホップ方程式）として書けるのかを勾配計算から直接導出します。
ベクトル形式の定式化
観測信号を長さ \(M\) のタップ遅延線でベクトル化します。
\[ \mathbf{x}[n] = \bigl[y[n],\ y[n-1],\ \ldots,\ y[n-M+1]\bigr]^\top \]フィルタ係数ベクトル \(\mathbf{w} = [w_0, w_1, \ldots, w_{M-1}]^\top\) による推定値は
\[ \hat{s}[n] = \mathbf{w}^\top \mathbf{x}[n] \]であり、これは式(2)の積分をタップ数 \(M\) で打ち切った離散近似に相当します。MSEコスト関数は
\[ J(\mathbf{w}) = E\left[\bigl(s[n] - \mathbf{w}^\top \mathbf{x}[n]\bigr)^2\right] \tag{4} \]展開と勾配ゼロ条件
2乗を展開すると
\[ J(\mathbf{w}) = E[s[n]^2] - 2\mathbf{w}^\top E[\mathbf{x}[n]s[n]] + \mathbf{w}^\top E[\mathbf{x}[n]\mathbf{x}[n]^\top]\mathbf{w} \]ここで自己相関行列 \(\mathbf{R}_{xx} = E[\mathbf{x}[n]\mathbf{x}[n]^\top]\) （本記事の \(R_{yy}\) の行列版）、相互相関ベクトル \(\mathbf{r}_{xd} = E[\mathbf{x}[n]s[n]]\) （\(R_{sy}\) の行列版）、信号分散 \(\sigma_s^2 = E[s[n]^2]\) とおくと、コスト関数は \(\mathbf{w}\) の2次形式になります。
\[ J(\mathbf{w}) = \sigma_s^2 - 2\mathbf{w}^\top \mathbf{r}_{xd} + \mathbf{w}^\top \mathbf{R}_{xx}\mathbf{w} \tag{5} \]\(\mathbf{w}\) で微分すると
\[ \nabla_{\mathbf{w}} J(\mathbf{w}) = -2\mathbf{r}_{xd} + 2\mathbf{R}_{xx}\mathbf{w} \]勾配をゼロと置くと、正規方程式（ウィーナー・ホップ方程式の離散版）
\[ \mathbf{R}_{xx}\mathbf{w} = \mathbf{r}_{xd} \tag{6} \]が得られます。これを連続時間・無限タップに一般化したものが、直交性原理から導かれる次の畳み込み方程式です。
\[R_{yy}(\tau) * h(\tau) = R_{sy}(\tau) \tag{7}\]ここで \(R_{yy}\) は観測信号の自己相関、\(R_{sy}\) は所望信号と観測信号の相互相関です。式(6)と式(7)は同じ最適条件（推定誤差 \(s[n]-\hat{s}[n]\) が観測データ \(\mathbf{x}[n]\) のどの成分とも無相関になる「直交性原理」）を、有限次元と無限次元でそれぞれ表現したものです。
\(J(\mathbf{w})\) の凸性と\(\mathbf{R}_{xx}\) の半正定値性
式(6)が本当にMSEを最小化するかを確認するには、\(J(\mathbf{w})\) が凸であることを示す必要があります。勾配を再度 \(\mathbf{w}\) で微分するとヘッセ行列は
\[ \nabla^2_{\mathbf{w}} J(\mathbf{w}) = 2\mathbf{R}_{xx} \]となります。任意のベクトル \(\mathbf{a} \neq \mathbf{0}\) に対して
\[ \mathbf{a}^\top \mathbf{R}_{xx}\mathbf{a} = E\bigl[(\mathbf{a}^\top \mathbf{x}[n])^2\bigr] \ge 0 \]が成り立つため、\(\mathbf{R}_{xx}\) は常に半正定値です（2乗の期待値は非負なので）。したがってヘッセ行列 \(2\mathbf{R}_{xx}\) も半正定値であり、\(J(\mathbf{w})\) は \(\mathbf{w}\) について凸関数です。凸関数の停留点は大域最小点であるため、式(6)を満たす \(\mathbf{w}\) は必ずMSEを大域的に最小化します。
さらに \(\mathbf{R}_{xx}\) が正定値（固有値がすべて厳密に正）であれば最小点は一意に定まります。逆に \(\mathbf{x}[n]\) の成分間に線形従属な関係がある場合（後述する純粋な正弦波が典型例）は \(\mathbf{R}_{xx}\) が特異になり、解が一意でなくなったり数値的に不安定になったりします。
周波数領域での解の導出
式(7)は時間領域の畳み込み方程式であり、両辺をフーリエ変換すると畳み込み定理より積になります。
\[ \mathcal{F}\{R_{yy}(\tau) * h(\tau)\} = \mathcal{F}\{R_{yy}(\tau)\} \cdot \mathcal{F}\{h(\tau)\} \]ウィーナー・ヒンチンの定理より、自己相関関数（および相互相関関数）のフーリエ変換はパワースペクトル密度（およびクロススペクトル密度）に一致します。
\[ \mathcal{F}\{R_{yy}(\tau)\} = S_{yy}(f), \qquad \mathcal{F}\{R_{sy}(\tau)\} = S_{sy}(f), \qquad \mathcal{F}\{h(\tau)\} = H(f) \]したがって式(7)は周波数領域で
\[ S_{yy}(f) \, H(f) = S_{sy}(f) \]という、周波数ビンごとに独立なスカラー方程式に分解されます。両辺を \(S_{yy}(f)\) で割ると
\[H(f) = \frac{S_{sy}(f)}{S_{yy}(f)} \tag{8}\]が得られます。時間領域では \(M\) 元（あるいは無限次元）の連立方程式だった正規方程式が、周波数領域では周波数ごとの単純な除算に変わる点が、ウィーナーフィルタが周波数領域で簡潔に解ける理由です。
信号とノイズが無相関の場合:
\[H(f) = \frac{S_{ss}(f)}{S_{ss}(f) + S_{nn}(f)} = \frac{\text{SNR}(f)}{\text{SNR}(f) + 1} \tag{9}\]ここで \(S_{ss}(f)\) は信号のパワースペクトル、\(S_{nn}(f)\) はノイズのパワースペクトル、\(\text{SNR}(f) = S_{ss}(f) / S_{nn}(f)\) は周波数ごとの信号対雑音比です。これは \(y=s+n\) 、\(E[sn]=0\) のとき \(S_{yy}(f)=S_{ss}(f)+S_{nn}(f)\) 、\(S_{sy}(f)=S_{ss}(f)\) となることから式(8)より直ちに導かれます。
直感的な理解: SNRが高い周波数帯域はそのまま通し、SNRが低い帯域は減衰させます。
数値実験：ウィーナー・ホップ方程式の求解とMMSE最適性の検証
式(5)(6)の主張——正規方程式を解けば本当にMSEが最小化される——を、実際にnumpyで確認します。既存のPython実装と同じ信号（\(fs=1000\)
、\(5\,\mathrm{Hz}+20\,\mathrm{Hz}\)
の正弦波合成、標準偏差 \(0.8\)
の加法性白色ノイズ、np.random.seed(42)）から、タップ数 \(M=8\)
のタップ遅延線 \(\mathbf{x}[n]\)
と所望信号 \(d[n]=s[n]\)
（クリーンな信号）を構成し、経験的な \(\mathbf{R}_{xx}\)
、\(\mathbf{r}_{xd}\)
を計算して np.linalg.solve で式(6)を解きました。
import numpy as np
np.random.seed(42)
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*20*t)
noise = np.random.normal(0, 0.8, len(t))
observed = signal + noise
M = 8
N = len(observed)
X = np.zeros((N - M + 1, M))
for i in range(M):
X[:, i] = observed[M-1-i : N-i]
d = signal[M-1:]
R = (X.T @ X) / X.shape[0]
r = (X.T @ d) / X.shape[0]
w_opt = np.linalg.solve(R, r)
実行結果は次の通りです。
- 解 \(\mathbf{w}_{\text{opt}} = [0.1507,\ 0.1398,\ 0.1276,\ 0.1134,\ 0.1054,\ 0.0890,\ 0.0744,\ 0.0599]\) 、最小コスト \(J(\mathbf{w}_{\text{opt}}) = 0.088442\)
- 勾配ゼロの確認: \(\|\nabla J(\mathbf{w}_{\text{opt}})\|_\infty = 2.22\times10^{-16}\) （倍精度の丸め誤差レベルで実質ゼロ）
- \(\mathbf{R}_{xx}\) の正定値性: 固有値の最小値 \(0.5708\) 、最大値 \(5.7854\) （すべて正 \(\Rightarrow\) 正定値 \(\Rightarrow\) 解は一意）
- 最適性の直接検証: \(\mathbf{w}_{\text{opt}}\) に標準偏差 \(0.05\) のガウス摂動を加えた \(2{,}000\) 回の試行すべて（\(100.00\%\) ）で \(J(\mathbf{w}_{\text{opt}}+\boldsymbol{\delta}) \ge J(\mathbf{w}_{\text{opt}})\) となり、コスト増加量の最小値ですら \(1.69\times10^{-3} > 0\) でした
勾配が数値的にゼロであること、\(\mathbf{R}_{xx}\) が正定値であること、そしてどの方向に摂動を与えてもコストが増加しないことの3点が、式(6)の解が大域最小点であるという上の証明を実データで裏付けています。
Python実装
1次元信号のノイズ除去
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq
np.random.seed(42)
# --- 信号生成 ---
fs = 1000 # サンプリング周波数
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise = np.random.normal(0, 0.8, len(t))
observed = signal + noise
# --- ウィーナーフィルタ（周波数領域） ---
N = len(observed)
Y = fft(observed)
freqs = fftfreq(N, 1/fs)
# パワースペクトル推定
S_yy = np.abs(Y)**2 / N
noise_power = np.var(noise) * N # ノイズのパワー（一様）
S_nn = noise_power / N * np.ones(N)
S_ss = np.maximum(S_yy - S_nn, 0) # 信号パワーの推定
# ウィーナーフィルタ
H = S_ss / (S_ss + S_nn + 1e-10)
S_filtered = ifft(Y * H).real
# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(12, 8))
axes[0].plot(t, signal, 'b-', alpha=0.7, label='Original')
axes[0].set_title('Original Signal')
axes[0].set_ylabel('Amplitude')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[1].plot(t, observed, 'gray', alpha=0.5, label='Noisy')
axes[1].set_title('Observed (Signal + Noise)')
axes[1].set_ylabel('Amplitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[2].plot(t, signal, 'b-', alpha=0.3, label='Original')
axes[2].plot(t, S_filtered, 'r-', alpha=0.8, label='Wiener filtered')
axes[2].set_title('Wiener Filter Output')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
mse_noisy = np.mean((signal - observed)**2)
mse_filtered = np.mean((signal - S_filtered)**2)
print(f"MSE (noisy): {mse_noisy:.4f}")
print(f"MSE (filtered): {mse_filtered:.4f}")
print(f"Improvement: {(1 - mse_filtered/mse_noisy)*100:.1f}%")
実行すると MSE (noisy): 0.6133、MSE (filtered): 0.1284、Improvement: 79.1% という結果が得られます。ノイズ分散既知（真値を使用）という理想条件下では、周波数領域ウィーナーフィルタがMSEを約8割削減できることが確認できました。
フィルタの周波数応答
plt.figure(figsize=(10, 4))
positive_freqs = freqs[:N//2]
plt.plot(positive_freqs, H[:N//2], 'g-', linewidth=2)
plt.xlabel('Frequency [Hz]')
plt.ylabel('|H(f)|')
plt.title('Wiener Filter Frequency Response')
plt.grid(True, alpha=0.3)
plt.xlim(0, 100)
plt.tight_layout()
plt.show()
因果性の制約：非因果的フィルタとペイリー・ウィーナー条件
式(8)(9)で得た \(H(f)\)
は、一般には非因果的（未来の観測サンプルも使う）フィルタです。式(8)の解は周波数領域の代数方程式を満たすだけで、対応するインパルス応答 \(h(\tau)\)
が \(\tau<0\)
で恒等的にゼロである保証はありません。実際、本記事のPython実装（fft→乗算→ifftのバッチ処理）は信号全体を使うため非因果的であり、リアルタイム処理には使えません。
リアルタイム用途では、因果的ウィーナーフィルタ（\(h(\tau)=0\ (\tau<0)\) ）を別途構成する必要があります。標準的な手順は次のスペクトル分解（spectral factorization）です。
- \(S_{yy}(f)\) を、最小位相（=因果的にフィルタ・逆フィルタできる）成分 \(L(f)\) を使って \(S_{yy}(f) = \sigma^2\, L(f)\, L^\ast(f)\) の形に分解する。
- 因果的ウィーナーフィルタは
で与えられる。ここで \([\cdot]_+\) は逆フーリエ変換した結果の非負時刻（\(\tau \ge 0\) ）成分だけを残す演算を表す。
このスペクトル分解が可能であるための必要十分条件がペイリー・ウィーナーの条件です。
\[ \int_{-\infty}^{\infty} \frac{\bigl|\ln S_{yy}(f)\bigr|}{1+f^2}\,df < \infty \]直感的には、\(S_{yy}(f)\) がある帯域で恒等的にゼロになる（＝有限区間で完全に無音になる）ような病的なスペクトルでない限りこの条件は満たされます。実務上は、加法性ノイズがあれば \(S_{yy}(f) = S_{ss}(f)+S_{nn}(f) \ge S_{nn}(f) > 0\) となりこの条件は自動的に満たされることが多く、条件そのものよりも「因果的フィルタは非因果的フィルタよりMSEが大きい（または等しい）」という点が実務上重要です。非因果的フィルタは未来の情報を使える分だけ必ず有利（またはタイ）であり、因果的フィルタはその劣化を受け入れた近似解になります。カルマンフィルタは、この因果的ウィーナーフィルタを状態空間モデル上で逐次計算する具体的な手続きに相当します（後述）。
悪条件な自己相関行列と正則化
式(6)の凸性の議論で触れたとおり、\(\mathbf{x}[n]\) の成分が線形従属に近いと \(\mathbf{R}_{xx}\) は特異または悪条件（ill-conditioned）になり、解 \(\mathbf{w}=\mathbf{R}_{xx}^{-1}\mathbf{r}_{xd}\) が不安定になります。代表的な悪条件ケースが決定論的な正弦波信号です。純粋な正弦波はわずか2つの自由度（振幅・位相）しか持たないため、タップ数 \(M>2\) のFIRフィルタでは自己相関行列の階数が理論上たかだか2にしかなりません。
import numpy as np
fs = 1000
t = np.arange(0, 1, 1/fs)
pure_sine = np.sin(2*np.pi*10*t) # ノイズなしの純粋な正弦波
M = 5 # タップ数(2より大きい)
N = len(pure_sine)
X = np.zeros((N - M + 1, M))
for i in range(M):
X[:, i] = pure_sine[M-1-i : N-i]
R = (X.T @ X) / X.shape[0]
d = pure_sine[M-1:]
r = (X.T @ d) / X.shape[0]
cond = np.linalg.cond(R)
eigvals = np.linalg.eigvalsh(R)
実行すると、\(\text{cond}(\mathbf{R}_{xx}) = 7.13\times10^{17}\)
（倍精度の限界に近い極端な悪条件）、固有値は \([\,{-0.0},\ 0.0,\ 0.0,\ 0.0196,\ 2.4902\,]\)
となり、理論通り3つの固有値がほぼゼロ（数値誤差でわずかに負にもなる）でした。この \(\mathbf{R}_{xx}\)
に対して np.linalg.solve(R, r) を実行すると、実際に LinAlgError: Singular matrix で失敗します。np.linalg.lstsq を使えば最小ノルム解 \(\mathbf{w}=[0.5981,\ 0.3998,\ 0.2000,\ -0.0006,\ -0.2013]\)
（\(\|\mathbf{w}\|=0.7734\)
）は得られますが、この解は観測データへの微小な摂動に対して敏感です。実際、観測信号に標準偏差 \(10^{-6}\)
の無視できるほど小さいノイズを加えて \(\mathbf{R}_{xx}\)
・\(\mathbf{r}_{xd}\)
を作り直すだけで（100試行）、解の変化量 \(\|\Delta\mathbf{w}\|\)
は \(5.5\text{–}7.5\times10^{-5}\)
のオーダーに達しました（入力の摂動レベルの約60倍に増幅されています）。
対処法はリッジ回帰と同型の正則化です。\(\mathbf{R}_{xx}\) の対角に微小な正の定数 \(\lambda\) を加えた \(\mathbf{R}_{xx}+\lambda I\) を用いて解きます。
| \(\lambda\) | \(\text{cond}(\mathbf{R}_{xx}+\lambda I)\) | \(\|\mathbf{w}\|\) | 微小摂動への感度 \(\|\Delta\mathbf{w}\|\) |
|---|---|---|---|
| \(0\) | \(7.13\times10^{17}\) （特異） | 解けない | — |
| \(10^{-8}\) | \(2.49\times10^{8}\) | \(0.7734\) | \(5.5\text{–}7.5 \times 10^{-5}\) |
| \(10^{-4}\) | \(2.49\times10^{4}\) | \(0.7707\) | — |
| \(10^{-2}\) | \(2.50\times10^{2}\) | \(0.6099\) | \(2.8\times10^{-9}\text{–}4.4\times10^{-8}\) |
| \(10^{-1}\) | \(2.59\times10^{1}\) | \(0.4406\) | — |
\(\lambda=10^{-2}\) まで正則化を強めると条件数は \(2.5\times10^{2}\) まで改善し、同じ \(10^{-6}\) スケールの微小摂動に対する解の変化量は \(2.8\times10^{-9}\text{–}4.4\times10^{-8}\) と、入力の摂動よりさらに小さいオーダーまで収まります（\(\lambda=10^{-8}\) のときと比べて感度が3桁以上改善）。この結果からわかる実務上の指針は、「\(\mathbf{R}_{xx}\) の条件数を確認し、悪条件が疑われる場合（ほぼ周期的な信号、タップ数が信号の実効自由度より多い場合など）は正則化パラメータ \(\lambda\) を、性能劣化とのトレードオフを見ながら適切に選ぶ」ことです。
ノイズ統計量の推定誤差がフィルタ性能に与える影響
式(9)の \(H(f)=S_{ss}/(S_{ss}+S_{nn})\) は真のパワースペクトルを前提としますが、実際には \(S_{nn}\) （ノイズパワー）は有限データから推定するしかなく、推定誤差がフィルタ性能に直接影響します。広帯域信号（AR(1)過程 \(s_k = 0.9\,s_{k-1}+w_k\) 、\(w_k\sim\mathcal{N}(0,0.05)\) ）に分散 \(0.25\) の白色ノイズを加えた合成信号で、ノイズ分散の推定値を真値から意図的にずらして周波数領域ウィーナーフィルタを構成し、真の信号に対するMSEを測定しました。
| ノイズ分散の推定誤差 | MSE |
|---|---|
| \(-50\%\) | \(0.142980\) |
| \(-30\%\) | \(0.122002\) |
| \(0\%\) （真値） | \(0.101663\) |
| \(+30\%\) | \(0.089713\) |
| \(+50\%\) | \(0.084168\) |
| \(+100\%\) | \(0.076571\) |
| \(+200\%\) | \(0.076315\) （最小） |
| \(+300\%\) | \(0.083668\) |
| \(+500\%\) | \(0.101872\) |
| \(+1000\%\) | \(0.148345\) |
過小評価（\(-50\%\sim-30\%\) ）はMSEを単調に悪化させます。ノイズパワーを過小に見積もると \(H(f)\) が過大評価され、残留ノイズが通過してしまうためで、直感通りの挙動です。一方で意外なのは、ノイズ分散を過大に見積もった場合（\(+30\%\sim+200\%\) ）に、真値を使うよりMSEが改善している点です。これは、単一時刻のスナップショットから計算した周期グラム \(S_{yy}=|Y|^2/N\) 自体が分散の大きい（データ長を伸ばしても分散が減らない）不安定な推定量であること（詳細は 窓関数とパワースペクトル密度（PSD）の理論とPython実装 を参照）に起因します。ノイズパワーをやや多めに見積もる「過減算（over-subtraction）」は、この周期グラムの分散に由来するゲイン \(H(f)\) のばらつきを抑える方向に働き、音声強調分野のスペクトル減算法でも知られた経験則です。しかし過大評価がさらに進む（\(+300\%\) 以降）と、今度はバイアス（信号成分そのものを誤って減衰させる副作用）が支配的になりMSEは単調に悪化します。実験では \(+200\%\) 付近でMSEが最小（\(0.076315\) ）となるU字型の挙動が確認でき、「真のノイズパワーを使うことが常に最善とは限らない」という、ノイズ統計量推定に伴う実務上の教訓が定量的に裏付けられました。
カルマンフィルタとの関係
| 特徴 | ウィーナーフィルタ | カルマンフィルタ |
|---|---|---|
| 前提 | 定常過程 | 非定常でも可 |
| 処理方式 | バッチ（周波数領域） | 逐次（時間領域） |
| 最適性 | MSE最小（定常） | MSE最小（線形ガウス） |
| 適用範囲 | 信号処理 | 状態推定・制御 |
ウィーナーフィルタは定常過程における最適フィルタであり、 カルマンフィルタ はその非定常・逐次処理への拡張と見なせます。両者は無関係な理論ではなく、定常状態のカルマンフィルタは因果的ウィーナーフィルタそのものです。これを、状態空間モデル \(s_k = a\,s_{k-1}+w_k\) （\(w_k\sim\mathcal{N}(0,q)\) ）、観測 \(z_k = s_k+v_k\) （\(v_k\sim\mathcal{N}(0,r)\) 、\(a=0.9,\ q=1.0,\ r=2.0\) ）で数値的に検証します。
まず scipy.linalg.solve_discrete_are で離散代数リカッチ方程式を解き、定常カルマンゲインを求めます。
from scipy.linalg import solve_discrete_are
import numpy as np
a, q, r = 0.9, 1.0, 2.0
P_pred = solve_discrete_are(np.array([[a]]), np.array([[1.0]]),
np.array([[q]]), np.array([[r]]))[0, 0]
K_inf = P_pred / (P_pred + r)
P_filt = (1 - K_inf) * P_pred
結果は \(P_{k|k-1}=1.757791\) （予測誤差分散）、定常カルマンゲイン \(K_\infty=0.467772\) 、更新後誤差分散 \(P_{k|k}=0.935545\) でした。
次に、この状態空間モデルが定める真の自己相関 \(R_{ss}(k)=\frac{q}{1-a^2}a^{|k|}\) 、\(R_{zz}(k)=R_{ss}(k)+r\,\delta_{k,0}\) 、\(R_{sz}(k)=R_{ss}(k)\) から、タップ数 \(M\) を増やしながら因果的FIRウィーナーフィルタ（式(6)の \(\mathbf{R}_{xx}\mathbf{w}=\mathbf{r}_{xd}\) ）を解き、理論MMSE \(=\sigma_s^2-\mathbf{w}^\top\mathbf{r}_{xd}\) を計算しました。
| \(M\) | MMSE(\(M\) ) | \(w_0\) | \(w_1/w_0\) |
|---|---|---|---|
| 1 | 1.449275 | 0.724638 | — |
| 2 | 1.041667 | 0.520833 | 0.600000 |
| 4 | 0.940979 | 0.470489 | 0.485863 |
| 8 | 0.935560 | 0.467780 | 0.479024 |
| 16 | 0.935545 | 0.467772 | 0.479005 |
| 30 | 0.935545 | 0.467772 | 0.479005 |
タップ数を増やすとMMSE(\(M\) )は単調に減少し、\(M=16\) 以降 \(0.935545\) に収束します。これは定常カルマンフィルタの更新後誤差分散 \(P_{k|k}=0.935545\) に小数点以下6桁まで一致します。さらに、収束後のFIR係数の先頭 \(w_0=0.467772\) はカルマンゲイン \(K_\infty=0.467772\) に、隣接タップ比 \(w_1/w_0=0.479005\) はカルマン再帰の等価な極 \((1-K_\infty)a=0.479005\) に、それぞれ完全に一致します。これは理論的にも当然で、定常カルマン再帰 \(\hat{s}_{k|k}=(1-K_\infty)a\,\hat{s}_{k-1|k-1}+K_\infty z_k\) を逐次展開すると
\[ \hat{s}_{k|k} = K_\infty \sum_{j=0}^{\infty} \bigl[(1-K_\infty)a\bigr]^j z_{k-j} \]という無限タップの因果的ウィーナーフィルタの形になり、その係数は等比数列 \(w_j = K_\infty\bigl[(1-K_\infty)a\bigr]^j\) に一致するからです。
最後に、\(N=5000\)
サンプルの独立な系列を30回生成し（np.random.default_rng(seed)、\(\text{seed}=0,\ldots,29\)
）、定常カルマンフィルタと \(M=40\)
のFIRウィーナーフィルタを同一データに適用してモンテカルロ平均MSEを比較しました。結果は両者とも \(0.934818 \pm 0.003926\)
（標準誤差）で完全に一致し、理論値 \(P_{k|k}=0.935545\)
とも1標準誤差以内で整合しました。以上より、「定常状態でカルマンゲインがウィーナーフィルタに収束する」という関係が、解析的一致・数値実験の両面で確認できました。
適応ウィーナーフィルタ：LMS・RLSによる逐次近似
ここまでの議論は、\(\mathbf{R}_{xx}\) と \(\mathbf{r}_{xd}\) （あるいは \(S_{yy}(f)\) 、\(S_{sy}(f)\) ）が既知であることを前提としてきました。しかし実際の環境では信号・ノイズの統計量は未知、あるいは非定常（時間とともに変化）であることが多く、式(6)を直接解くバッチ処理は使えません。 適応フィルタ（LMS/RLS）の理論とPython実装 、 LMS／NLMSアルゴリズムの理論とPython実装 、 RLS（逐次最小二乗法）アルゴリズムの理論とPython実装 で詳述したLMS・RLSは、観測サンプルから逐次的に式(6)の解に近づく適応アルゴリズムであり、本記事のウィーナー解を実データのみから（統計量を明示的に推定せずに）実現する手段と位置づけられます。
- LMSは瞬時勾配 \(\nabla J \approx -2e[n]\mathbf{x}[n]\) （真の期待値の代わりに単一サンプルの積を使う確率的近似）で係数を更新し、平均的にウィーナー解 \(\mathbf{w}^\ast=\mathbf{R}^{-1}\mathbf{p}\) に線形収束します。
- RLSは式(6)と同型の正規方程式 \(\mathbf{R}(n)\mathbf{w}(n)=\mathbf{r}(n)\) を、行列反転補題（Sherman-Morrison公式）によって毎時刻 \(O(M^2)\) で厳密に更新し、二次収束（入力の統計に依存せず最短時間で正規方程式の解に到達）します。実測では、雑音床の2倍のMSEに到達するまでのサンプル数がLMS 57サンプルに対しRLS 27サンプルと、約半分に短縮されることが確認されています。
言い換えると、本記事で解いた \(\mathbf{R}_{xx}\mathbf{w}=\mathbf{r}_{xd}\) はあくまで統計量が既知の場合の「正解」を与えるものであり、LMS・RLSはその正解を未知環境下でオンラインに近似するための実用的アルゴリズムです。定常性が崩れる環境（\(\mathbf{R}_{xx}\) や \(\mathbf{r}_{xd}\) が時間変化する場合）では、忘却係数を持つRLSやNLMSが有効であり、その意味で適応フィルタは「統計量未知・非定常」というウィーナーフィルタの前提の欠落を、カルマンフィルタとは異なる（モデルベースではなくデータ駆動の）アプローチで補完する手法群と言えます。
関連記事
- 適応フィルタの基礎理論とデジタル信号処理への応用 - 本記事のWiener-Hopf方程式（式(6)）を、MMSE最適化問題として勾配計算から導出した理論的な出発点です。悪条件な\(R\) ・LMSとの収束関係・非定常性という3つのエッジケースも数値検証しています。
- バターワースフィルタの設計と実装 - カットオフ周波数を手動で設定するバターワースとの対比で、ウィーナーの自動的な帯域設計が際立ちます。
- ローパスフィルタの設計と比較 - 移動平均・バターワース・チェビシェフを周波数応答で比較した古典的設計と、SNRから帯域を決めるウィーナー設計の対比が理解できます。
- チェビシェフフィルタの設計原理とPython実装 - 等リップルで急峻な遮断を実現する古典的IIRフィルタ。ウィーナーの「SNRに従って自動的に帯域を絞る」設計思想と対比できます。
- FIRフィルタとIIRフィルタの比較 - ウィーナーフィルタと古典的フィルタの設計思想の違いを理解できます。
- FFTの基礎：高速フーリエ変換の理論と実装 - ウィーナーフィルタの実装に不可欠なFFTの解説です。
- カルマンフィルタの理論とPython実装 - ウィーナーフィルタの非定常・逐次処理への拡張であるカルマンフィルタを解説しています。
- 時系列データの異常検知 - ウィーナーフィルタの応用先の一つである異常検知手法を解説しています。
- 指数移動平均（EMA）フィルタの周波数特性 - EMAは簡便なローパスフィルタですが、ウィーナーフィルタはSNRに基づく最適な帯域設計です。
- 移動平均フィルタの種類と比較 - 移動平均とウィーナーフィルタの設計思想の違いを理解できます。
- 適応フィルタ（LMS/RLS）の理論とPython実装 - ウィーナーフィルタの最適解を逐次的に近似するLMSアルゴリズムと、より高速に収束するRLSアルゴリズムを解説しています。
- 窓関数とパワースペクトル密度（PSD）の理論とPython実装 - ウィーナーフィルタの設計に必要な信号・ノイズのPSD推定を、Welch法と窓関数の選択を含めて詳しく解説しています。
- 相補フィルタ（Complementary Filter）の理論とPython実装 - 高周波・低周波で異なる信号源を融合する相補フィルタは、SNRの周波数依存性に基づくウィーナーフィルタの直感的な特殊ケースとして理解できます。
- LMS／NLMSアルゴリズムの理論とPython実装 - 本記事のWiener解 \(\mathbf{R}^{-1}\mathbf{p}\) を逐次的に近似する適応アルゴリズム。バッチ処理 vs 逐次処理の対比が理解できます。
- RLS（逐次最小二乗法）アルゴリズムの理論とPython実装 - 本記事の正規方程式を行列反転補題で逐次的に厳密解にする適応アルゴリズム。LMSより高速な二次収束、カルマンフィルタとの数値的な完全一致を確認しています。
- 短時間フーリエ変換（STFT）の理論とPython実装 - 時変な信号にウィーナーを適用する場合の前処理。フレームごとに PSD を推定して時変 Wiener gain を構成する応用が広がります。
- デジタルフィルタ設計指針ハブ：選定の3軸と特性比較 — 本記事のウィーナーフィルタを他のフィルタ群と比較しながら選定する際の指針として活用できます。
参考文献
- Haykin, S. (2014). Adaptive Filter Theory (5th ed.). Pearson. Chapters 2-3.
- Wiener, N. (1949). Extrapolation, Interpolation, and Smoothing of Stationary Time Series. MIT Press.
- Vaseghi, S. V. (2008). Advanced Digital Signal Processing and Noise Reduction (4th ed.). Wiley.
- Kailath, T. (1974). A View of Three Decades of Linear Filtering Theory. IEEE Transactions on Information Theory, 20(2), 146-181.