ラジオのチューニングを合わせたとき、背景にサーっというノイズが混じって聞こえることがあります。あるいは、古いレコードの音楽にはプチプチという針音がつきものです。このように、私たちが本当に聞きたい信号は、常に何らかの雑音に汚染されています。では、雑音に埋もれた信号をできるだけ正確に復元するには、どのようなフィルタを設計すればよい のでしょうか。
この問いに対する最も基本的かつ強力な答えの1つが ウィーナフィルタ（Wiener filter） です。1940年代にノーバート・ウィーナーが定式化したこのフィルタは、最小平均二乗誤差（Minimum Mean Square Error, MMSE） という基準のもとで「最適」であることが数学的に証明されています。つまり、信号と雑音の統計的性質（パワースペクトル）がわかっている限り、これ以上良い線形フィルタは存在しないのです。
ウィーナフィルタを理解すると、次のような分野で強力な理論的基盤を手に入れることができます。
- 通信工学: 受信信号からチャネルの歪みと雑音を同時に除去するチャネル等化器の設計原理がわかります
- 音声処理: マイクで収録した音声から環境雑音を除去するスペクトルサブトラクション法の理論的根拠を理解できます
- 画像処理: ぼけた画像やノイズの乗った画像を復元する逆フィルタリングの基盤となります
- 制御工学: カルマンフィルタの原型であり、最適推定理論への入口として位置づけられます
本記事の内容
- MMSE基準の定義と「なぜ二乗誤差なのか」の直感
- ウィーナー・ホップ方程式の導出（時間領域）
- 周波数領域でのウィーナフィルタの閉じた形の解
- 因果的ウィーナフィルタとスペクトル分解
- FIRウィーナフィルタと正規方程式（ウィーナー・ホップ方程式の離散版）
- LMSアルゴリズムとの関係
- 応用例（雑音除去・チャネル等化・予測）
- Pythonによるウィーナフィルタの実装と雑音除去デモ
- LMSフィルタとの比較実験
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
MMSE基準とは — 「最も正確な復元」を数学で定義する
なぜ二乗誤差を使うのか
フィルタの「良さ」を議論するには、まず「良さ」を定量化する基準が必要です。日常的な直感で考えてみましょう。射撃の練習をしていて、的の中心からどれだけずれたかを測りたいとします。右に5cm外れたショットと左に5cm外れたショットを単純に足し合わせると、偏差がプラスとマイナスで打ち消し合ってしまい「平均的には完璧だ」という誤った結論になります。
この問題を回避する最も自然な方法は、偏差を二乗してから平均をとる ことです。二乗すれば常に非負になるため、正と負が打ち消し合う心配がありません。さらに、二乗誤差には微分可能であるという数学的な利点があり、最適化計算が扱いやすくなります。
MMSE基準の数学的定義
信号処理の文脈では、次のような設定を考えます。私たちが本当に知りたい 所望信号（desired signal） を $d(n)$ とし、手元にある 観測信号 を $x(n)$ とします。観測信号は、所望信号に雑音 $v(n)$ が加わったものであることが多いです。
$$ x(n) = d(n) + v(n) $$
フィルタ $h(n)$ を $x(n)$ に適用して、所望信号の推定値 $\hat{d}(n)$ を得ます。
$$ \hat{d}(n) = \sum_{k=-\infty}^{\infty} h(k) \, x(n – k) = (h * x)(n) $$
ここで $*$ は畳み込み演算を表します。推定値と真の所望信号との誤差を $e(n)$ とおくと、
$$ e(n) = d(n) – \hat{d}(n) $$
です。MMSE基準とは、この誤差の二乗の期待値（統計的平均）を最小化することです。
$$ J = E\bigl[|e(n)|^2\bigr] = E\bigl[|d(n) – \hat{d}(n)|^2\bigr] $$
この $J$ を コスト関数 または 目的関数 と呼びます。$J$ を最小にするフィルタ $h(n)$ が、MMSE基準のもとで 最適フィルタ すなわちウィーナフィルタです。
ここで重要な前提があります。上の定式化では、$d(n)$、$x(n)$、$v(n)$ がすべて 広義定常（wide-sense stationary, WSS） であることを仮定します。つまり、平均と自己相関関数が時刻に依存しないということです。この仮定のおかげで、最適フィルタはインパルス応答 $h(n)$ が時不変なものとして求められます。
では、この最適化問題を実際に解くとどうなるのか、次のセクションで導出していきましょう。
ウィーナー・ホップ方程式の導出
直交性原理からの出発
コスト関数 $J = E[|e(n)|^2]$ を最小化するために、最も本質的な道具は 直交性原理（orthogonality principle） です。
直交性原理とは、直感的に言えば「最適推定を行ったとき、推定誤差はフィルタへの入力と無相関になる」という主張です。幾何学的には、ベクトル空間への射影を考えるとわかりやすいでしょう。3次元空間のベクトルを2次元平面に射影したとき、射影誤差（残差ベクトル）は平面に直交します。この「残差が直交する」という性質が、射影が「最も近い点」であることの本質なのです。MMSE推定も全く同じ構造を持っています。
数学的に述べると、最適フィルタのもとで推定誤差 $e(n)$ は、フィルタ入力 $x(n-m)$（任意の $m$）と直交します。
$$ E\bigl[e(n) \, x^*(n – m)\bigr] = 0, \quad \forall m $$
ここで $^*$ は複素共役を表します（実数信号の場合は不要です）。
直交性原理の証明
なぜ直交性原理が成り立つのかを簡潔に示しましょう。コスト関数をフィルタ係数 $h(k)$ で微分して、それをゼロとおきます。
まず、誤差を明示的に書き下します。
$$ e(n) = d(n) – \sum_{k=-\infty}^{\infty} h(k) \, x(n-k) $$
コスト関数 $J$ を $h(m)$ で偏微分します。$J = E[|e(n)|^2] = E[e(n) \, e^*(n)]$ であることに注意すると、
$$ \frac{\partial J}{\partial h(m)} = E\left[\frac{\partial e(n)}{\partial h(m)} \, e^*(n)\right] + E\left[e(n) \, \frac{\partial e^*(n)}{\partial h(m)}\right] $$
実数信号の場合を考えると（複素の場合はウィルティンガー微分を使いますが本質は同じです）、
$$ \frac{\partial e(n)}{\partial h(m)} = -x(n – m) $$
なので、最適条件 $\frac{\partial J}{\partial h(m)} = 0$ は次のようになります。
$$ E\bigl[e(n) \, x(n – m)\bigr] = 0, \quad \forall m $$
これが直交性原理です。コスト関数の極値条件が、推定誤差とフィルタ入力の直交関係として表現されるわけです。
ウィーナー・ホップ方程式の導出
直交性原理の式に $e(n) = d(n) – \hat{d}(n)$ を代入して展開しましょう。
$$ E\bigl[(d(n) – \hat{d}(n)) \, x(n – m)\bigr] = 0 $$
$\hat{d}(n) = \sum_k h(k) \, x(n-k)$ を代入すると、
$$ E\bigl[d(n) \, x(n-m)\bigr] – \sum_{k=-\infty}^{\infty} h(k) \, E\bigl[x(n-k) \, x(n-m)\bigr] = 0 $$
ここで、広義定常過程の相互相関関数と自己相関関数を導入します。所望信号 $d(n)$ と観測信号 $x(n)$ の相互相関関数を
$$ R_{dx}(m) = E\bigl[d(n) \, x(n – m)\bigr] $$
観測信号の自己相関関数を
$$ R_{xx}(l) = E\bigl[x(n) \, x(n – l)\bigr] $$
と定義すると、上の式は次のように書けます。$k$ を $l$ に置き換え、遅延の差 $l – m$ に注目すると、
$$ R_{dx}(m) = \sum_{k=-\infty}^{\infty} h(k) \, R_{xx}(m – k), \quad \forall m $$
これが ウィーナー・ホップ方程式（Wiener-Hopf equation） です。右辺は $h$ と $R_{xx}$ の畳み込みになっていることに注意してください。つまり、
$$ R_{dx}(m) = (h * R_{xx})(m) $$
と書くこともできます。
ウィーナー・ホップ方程式は「最適フィルタのインパルス応答 $h(k)$ が満たすべき条件」を与える方程式です。この方程式が解ければ、ウィーナフィルタが得られるわけです。
時間領域でこの畳み込み方程式を直接解くのは一般に困難です。しかし、畳み込みはフーリエ変換で積に変わるという性質を利用すれば、周波数領域で非常にエレガントな解が得られます。次のセクションで見ていきましょう。
周波数領域でのウィーナフィルタ
畳み込み定理による変換
ウィーナー・ホップ方程式の両辺をフーリエ変換します。離散時間系であれば離散時間フーリエ変換（DTFT）を、連続時間系であればフーリエ変換を使いますが、ここでは連続周波数 $f$（あるいは角周波数 $\omega$）での議論を行います。
畳み込みの性質から、時間領域での畳み込みは周波数領域での積に対応します。
$$ S_{dx}(f) = H(f) \, S_{xx}(f) $$
ここで、$S_{dx}(f)$ は $d(n)$ と $x(n)$ のクロスパワースペクトル密度、$S_{xx}(f)$ は $x(n)$ のパワースペクトル密度、$H(f)$ はウィーナフィルタの周波数応答です。これらはそれぞれ、$R_{dx}(m)$、$R_{xx}(m)$、$h(k)$ のフーリエ変換です。
$H(f)$ について解くと、
$$ \boxed{H_{\mathrm{opt}}(f) = \frac{S_{dx}(f)}{S_{xx}(f)}} $$
これが 非因果的（non-causal）ウィーナフィルタ の周波数領域での解です。「非因果的」と呼ぶのは、フィルタのインパルス応答 $h(k)$ が $k < 0$ の部分も含みうる、すなわち未来の入力も使える理想的な状況を想定しているからです。
加法性雑音モデルでの具体的な形
特に重要な特殊ケースとして、$x(n) = d(n) + v(n)$（信号 + 雑音）であり、信号 $d(n)$ と雑音 $v(n)$ が互いに無相関であるという仮定を追加しましょう。
この仮定のもとでは、クロススペクトルとパワースペクトルが次のように分解されます。
信号と雑音が無相関なので、$E[d(n) \, v(n-m)] = 0$ です。まず、$S_{xx}(f)$ を展開します。$x = d + v$ なので、
$$ R_{xx}(m) = E[(d(n) + v(n))(d(n-m) + v(n-m))] $$
信号と雑音の無相関性を使って展開すると、
$$ R_{xx}(m) = R_{dd}(m) + R_{vv}(m) $$
したがって、フーリエ変換をとると、
$$ S_{xx}(f) = S_{dd}(f) + S_{vv}(f) $$
次に、$S_{dx}(f)$ を計算します。
$$ R_{dx}(m) = E[d(n) \, x(n-m)] = E[d(n)(d(n-m) + v(n-m))] = R_{dd}(m) + 0 = R_{dd}(m) $$
よって、$S_{dx}(f) = S_{dd}(f)$ です。
これらをウィーナフィルタの式に代入すると、
$$ \boxed{H_{\mathrm{opt}}(f) = \frac{S_{dd}(f)}{S_{dd}(f) + S_{vv}(f)}} $$
この式の物理的意味は極めて明快です。
- 信号のパワーが雑音に比べて大きい周波数帯域 では $S_{dd}(f) \gg S_{vv}(f)$ なので $H_{\mathrm{opt}}(f) \approx 1$ となり、信号をそのまま通します
- 雑音のパワーが信号に比べて大きい周波数帯域 では $S_{dd}(f) \ll S_{vv}(f)$ なので $H_{\mathrm{opt}}(f) \approx 0$ となり、その帯域を抑圧します
つまり、ウィーナフィルタは 周波数ごとにSN比（信号対雑音比）に応じた重みづけを行う フィルタなのです。SN比が高い帯域は通し、低い帯域は遮断する — この振る舞いは直感にも完全に合致します。
ところで、この非因果的ウィーナフィルタは未来の入力も使えるという理想的な仮定に基づいていました。リアルタイム処理のように、現在と過去のデータだけでフィルタリングしなければならない場合はどうすればよいでしょうか。次のセクションでは、この制約を課した 因果的ウィーナフィルタ を扱います。
因果的ウィーナフィルタとスペクトル分解
因果性制約とは
実時間システムでは、フィルタは未来の入力を利用できません。数学的には、インパルス応答 $h(k) = 0 \; (k < 0)$ という制約を課すことになります。この制約のもとでの最適フィルタを 因果的ウィーナフィルタ と呼びます。
非因果的な場合のように単純に $S_{dx}/S_{xx}$ とするだけでは因果性が保証されません。なぜなら、$S_{dx}(f)/S_{xx}(f)$ を逆フーリエ変換しても、一般に $k < 0$ の成分を持つインパルス応答が得られるからです。
スペクトル分解法
因果的ウィーナフィルタを求めるための古典的手法が スペクトル分解（spectral factorization） です。
パワースペクトル密度 $S_{xx}(f)$ は非負実数値をとるので、次のようにスペクトル分解できます。
$$ S_{xx}(f) = \Phi_x^+(f) \, \Phi_x^-(f) $$
ここで $\Phi_x^+(f)$ は「因果的な部分」（逆フーリエ変換が $k \geq 0$ でのみ非ゼロ）、$\Phi_x^-(f)$ は「反因果的な部分」（逆フーリエ変換が $k \leq 0$ でのみ非ゼロ）です。$\Phi_x^-(f) = [\Phi_x^+(f)]^*$ の関係があります。
この分解を使うと、ウィーナー・ホップ方程式は次のように変形できます。
$$ H_{\mathrm{opt}}(f) = \frac{S_{dx}(f)}{S_{xx}(f)} = \frac{S_{dx}(f)}{\Phi_x^+(f) \, \Phi_x^-(f)} $$
因果性制約を満たす解を得るために、右辺を以下の手順で処理します。
ステップ1: $S_{dx}(f) / \Phi_x^-(f)$ を計算し、これを逆フーリエ変換します。
ステップ2: 逆フーリエ変換の結果から、$k \geq 0$ の部分だけを取り出します。この操作を $[\cdot]_+$ と書きます。
ステップ3: 再びフーリエ変換して $\Phi_x^+(f)$ で割ると、因果的ウィーナフィルタが得られます。
$$ \boxed{H_{\mathrm{causal}}(f) = \frac{1}{\Phi_x^+(f)} \left[\frac{S_{dx}(f)}{\Phi_x^-(f)}\right]_+} $$
$[\cdot]_+$ の操作が因果性を保証する鍵です。非因果的な成分を「切り捨てる」ことで、過去と現在の入力だけで実現可能なフィルタを得るわけです。ただし、この切り捨てによって非因果的ウィーナフィルタよりも性能（MMSE）は悪化します。これは、未来の情報を使えないことの代償です。
スペクトル分解法は理論的にはエレガントですが、実用上は解析的にスペクトル分解を実行できるケースが限られます。そこで、より実用的なアプローチとして、有限長のフィルタ係数を直接求める FIRウィーナフィルタ が広く使われています。次にこの方法を見ていきましょう。
FIRウィーナフィルタと正規方程式
有限長フィルタの定式化
実際のディジタル信号処理では、無限長のインパルス応答を実装することはできません。そこで、フィルタの長さを $M$ に制限した FIRウィーナフィルタ を考えます。
フィルタの出力は次のように表せます。
$$ \hat{d}(n) = \sum_{k=0}^{M-1} h(k) \, x(n-k) = \bm{h}^T \bm{x}(n) $$
ここで、フィルタ係数ベクトルと入力ベクトルを次のように定義しました。
$$ \bm{h} = \begin{pmatrix} h(0) \\ h(1) \\ \vdots \\ h(M-1) \end{pmatrix}, \quad \bm{x}(n) = \begin{pmatrix} x(n) \\ x(n-1) \\ \vdots \\ x(n-M+1) \end{pmatrix} $$
コスト関数はこれまでと同じく $J = E[|e(n)|^2]$ です。
正規方程式の導出
コスト関数を展開してみましょう。
$$ J = E\bigl[|d(n) – \bm{h}^T \bm{x}(n)|^2\bigr] $$
$\bm{h}$ について微分してゼロとおきます。ここで、$\nabla_{\bm{h}} J = 0$ を計算するために、$J$ をもう少し展開します。
$$ J = E[d^2(n)] – 2\bm{h}^T E[\bm{x}(n) \, d(n)] + \bm{h}^T E[\bm{x}(n) \bm{x}^T(n)] \bm{h} $$
ここで、観測信号の 自己相関行列 と 相互相関ベクトル を次のように定義します。
$$ \bm{R}_{xx} = E[\bm{x}(n) \bm{x}^T(n)], \quad \bm{r}_{dx} = E[\bm{x}(n) \, d(n)] $$
自己相関行列 $\bm{R}_{xx}$ は $M \times M$ の対称正定値行列で、$(i,j)$ 成分は $R_{xx}(i-j)$ です。つまり テプリッツ行列 になっています。このテプリッツ構造は、広義定常過程の自己相関が遅延の差にのみ依存するという性質の直接的な帰結です。
相互相関ベクトル $\bm{r}_{dx}$ の第 $i$ 成分は $R_{dx}(i)$ です。
コスト関数を $\bm{h}$ で微分してゼロとおくと、
$$ \nabla_{\bm{h}} J = -2\bm{r}_{dx} + 2\bm{R}_{xx} \bm{h} = \bm{0} $$
整理すると、
$$ \boxed{\bm{R}_{xx} \bm{h}_{\mathrm{opt}} = \bm{r}_{dx}} $$
これが 正規方程式（normal equation） であり、ウィーナー・ホップ方程式の離散・有限長版です。英語では Wiener-Hopf equation in matrix form とも呼ばれます。
正規方程式の解
$\bm{R}_{xx}$ が正定値であれば逆行列が存在し、最適フィルタ係数は一意に定まります。
$$ \bm{h}_{\mathrm{opt}} = \bm{R}_{xx}^{-1} \bm{r}_{dx} $$
最小平均二乗誤差の値
最適フィルタを使ったときの最小コスト $J_{\min}$ も重要な量です。$\bm{h}_{\mathrm{opt}} = \bm{R}_{xx}^{-1} \bm{r}_{dx}$ をコスト関数に代入すると、
$$ J_{\min} = E[d^2(n)] – \bm{r}_{dx}^T \bm{R}_{xx}^{-1} \bm{r}_{dx} $$
$E[d^2(n)] = R_{dd}(0)$ は所望信号の電力です。第2項 $\bm{r}_{dx}^T \bm{R}_{xx}^{-1} \bm{r}_{dx}$ は、フィルタが「説明できた」誤差の量に相当します。$J_{\min} \geq 0$ であることは自明ですが、逆に $J_{\min} = 0$ となるのは所望信号が観測信号の線形結合で完全に表現できる場合のみです。
テプリッツ構造の活用 — レビンソン・ダービンアルゴリズム
$\bm{R}_{xx}$ がテプリッツ行列であるという特殊構造を活用すると、一般的なガウスの消去法（$O(M^3)$）よりも高速に正規方程式を解くことができます。レビンソン・ダービンアルゴリズム（Levinson-Durbin algorithm） は $O(M^2)$ の計算量で解を求めるアルゴリズムであり、実用上よく使われます。
正規方程式はウィーナフィルタの「バッチ解」を与えるものです。つまり、$\bm{R}_{xx}$ と $\bm{r}_{dx}$ を事前に知っている（あるいはデータから推定できる）場合に使えます。しかし、信号の統計的性質が未知あるいは時変である場合はどうすればよいでしょうか。この課題に対処するのが、適応フィルタリングのアプローチです。次のセクションでLMSアルゴリズムとの関係を見ていきましょう。
LMSアルゴリズムとの関係
勾配降下法としてのLMS
LMS（Least Mean Squares）アルゴリズムは、ウィーナフィルタの最適解に 適応的に 近づいていくアルゴリズムです。LMSアルゴリズムの記事で詳しく解説しましたが、ここではウィーナフィルタとの関係を明確にしておきます。
コスト関数 $J = E[|e(n)|^2]$ は $\bm{h}$ の二次関数であり、その等高線は楕円になります（$\bm{R}_{xx}$ が正定値のため）。この二次関数の勾配は、
$$ \nabla_{\bm{h}} J = -2\bm{r}_{dx} + 2\bm{R}_{xx} \bm{h} $$
です。最急降下法（steepest descent）はこの勾配に比例してフィルタ係数を更新します。
$$ \bm{h}(n+1) = \bm{h}(n) – \mu \nabla_{\bm{h}} J = \bm{h}(n) + 2\mu \bigl(\bm{r}_{dx} – \bm{R}_{xx} \bm{h}(n)\bigr) $$
ここで $\mu$ はステップサイズ（学習率）です。この更新式は、$\bm{R}_{xx}$ と $\bm{r}_{dx}$ がわかっている場合に使えますが、実際にはこれらの統計量は未知であることがほとんどです。
確率的近似 — 瞬時推定値の利用
LMSアルゴリズムのアイデアは、真の勾配 $\nabla_{\bm{h}} J$ を 瞬時推定値 で置き換えることです。具体的には、期待値を1サンプルでの瞬時値に置き換えます。
$$ E[\bm{x}(n) \, d(n)] \approx \bm{x}(n) \, d(n) $$
$$ E[\bm{x}(n) \bm{x}^T(n)] \bm{h}(n) \approx \bm{x}(n) \bm{x}^T(n) \bm{h}(n) = \bm{x}(n) \hat{d}(n) $$
これを代入すると、瞬時勾配推定値は
$$ \hat{\nabla} J = -2\bm{x}(n) \bigl(d(n) – \hat{d}(n)\bigr) = -2\bm{x}(n) \, e(n) $$
となり、LMSの更新式が得られます。
$$ \boxed{\bm{h}(n+1) = \bm{h}(n) + 2\mu \, e(n) \, \bm{x}(n)} $$
ウィーナフィルタとLMSの比較
両者の関係を整理すると、以下のようになります。
| 特性 | ウィーナフィルタ | LMSアルゴリズム |
|---|---|---|
| 解法 | 閉じた形の解（正規方程式） | 反復的更新 |
| 必要な情報 | $\bm{R}_{xx}$, $\bm{r}_{dx}$（統計量が既知） | 入力 $\bm{x}(n)$ と所望出力 $d(n)$ のサンプル列 |
| 計算量 | $O(M^2)$（レビンソン・ダービン）～$O(M^3)$ | 1サンプルあたり $O(M)$ |
| 時変環境への追従 | 不可（定常性仮定） | 可能（適応的に追従） |
| 収束後の定常状態 | 最適解（MMSE解） | 最適解の周辺で揺らぐ（ミスアジャストメント） |
LMSは「ウィーナフィルタの近似解を逐次的に求めるアルゴリズム」と位置づけることができます。統計量が既知で定常環境であればウィーナフィルタの方が直接的で最適ですが、未知かつ時変の環境ではLMSのような適応アルゴリズムが威力を発揮します。
ウィーナフィルタの理論がわかったところで、具体的にどのような場面で活用されるのか、代表的な応用を見ていきましょう。
ウィーナフィルタの応用
雑音除去（ノイズリダクション）
最も直接的な応用は、信号に加法性雑音が重畳した状況での雑音除去です。観測信号 $x(n) = d(n) + v(n)$ から所望信号 $d(n)$ を推定する問題で、先ほど導出した
$$ H_{\mathrm{opt}}(f) = \frac{S_{dd}(f)}{S_{dd}(f) + S_{vv}(f)} $$
をそのまま適用します。実用上は、$S_{dd}(f)$ と $S_{vv}(f)$ をデータから推定する必要があります。音声処理ではVAD（Voice Activity Detection）を使って無音区間を検出し、無音区間のパワースペクトルを $S_{vv}(f)$ の推定値として使うことが一般的です。
チャネル等化（イコライゼーション）
ディジタル通信では、送信信号 $s(n)$ がチャネル（伝送路）を通過する際に歪みと雑音が加わります。チャネルのインパルス応答を $c(n)$ とすると、受信信号は
$$ x(n) = (c * s)(n) + v(n) = \sum_k c(k) \, s(n-k) + v(n) $$
となります。チャネル等化器は、受信信号 $x(n)$ から元の送信信号 $s(n)$ を（遅延を除いて）復元するフィルタです。所望信号を $d(n) = s(n – \Delta)$（$\Delta$ は許容される遅延）とおけば、ウィーナフィルタの枠組みで等化器を設計できます。
チャネルの周波数応答を $C(f)$ とすると、ウィーナ等化器の周波数応答は
$$ H_{\mathrm{eq}}(f) = \frac{C^*(f)}{|C(f)|^2 + S_{vv}(f)/S_{ss}(f)} $$
と書けます。これは雑音がない極限（$S_{vv} \to 0$）でチャネルの逆フィルタ $1/C(f)$ に一致し、雑音が大きい帯域では振幅を抑えるという、直感的にも納得できる振る舞いです。
信号予測（ウィーナ予測器）
ウィーナフィルタの枠組みは、所望信号を「未来の自分自身」に設定することで、信号予測にも使えます。つまり、$d(n) = x(n + \Delta)$（$\Delta > 0$）とおけば、$\Delta$ ステップ先の予測値を出力するフィルタが設計できます。
この場合、相互相関が $R_{dx}(m) = R_{xx}(m + \Delta)$ となるだけで、正規方程式の構造は全く同じです。ウィーナ予測器は、線形予測符号化（LPC）や株価予測、気象予測など、さまざまな分野で応用の基盤となっています。
これらの応用を踏まえた上で、実際にPythonでウィーナフィルタを実装し、雑音除去の効果を目で確認してみましょう。
Pythonによるウィーナフィルタの実装
周波数領域ウィーナフィルタによる雑音除去
まず、最も基本的な非因果的ウィーナフィルタを周波数領域で実装します。正弦波の合成信号にホワイトノイズを加え、ウィーナフィルタで復元するデモを行います。
import numpy as np
import matplotlib.pyplot as plt
# --- 信号生成 ---
np.random.seed(42)
N = 1024 # サンプル数
fs = 256 # サンプリング周波数 [Hz]
t = np.arange(N) / fs
# 所望信号: 2つの正弦波の合成
d = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
# 雑音: ホワイトガウスノイズ
noise_power = 1.0
v = np.sqrt(noise_power) * np.random.randn(N)
# 観測信号
x = d + v
# --- 周波数領域ウィーナフィルタ ---
# パワースペクトル密度の推定（ピリオドグラム）
D = np.fft.fft(d)
X = np.fft.fft(x)
V = np.fft.fft(v)
# 真のパワースペクトル密度（理想的な場合）
Sdd = np.abs(D) ** 2 / N
Svv = np.abs(V) ** 2 / N
# ウィーナフィルタの周波数応答
H_wiener = Sdd / (Sdd + Svv)
# フィルタリング
D_hat = H_wiener * X
d_hat = np.real(np.fft.ifft(D_hat))
# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
axes[0].plot(t, d, color='#00d4ff', linewidth=0.8)
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Desired Signal d(n)')
axes[0].set_xlim([0, 2])
axes[1].plot(t, x, color='#ff6b6b', linewidth=0.5, alpha=0.8)
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Observed Signal x(n) = d(n) + v(n)')
axes[1].set_xlim([0, 2])
axes[2].plot(t, d, color='#00d4ff', linewidth=0.8, alpha=0.5, label='Original d(n)')
axes[2].plot(t, d_hat, color='#ffd93d', linewidth=0.8, label='Wiener estimate')
axes[2].set_ylabel('Amplitude')
axes[2].set_xlabel('Time [s]')
axes[2].set_title('Wiener Filter Output vs Original')
axes[2].set_xlim([0, 2])
axes[2].legend()
plt.tight_layout()
plt.savefig('wiener_filter_denoising.png', dpi=150, bbox_inches='tight')
plt.show()
上のコードでは、5 Hzと20 Hzの正弦波を合成した所望信号にホワイトノイズを加え、周波数領域でウィーナフィルタを適用しています。上段のグラフが元の所望信号、中段がノイズに埋もれた観測信号、下段がウィーナフィルタによる復元結果です。
結果を見ると、2つの重要なことが読み取れます。
- ウィーナフィルタの出力は元の信号をかなり正確に復元しています — ノイズに完全に埋もれていた正弦波が明瞭に再現されており、フィルタの雑音除去効果が視覚的に確認できます
- 高周波成分（20 Hz）の振幅がわずかに減衰しています — これは20 Hz付近でのSN比が5 Hz付近よりも低い（信号振幅が半分）ため、ウィーナフィルタがその帯域でやや強く抑圧したことを反映しています
次に、ウィーナフィルタの周波数応答を直接確認してみましょう。
import numpy as np
import matplotlib.pyplot as plt
# 信号パラメータ（前のコードと同じ）
np.random.seed(42)
N = 1024
fs = 256
t = np.arange(N) / fs
d = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise_power = 1.0
v = np.sqrt(noise_power) * np.random.randn(N)
x = d + v
# パワースペクトル密度
D = np.fft.fft(d)
V = np.fft.fft(v)
Sdd = np.abs(D) ** 2 / N
Svv = np.abs(V) ** 2 / N
# ウィーナフィルタの周波数応答
H_wiener = Sdd / (Sdd + Svv)
# 周波数軸
freqs = np.fft.fftfreq(N, d=1/fs)
mask = freqs >= 0 # 正の周波数のみ表示
fig, axes = plt.subplots(2, 1, figsize=(10, 6))
# パワースペクトル密度の比較
axes[0].semilogy(freqs[mask], Sdd[mask], color='#00d4ff', label='$S_{dd}(f)$ (signal)')
axes[0].semilogy(freqs[mask], Svv[mask], color='#ff6b6b', alpha=0.6, label='$S_{vv}(f)$ (noise)')
axes[0].set_ylabel('Power Spectral Density')
axes[0].set_title('Power Spectra of Signal and Noise')
axes[0].legend()
axes[0].set_xlim([0, fs/2])
# ウィーナフィルタの周波数応答
axes[1].plot(freqs[mask], H_wiener[mask], color='#ffd93d', linewidth=1.0)
axes[1].set_ylabel('$|H_{opt}(f)|$')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_title('Wiener Filter Frequency Response')
axes[1].set_xlim([0, fs/2])
axes[1].set_ylim([-0.05, 1.05])
plt.tight_layout()
plt.savefig('wiener_filter_frequency_response.png', dpi=150, bbox_inches='tight')
plt.show()
上段のグラフはパワースペクトル密度の比較です。5 Hzと20 Hzに信号成分のピークが立ち、それ以外の帯域では雑音が支配的であることがわかります。下段のウィーナフィルタの周波数応答を見ると、まさに信号のピーク周波数付近でゲインが1に近く、雑音が支配的な帯域ではゲインが0に近い値を取っています。これは先ほど理論的に予想した通りの振る舞いであり、「周波数ごとのSN比に応じた重みづけ」が実現されていることが視覚的に確認できます。
FIRウィーナフィルタの実装
次に、正規方程式を直接解くFIRウィーナフィルタを実装します。
import numpy as np
import matplotlib.pyplot as plt
def fir_wiener_filter(x, d, M):
"""
FIRウィーナフィルタの係数を正規方程式から求める
Parameters
----------
x : ndarray - 観測信号
d : ndarray - 所望信号
M : int - フィルタ長
Returns
-------
h_opt : ndarray - 最適フィルタ係数
"""
N = len(x)
# 自己相関行列 R_xx の推定
R_xx = np.zeros((M, M))
for i in range(M):
for j in range(M):
# サンプル平均で推定
total = 0.0
count = 0
for n in range(max(i, j), N):
total += x[n - i] * x[n - j]
count += 1
R_xx[i, j] = total / count
# 相互相関ベクトル r_dx の推定
r_dx = np.zeros(M)
for i in range(M):
total = 0.0
count = 0
for n in range(i, N):
total += d[n] * x[n - i]
count += 1
r_dx[i] = total / count
# 正規方程式を解く: R_xx @ h_opt = r_dx
h_opt = np.linalg.solve(R_xx, r_dx)
return h_opt
def apply_fir_filter(x, h):
"""FIRフィルタを適用する"""
M = len(h)
N = len(x)
y = np.zeros(N)
for n in range(M - 1, N):
y[n] = np.dot(h, x[n - M + 1:n + 1][::-1])
return y
# --- 信号生成 ---
np.random.seed(42)
N = 1024
fs = 256
t = np.arange(N) / fs
d = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise_power = 1.0
v = np.sqrt(noise_power) * np.random.randn(N)
x = d + v
# --- 異なるフィルタ長で実験 ---
filter_lengths = [8, 32, 128]
fig, axes = plt.subplots(len(filter_lengths) + 1, 1, figsize=(10, 10), sharex=True)
axes[0].plot(t, x, color='#ff6b6b', linewidth=0.5, alpha=0.6, label='Observed x(n)')
axes[0].plot(t, d, color='#00d4ff', linewidth=0.8, label='Desired d(n)')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Observed and Desired Signals')
axes[0].legend()
axes[0].set_xlim([0, 2])
for idx, M in enumerate(filter_lengths):
h_opt = fir_wiener_filter(x, d, M)
d_hat = apply_fir_filter(x, h_opt)
# MSE計算（フィルタ遅延分を除く）
mse = np.mean((d[M:] - d_hat[M:]) ** 2)
axes[idx + 1].plot(t, d, color='#00d4ff', linewidth=0.8, alpha=0.4, label='Original')
axes[idx + 1].plot(t, d_hat, color='#ffd93d', linewidth=0.8, label=f'FIR Wiener (M={M})')
axes[idx + 1].set_ylabel('Amplitude')
axes[idx + 1].set_title(f'FIR Wiener Filter (M={M}), MSE={mse:.4f}')
axes[idx + 1].legend()
axes[idx + 1].set_xlim([0, 2])
axes[-1].set_xlabel('Time [s]')
plt.tight_layout()
plt.savefig('fir_wiener_filter_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
このコードでは、フィルタ長 $M$ を8、32、128の3通りで実験しています。結果から以下のことが読み取れます。
- フィルタ長が短い（$M = 8$）場合は復元精度が低く、雑音が残っています — 8タップでは信号の周波数構造を十分に捉えられないためです。信号の最低周波数である5 Hzの1周期は $f_s / 5 = 51.2$ サンプルであり、8タップではこの周期成分を正確に表現できません
- フィルタ長を増やすと（$M = 32, 128$）MSEが減少し、復元が改善されます — より多くの遅延タップを使うことで、信号の自己相関構造をより精密に反映できるようになります
- ただし、フィルタ長を増やしすぎると計算量が増加し、有限データ長のもとでは相関関数の推定精度が落ちるトレードオフがあります — 実用ではデータ長やSN比を考慮して適切なフィルタ長を選ぶ必要があります
LMSフィルタとの比較
最後に、ウィーナフィルタ（正規方程式の解）とLMSアルゴリズムを同じ条件で比較します。ウィーナフィルタが「一発で最適解を得る」のに対し、LMSは反復的に近づいていく様子を確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
def lms_filter(x, d, M, mu):
"""
LMSアルゴリズムによる適応フィルタ
Parameters
----------
x : ndarray - 観測信号
d : ndarray - 所望信号
M : int - フィルタ長
mu : float - ステップサイズ
Returns
-------
y : ndarray - フィルタ出力
e : ndarray - 誤差信号
h_history : ndarray - 係数の履歴
"""
N = len(x)
h = np.zeros(M)
y = np.zeros(N)
e = np.zeros(N)
h_history = np.zeros((N, M))
for n in range(M - 1, N):
x_vec = x[n - M + 1:n + 1][::-1]
y[n] = np.dot(h, x_vec)
e[n] = d[n] - y[n]
h = h + 2 * mu * e[n] * x_vec
h_history[n] = h.copy()
return y, e, h_history
# --- 信号生成 ---
np.random.seed(42)
N = 4096 # LMSの収束を見るため長めのデータ
fs = 256
t = np.arange(N) / fs
d = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise_power = 1.0
v = np.sqrt(noise_power) * np.random.randn(N)
x = d + v
# --- フィルタ設計 ---
M = 64 # フィルタ長
# ウィーナフィルタ（正規方程式）
R_xx = np.zeros((M, M))
for i in range(M):
for j in range(M):
shifts = np.arange(max(i, j), N)
R_xx[i, j] = np.mean(x[shifts - i] * x[shifts - j])
r_dx = np.zeros(M)
for i in range(M):
shifts = np.arange(i, N)
r_dx[i] = np.mean(d[shifts] * x[shifts - i])
h_wiener = np.linalg.solve(R_xx, r_dx)
d_wiener = np.convolve(x, h_wiener, mode='full')[:N]
# LMSフィルタ
mu = 0.0005 # ステップサイズ
d_lms, e_lms, h_lms_history = lms_filter(x, d, M, mu)
# --- 学習曲線（MSEの推移） ---
# ウィーナフィルタのMSE（定数）
wiener_mse = np.mean((d[M:] - d_wiener[M:]) ** 2)
# LMSの学習曲線（移動平均で平滑化）
window = 64
e_lms_sq = e_lms ** 2
lms_learning_curve = np.convolve(e_lms_sq, np.ones(window)/window, mode='valid')
fig, axes = plt.subplots(3, 1, figsize=(10, 10))
# 学習曲線
axes[0].plot(lms_learning_curve, color='#ff6b6b', linewidth=0.8, label='LMS MSE (smoothed)')
axes[0].axhline(y=wiener_mse, color='#00d4ff', linestyle='--', linewidth=1.5,
label=f'Wiener MSE = {wiener_mse:.4f}')
axes[0].set_xlabel('Sample index')
axes[0].set_ylabel('MSE')
axes[0].set_title('Learning Curve: LMS vs Wiener Filter')
axes[0].legend()
axes[0].set_yscale('log')
# 収束後の比較（末尾2秒間）
t_end = t[-512:]
d_end = d[-512:]
d_wiener_end = d_wiener[-512:]
d_lms_end = d_lms[-512:]
axes[1].plot(t_end, d_end, color='#00d4ff', linewidth=0.8, alpha=0.5, label='Original d(n)')
axes[1].plot(t_end, d_wiener_end, color='#ffd93d', linewidth=0.8, label='Wiener')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Wiener Filter Output (last 2s)')
axes[1].legend()
axes[2].plot(t_end, d_end, color='#00d4ff', linewidth=0.8, alpha=0.5, label='Original d(n)')
axes[2].plot(t_end, d_lms_end, color='#ff6b6b', linewidth=0.8, label='LMS')
axes[2].set_ylabel('Amplitude')
axes[2].set_xlabel('Time [s]')
axes[2].set_title('LMS Filter Output (last 2s, after convergence)')
axes[2].legend()
plt.tight_layout()
plt.savefig('wiener_vs_lms_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
# 定量比較
lms_mse_final = np.mean(e_lms[-512:] ** 2)
print(f'Wiener filter MSE: {wiener_mse:.6f}')
print(f'LMS filter MSE (final 512 samples): {lms_mse_final:.6f}')
この比較実験から、以下の点が確認できます。
- 学習曲線（上段）を見ると、LMSのMSEは初期の大きな値からウィーナフィルタのMSE（破線）に向かって徐々に収束しています — LMSは反復を重ねるごとにウィーナフィルタの最適解に近づいていく様子が明確です。ただし、完全に破線の値に一致するのではなく、わずかに上をふらつきます。これが ミスアジャストメント（misadjustment） と呼ばれる現象で、瞬時勾配推定のノイズに起因します
- ウィーナフィルタの出力（中段）は最初から安定して高品質な復元を行っています — 統計量が既知である（あるいは十分なデータから正確に推定できる）場合、正規方程式を解くことで即座に最適フィルタが得られるというメリットが視覚的に明らかです
- LMSの出力（下段）は収束後にはウィーナフィルタに匹敵する品質を示しますが、収束に時間がかかります — ステップサイズ $\mu$ を大きくすれば収束は速くなりますが、定常状態でのミスアジャストメントが増大するというトレードオフがあります
ウィーナフィルタの誤差性能面（MMSE曲面）
ウィーナフィルタの最適性をより深く理解するために、コスト関数 $J(\bm{h})$ の形状を可視化してみましょう。フィルタ長 $M = 2$ の場合は、2次元のコスト関数面（誤差性能面）を3Dプロットで確認できます。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# --- 信号生成 ---
np.random.seed(42)
N = 2048
fs = 256
t = np.arange(N) / fs
d = np.sin(2 * np.pi * 5 * t)
noise_power = 0.5
v = np.sqrt(noise_power) * np.random.randn(N)
x = d + v
# --- M=2 のウィーナフィルタ ---
M = 2
R_xx = np.zeros((M, M))
for i in range(M):
for j in range(M):
shifts = np.arange(max(i, j), N)
R_xx[i, j] = np.mean(x[shifts - i] * x[shifts - j])
r_dx = np.zeros(M)
for i in range(M):
shifts = np.arange(i, N)
r_dx[i] = np.mean(d[shifts] * x[shifts - i])
h_opt = np.linalg.solve(R_xx, r_dx)
J_min = np.mean(d ** 2) - r_dx @ h_opt
# --- コスト関数の計算 ---
h0_range = np.linspace(h_opt[0] - 1.0, h_opt[0] + 1.0, 100)
h1_range = np.linspace(h_opt[1] - 1.0, h_opt[1] + 1.0, 100)
H0, H1 = np.meshgrid(h0_range, h1_range)
J_surface = np.zeros_like(H0)
Rdd_0 = np.mean(d ** 2)
for i in range(H0.shape[0]):
for j in range(H0.shape[1]):
h_vec = np.array([H0[i, j], H1[i, j]])
J_surface[i, j] = Rdd_0 - 2 * r_dx @ h_vec + h_vec @ R_xx @ h_vec
# --- 3Dプロット ---
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(H0, H1, J_surface, cmap='viridis', alpha=0.7, edgecolor='none')
ax.scatter([h_opt[0]], [h_opt[1]], [J_min], color='red', s=100, zorder=5,
label=f'Optimal: h=({h_opt[0]:.3f}, {h_opt[1]:.3f})')
ax.set_xlabel('$h(0)$')
ax.set_ylabel('$h(1)$')
ax.set_zlabel('$J(h)$')
ax.set_title('Error Performance Surface (M=2)')
ax.legend()
plt.tight_layout()
plt.savefig('error_performance_surface.png', dpi=150, bbox_inches='tight')
plt.show()
print(f'Optimal filter coefficients: h = {h_opt}')
print(f'Minimum MSE: J_min = {J_min:.6f}')
このグラフからは、コスト関数 $J(\bm{h})$ が「お椀型」の二次曲面（楕円放物面）であることが明瞭にわかります。赤い点で示された最適解はお椀の底に位置しており、唯一のグローバル最小値 であることが確認できます。
二次曲面であるということは、局所最小値が存在せず、どのような初期値から出発しても勾配降下法（あるいはLMS）で必ずこの最小値に到達できることを意味しています。これはウィーナフィルタの問題がもつ「凸最適化」という構造的な良さであり、LMSアルゴリズムが安定して収束する理論的な保証の源泉です。
MMSE推定の統計的解釈
ここまでの導出を統計的な観点から再整理しておきましょう。この整理は、ウィーナフィルタの理論をカルマンフィルタやベイズ推定に発展させていく際の橋渡しになります。
線形MMSE推定としての位置づけ
ウィーナフィルタは、所望信号 $d(n)$ の 線形最小平均二乗推定量（Linear MMSE estimator, LMMSE） です。「線形」というのは、推定値 $\hat{d}(n)$ が観測信号 $x(n)$ の線形結合として表されるという意味です。
一般に、MMSE推定量（線形制約なし）は条件付き期待値 $E[d(n) | \bm{x}(n)]$ であり、これは非線形な関数になりうります。しかし、信号と雑音が同時正規分布に従う場合には、条件付き期待値は線形関数になるため、線形MMSE推定量は真のMMSE推定量と一致します。
つまり、ガウス信号とガウス雑音の場合、ウィーナフィルタは「全ての推定器（非線形を含む）のなかで最適」です。非ガウスの場合には、非線形推定器がウィーナフィルタを上回る可能性がありますが、ウィーナフィルタは依然として「全ての線形推定器のなかで最適」です。
直交性原理の幾何学的意味
直交性原理 $E[e(n) \, x(n-m)] = 0$ の幾何学的意味をもう一度振り返ります。確率変数の集合を内積 $\langle a, b \rangle = E[ab]$ をもつヒルベルト空間と見なすと、$\hat{d}(n)$ は観測信号 $\{x(n), x(n-1), \ldots\}$ が張る部分空間への $d(n)$ の直交射影です。
射影定理によれば、直交射影は「部分空間内の元のうち、$d(n)$ に最も近いもの」です。そして、射影の残差（推定誤差 $e(n)$）は部分空間に直交します。これが直交性原理の幾何学的実体です。
この見方は、最小二乗法、カルマンフィルタ、ベイズ線形回帰など、MMSE推定に関わるあらゆる分野で共通する数学的構造です。ウィーナフィルタの理論を深く理解しておくことで、これらの発展的なトピックへの接続がスムーズになります。
SN比とフィルタ性能の関係
ウィーナフィルタの性能は、入力のSN比に本質的に依存します。加法性雑音モデルにおけるウィーナフィルタの周波数応答は
$$ H_{\mathrm{opt}}(f) = \frac{S_{dd}(f)}{S_{dd}(f) + S_{vv}(f)} = \frac{1}{1 + \frac{1}{\mathrm{SNR}(f)}} $$
と書けます。ここで $\mathrm{SNR}(f) = S_{dd}(f) / S_{vv}(f)$ は周波数 $f$ におけるSN比です。
この式から、いくつかの重要な性質が読み取れます。
- $\mathrm{SNR}(f) \to \infty$ のとき $H_{\mathrm{opt}}(f) \to 1$: 完全な通過（雑音がないに等しい帯域）
- $\mathrm{SNR}(f) \to 0$ のとき $H_{\mathrm{opt}}(f) \to 0$: 完全な遮断（信号がないに等しい帯域）
- $\mathrm{SNR}(f) = 1$ のとき $H_{\mathrm{opt}}(f) = 0.5$: 半分の重みで通過
ウィーナフィルタは、各周波数でのSN比を「確信度」として利用し、確信度の高い帯域のみを残すという、非常に合理的な戦略を自動的に実行しているのです。
ウィーナフィルタの限界と拡張
ウィーナフィルタの限界
ウィーナフィルタは美しい理論ですが、実用上いくつかの制約を意識しておく必要があります。
1. 定常性の仮定: ウィーナフィルタは信号と雑音が広義定常であることを前提としています。音声のようにスペクトル特性が時間的に変化する信号では、この仮定が成り立ちません。この場合、短時間窓を使って局所的にパワースペクトルを推定し、窓ごとにウィーナフィルタを適用する 時変ウィーナフィルタ が使われます。
2. 統計量の推定: 真のパワースペクトル $S_{dd}(f)$ と $S_{vv}(f)$ は通常未知であり、データから推定する必要があります。推定精度がフィルタ性能に直結するため、パワースペクトルの推定手法（ピリオドグラム、ウェルチ法、自己回帰モデルなど）の選択が重要です。
3. 線形性の制約: ウィーナフィルタは線形フィルタであるため、信号と雑音の間に非線形な関係がある場合には最適ではありません。例えば、インパルス雑音（突発的な大きなスパイクノイズ）の除去には、メディアンフィルタのような非線形手法がより適しています。
カルマンフィルタへの発展
ウィーナフィルタの限界の多くを克服したのが カルマンフィルタ です。カルマンフィルタは、信号を状態空間モデル（状態方程式 + 観測方程式）で表現し、時刻ごとに最適推定を逐次更新します。
ウィーナフィルタが「周波数領域で全時刻のデータを一括処理」するのに対し、カルマンフィルタは「時間領域で1ステップずつ逐次処理」するという違いがあります。両者は定常・ガウスの条件下では等価な結果を与えますが、カルマンフィルタは非定常な系にも適用可能であり、計算効率でも優れています。
ウィーナフィルタの理論を理解しておくことは、カルマンフィルタの「なぜ最適なのか」を深く理解するための不可欠な基盤となります。
まとめ
本記事では、ウィーナフィルタの理論を基礎から丁寧に導出しました。
- MMSE基準 は、推定誤差の二乗平均を最小化するという、直感的で数学的にも扱いやすい最適性の基準です
- 直交性原理 から ウィーナー・ホップ方程式 $R_{dx}(m) = (h * R_{xx})(m)$ が導かれ、フーリエ変換によって周波数領域の閉じた解 $H_{\mathrm{opt}}(f) = S_{dx}(f) / S_{xx}(f)$ が得られます
- 加法性雑音モデルでは $H_{\mathrm{opt}}(f) = S_{dd}(f) / (S_{dd}(f) + S_{vv}(f))$ となり、周波数ごとのSN比に基づく最適な重みづけ として解釈できます
- 因果的ウィーナフィルタ はスペクトル分解によって導かれ、リアルタイム処理の制約下での最適解を与えます
- FIRウィーナフィルタ は正規方程式 $\bm{R}_{xx} \bm{h}_{\mathrm{opt}} = \bm{r}_{dx}$ を解くことで有限長の最適フィルタ係数が得られます
- LMSアルゴリズム はウィーナフィルタの最適解に適応的に収束するアルゴリズムであり、統計量が未知の環境で威力を発揮します
ウィーナフィルタは信号処理における最適推定の出発点であり、ここで得た知見はカルマンフィルタ、整合フィルタ、チャネル等化器など、多くの発展的トピックに直結します。
次のステップとして、以下の記事も参考にしてください。