ノイズキャンセリングイヤホンを装着すると、電車の中でも周囲の騒音がすっと消えて音楽が鮮明に聞こえます。この背後では、マイクが拾った外部ノイズをリアルタイムで分析し、それを打ち消す信号を生成する適応フィルタが動作しています。ここで重要なのは、騒音の特性は刻々と変化するという点です。電車の走行音、車内アナウンス、乗客の会話など、ノイズ環境は一定ではありません。あらかじめ設計した固定フィルタでは、こうした時変の環境に対応できないのです。
では、フィルタの係数を環境の変化に追従して自動的に更新する方法はないのでしょうか。この問いに対する答えが、本記事で扱う適応フィルタとその中心的アルゴリズムであるLMS（Least Mean Squares） です。
適応フィルタを理解すると、次のような幅広い応用が見えてきます。
- ノイズキャンセリング: 環境騒音の除去（イヤホン、補聴器、航空機コックピット）
- エコーキャンセリング: 電話回線やWeb会議での音響エコーの除去
- チャネル等化: 無線通信における伝搬路の歪みの補償
- システム同定: 未知のシステムの特性をリアルタイムで推定
- ビームフォーミング: アンテナアレイやマイクアレイにおける方向選択的な信号抽出
本記事の内容
- 適応フィルタの基本構造と動作原理
- ウィーナーフィルタの最適解（ウィーナー・ホップ方程式）
- 最急降下法による係数更新の導出
- LMSアルゴリズム: 確率的勾配降下とその収束条件
- ステップサイズ $\mu$ の設計と収束速度 vs 定常誤差のトレードオフ
- NLMS（正規化LMS）アルゴリズム
- Python実装: ノイズキャンセリング、チャネル等化、収束曲線
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
適応フィルタとは
固定フィルタの限界
ディジタルフィルタの記事で学んだFIRフィルタやIIRフィルタは、設計時にフィルタ係数を決定し、以降は固定して使います。これは、処理対象の信号の統計的性質（周波数スペクトル、ノイズの特性など）が事前にわかっていることが前提です。
しかし、現実の多くの場面では、この前提が成り立ちません。
- 無線通信のチャネル特性は、送信機や受信機の移動、建物の反射などで絶えず変化する
- 音響エコーの経路は、部屋の中で人が動くたびに変わる
- 生体信号のノイズ特性は、被験者の状態や体動によって変動する
こうした未知かつ時変の環境に対応するには、フィルタ係数を観測データに基づいて逐次更新する仕組みが必要です。これが適応フィルタの基本的な発想です。
適応フィルタの基本構造
適応フィルタの構成要素は4つです。
FIRフィルタ: $L$ 個の係数 $\bm{w} = [w_0, w_1, \dots, w_{L-1}]^T$ を持つ横断型（transversal）フィルタ。入力ベクトル $\bm{x}[n] = [x[n], x[n-1], \dots, x[n-L+1]]^T$ に対してフィルタ出力は
$$ y[n] = \bm{w}^T \bm{x}[n] = \sum_{k=0}^{L-1} w_k x[n-k] $$
所望信号（desired signal）: $d[n]$。フィルタ出力が近づくべき目標の信号です。
誤差信号: $e[n] = d[n] – y[n]$。所望信号とフィルタ出力の差です。
適応アルゴリズム: 誤差 $e[n]$ を用いて係数 $\bm{w}$ を更新する規則。誤差を小さくする方向に係数を調整し続けます。
適応フィルタの動作を一言でまとめると、「出力 $y[n]$ が所望信号 $d[n]$ にできるだけ近くなるように、係数 $\bm{w}$ を繰り返し修正する」ということです。この修正の最適な方向を教えてくれるのが、次に説明するウィーナーフィルタの理論です。
ここで重要な疑問が生まれます。「係数をどの値に設定すれば、誤差が最小になるのか？」この最適解がウィーナーフィルタです。
ウィーナーフィルタの最適解
コスト関数の定義
適応フィルタの目標は、誤差信号 $e[n] = d[n] – y[n]$ を何らかの意味で最小化することです。最も自然で数学的に扱いやすいのは、平均二乗誤差（MSE: Mean Squared Error） を最小化することです。
$$ J(\bm{w}) = E\left[|e[n]|^2\right] = E\left[|d[n] – \bm{w}^T\bm{x}[n]|^2\right] $$
ここで $E[\cdot]$ は期待値（統計的平均）を表します。$J(\bm{w})$ はコスト関数（目的関数）と呼ばれ、フィルタ係数 $\bm{w}$ の関数です。
このコスト関数が、なぜ「自然で扱いやすい」のかを直感的に理解しましょう。二乗することで正負の誤差が相殺しないようになり（絶対値と同様の効果）、さらに期待値を取ることで個々のサンプルの偶然の変動に左右されず、統計的に安定した指標になります。加えて、二乗関数は微分可能で凸関数なので、最適化が容易です。
コスト関数の展開
$J(\bm{w})$ を展開してみましょう。
$$ J(\bm{w}) = E\left[(d[n] – \bm{w}^T\bm{x}[n])^2\right] $$
積を展開すると
$$ J(\bm{w}) = E[d^2[n]] – 2\bm{w}^T E[d[n]\bm{x}[n]] + \bm{w}^T E[\bm{x}[n]\bm{x}^T[n]] \bm{w} $$
ここで、以下の量を定義します。
入力の自己相関行列:
$$ \bm{R} = E[\bm{x}[n]\bm{x}^T[n]] $$
$\bm{R}$ は $L \times L$ の正定値対称行列で、入力信号の統計的性質を集約しています。対角要素 $R_{ii} = E[x^2[n-i]]$ は各遅延サンプルのパワー、非対角要素 $R_{ij} = E[x[n-i]x[n-j]]$ は異なる時刻のサンプル間の相関を表します。
入力と所望信号の相互相関ベクトル:
$$ \bm{p} = E[d[n]\bm{x}[n]] $$
$\bm{p}$ は $L$ 次元のベクトルで、所望信号と入力信号の統計的な関連性を表しています。
所望信号のパワー:
$$ \sigma_d^2 = E[d^2[n]] $$
これらを用いると、コスト関数は
$$ J(\bm{w}) = \sigma_d^2 – 2\bm{w}^T\bm{p} + \bm{w}^T\bm{R}\bm{w} $$
と書けます。これは $\bm{w}$ に関する2次関数（放物面）です。$\bm{R}$ が正定値であるため、$J(\bm{w})$ は下に凸で、唯一の最小値を持ちます。この幾何学的な性質は、最適化の上で非常に都合が良いです。
ウィーナー・ホップ方程式の導出
$J(\bm{w})$ を最小にする最適係数 $\bm{w}_{opt}$ は、勾配がゼロになる点で得られます。$\bm{w}$ に関する勾配を計算します。
$$ \nabla_{\bm{w}} J(\bm{w}) = \frac{\partial J}{\partial \bm{w}} = -2\bm{p} + 2\bm{R}\bm{w} $$
ここで、$\sigma_d^2$ は $\bm{w}$ に依存しないので消え、$-2\bm{w}^T\bm{p}$ の勾配は $-2\bm{p}$、$\bm{w}^T\bm{R}\bm{w}$ の勾配は $2\bm{R}\bm{w}$ です。
勾配をゼロに設定すると
$$ -2\bm{p} + 2\bm{R}\bm{w}_{opt} = \bm{0} $$
$$ \boxed{\bm{R}\bm{w}_{opt} = \bm{p}} $$
これがウィーナー・ホップ方程式（Wiener-Hopf equation）です。$\bm{R}$ が正則（逆行列が存在する）なら、最適解は
$$ \bm{w}_{opt} = \bm{R}^{-1}\bm{p} $$
で一意に決まります。このフィルタをウィーナーフィルタと呼びます。
最小平均二乗誤差
最適係数 $\bm{w}_{opt}$ を $J(\bm{w})$ に代入すると、達成可能な最小MSEは
$$ J_{\min} = \sigma_d^2 – \bm{p}^T\bm{R}^{-1}\bm{p} $$
です。$J_{\min} > 0$ であれば、ウィーナーフィルタでも完全には所望信号を再現できないことを意味します。これは入力信号だけでは説明できない所望信号の成分（例えばフィルタ長が不足している、入力とは無関係なノイズがある等）が存在するためです。
ウィーナー解の限界
ウィーナーフィルタは理論的に最適ですが、実用上は2つの大きな問題があります。
問題1: 統計量が未知。$\bm{R}$ と $\bm{p}$ を計算するには、入力信号と所望信号の統計的性質を知る必要があります。しかし実際の環境では、これらは事前にはわかりません。
問題2: 行列の逆演算。$\bm{R}^{-1}$ の計算は $O(L^3)$ の計算量を要し、フィルタ長 $L$ が大きいと計算コストが膨大になります。さらに、環境が時変の場合は統計量が変化するたびに再計算が必要です。
これらの問題を解決するのが、次に説明する反復的な適応アルゴリズムです。行列の逆演算を行わず、データが到着するたびに少しずつ最適解に近づいていくアプローチです。
最急降下法
反復的最適化のアイデア
ウィーナー解を一気に求める代わりに、コスト関数の勾配（最も急な上り坂の方向）の逆方向に少しずつ係数を動かすという戦略を考えます。これが最急降下法（steepest descent method） です。
山の中で霧に囲まれて谷底（最小値）を探す場面を想像してください。周りが見えないので、足元の傾斜だけを頼りに、最も急な下り坂の方向に一歩ずつ進みます。平坦な場所（勾配ゼロ）に到達したら、そこが谷底です。$J(\bm{w})$ が下に凸なので、このような戦略で必ず唯一の最小値に到達できます。
係数更新式の導出
時刻 $n$ での係数ベクトルを $\bm{w}[n]$ とします。最急降下法の更新則は
$$ \bm{w}[n+1] = \bm{w}[n] – \frac{\mu}{2} \nabla_{\bm{w}} J(\bm{w}[n]) $$
ここで $\mu > 0$ はステップサイズ（学習率）パラメータで、各ステップでの移動量を制御します。$1/2$ の係数は勾配式の $2$ と相殺させるための慣例です。
先ほど求めた勾配 $\nabla_{\bm{w}} J = -2\bm{p} + 2\bm{R}\bm{w}$ を代入すると
$$ \bm{w}[n+1] = \bm{w}[n] – \mu(-\bm{p} + \bm{R}\bm{w}[n]) $$
$$ \boxed{\bm{w}[n+1] = \bm{w}[n] + \mu(\bm{p} – \bm{R}\bm{w}[n])} $$
この更新式は、現在の係数 $\bm{w}[n]$ と最適解 $\bm{w}_{opt}$ のずれ（$\bm{R}\bm{w}[n] – \bm{p} = \bm{R}(\bm{w}[n] – \bm{w}_{opt})$）に比例した修正を行っています。
収束条件
最急降下法が安定に収束するためには、ステップサイズ $\mu$ が適切な範囲にある必要があります。収束条件を導出するために、誤差ベクトル $\bm{v}[n] = \bm{w}[n] – \bm{w}_{opt}$ の時間発展を解析します。
更新式から $\bm{w}_{opt}$ を引くと
$$ \bm{v}[n+1] = \bm{v}[n] – \mu\bm{R}\bm{v}[n] = (\bm{I} – \mu\bm{R})\bm{v}[n] $$
$\bm{R}$ の固有値分解 $\bm{R} = \bm{Q}\bm{\Lambda}\bm{Q}^T$（$\bm{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_L)$）を用いて、変換ベクトル $\bm{v}'[n] = \bm{Q}^T\bm{v}[n]$ を定義すると
$$ \bm{v}'[n+1] = (\bm{I} – \mu\bm{\Lambda})\bm{v}'[n] $$
各成分は独立に
$$ v’_k[n+1] = (1 – \mu\lambda_k) v’_k[n] $$
と更新されます。$v’_k[n] = (1 – \mu\lambda_k)^n v’_k[0]$ なので、$n \to \infty$ で $v’_k[n] \to 0$ となるためには
$$ |1 – \mu\lambda_k| < 1 \quad \Longleftrightarrow \quad 0 < \mu < \frac{2}{\lambda_k} $$
これが全ての固有値 $\lambda_k$ で成り立つ必要があるので、最大固有値 $\lambda_{\max}$ に対する条件が最も厳しく
$$ \boxed{0 < \mu < \frac{2}{\lambda_{\max}}} $$
が安定性の必要十分条件です。
収束速度
各成分の収束の時定数は
$$ \tau_k = -\frac{1}{\ln|1 – \mu\lambda_k|} \approx \frac{1}{\mu\lambda_k} \quad (\mu\lambda_k \ll 1 \text{ のとき}) $$
最も遅い収束は最小固有値 $\lambda_{\min}$ に対応する成分で
$$ \tau_{\max} \approx \frac{1}{\mu\lambda_{\min}} $$
です。$\lambda_{\max} / \lambda_{\min}$（固有値スプレッド）が大きいと、一部の成分は速く収束するが別の成分は非常に遅く収束するという不均一な振る舞いが生じます。この固有値スプレッドは入力信号の相関構造に依存し、白色ノイズ入力（$\bm{R} = \sigma^2\bm{I}$）で最も小さく（= 1）、相関の強い入力で大きくなります。
最急降下法は正確な勾配 $\nabla J = -2\bm{p} + 2\bm{R}\bm{w}$ を使いますが、これには $\bm{R}$ と $\bm{p}$ の知識が必要です。しかし先述の通り、実際の環境ではこれらの統計量は未知です。そこで、真の勾配を瞬時値で推定するという大胆な近似が登場します。これがLMSアルゴリズムの核心です。
LMSアルゴリズム
確率的勾配降下の発想
Widrow と Hoff が1960年に提案したLMS（Least Mean Squares）アルゴリズム は、最急降下法の勾配
$$ \nabla J = -2\bm{p} + 2\bm{R}\bm{w} = -2E[d[n]\bm{x}[n]] + 2E[\bm{x}[n]\bm{x}^T[n]]\bm{w} $$
の期待値を、時刻 $n$ の瞬時値で置き換えます。
$$ E[d[n]\bm{x}[n]] \approx d[n]\bm{x}[n], \quad E[\bm{x}[n]\bm{x}^T[n]]\bm{w} \approx \bm{x}[n]\bm{x}^T[n]\bm{w}[n] = \bm{x}[n]y[n] $$
したがって、瞬時勾配の推定値は
$$ \hat{\nabla} J[n] = -2d[n]\bm{x}[n] + 2\bm{x}[n]y[n] = -2(d[n] – y[n])\bm{x}[n] = -2e[n]\bm{x}[n] $$
これを最急降下法の更新式に代入すると
$$ \bm{w}[n+1] = \bm{w}[n] – \frac{\mu}{2}(-2e[n]\bm{x}[n]) $$
$$ \boxed{\bm{w}[n+1] = \bm{w}[n] + \mu \, e[n] \, \bm{x}[n]} $$
これがLMSアルゴリズムの更新式です。
LMSアルゴリズムのステップ
LMSの全体的な手順をまとめます。
初期化: $\bm{w}[0] = \bm{0}$（ゼロベクトルで初期化するのが一般的）
各時刻 $n = 0, 1, 2, \dots$ で:
- フィルタ出力を計算: $y[n] = \bm{w}^T[n] \bm{x}[n]$
- 誤差信号を計算: $e[n] = d[n] – y[n]$
- 係数を更新: $\bm{w}[n+1] = \bm{w}[n] + \mu \, e[n] \, \bm{x}[n]$
わずか3行の式で表される、驚くほどシンプルなアルゴリズムです。各時刻の計算量は $O(L)$（$L$ はフィルタ長）で、ウィーナー解の $O(L^3)$ と比べると桁違いに効率的です。
LMSの直感的理解
LMSの更新式 $\bm{w}[n+1] = \bm{w}[n] + \mu \, e[n] \, \bm{x}[n]$ を直感的に理解しましょう。
誤差 $e[n]$ の役割: $e[n] > 0$（出力が所望信号より小さい）なら、ゲインを上げる方向に修正します。$e[n] < 0$ なら逆です。$e[n] = 0$ なら修正なし。つまり、「間違えた方向と反対に修正する」という自然な学習則です。
入力 $\bm{x}[n]$ の役割: 修正量は入力に比例します。入力が大きいサンプルほど出力に大きな影響を与えるため、そのときの誤差情報をより強く係数に反映させます。逆に、入力がゼロに近いときは係数をほとんど変更しません。
ステップサイズ $\mu$ の役割: 各ステップでの修正量のスケールを決めます。大きすぎると振動して収束しない、小さすぎると収束が遅いというトレードオフがあります。
LMSの収束条件
LMSの収束は確率的であり、最急降下法とは異なる解析が必要です。平均的な意味での収束（$E[\bm{w}[n]] \to \bm{w}_{opt}$）のためには
$$ \boxed{0 < \mu < \frac{2}{\text{tr}(\bm{R})} = \frac{2}{\sum_{k=1}^{L} \lambda_k} = \frac{2}{L \cdot \sigma_x^2}} $$
が必要です。ここで $\text{tr}(\bm{R})$ は $\bm{R}$ のトレース（対角要素の和 = 固有値の和）、$\sigma_x^2 = E[x^2[n]]$ は入力信号のパワーです。最後の等号は、入力が定常かつ各遅延サンプルのパワーが同じ場合に成立します。
この条件は最急降下法の条件 $\mu < 2/\lambda_{\max}$ よりも厳しいです。$\text{tr}(\bm{R}) = \sum \lambda_k \geq \lambda_{\max}$ なので、LMSの安定なステップサイズの範囲は最急降下法よりも狭くなります。瞬時勾配のノイズ（推定誤差）が原因です。
実用上の目安として、安定かつ良好な収束のためには
$$ \mu \approx \frac{1}{10 \cdot L \cdot \sigma_x^2} $$
程度に設定することが多いです。
ステップサイズのトレードオフ
LMSアルゴリズムにおいて、ステップサイズ $\mu$ は最も重要な設計パラメータです。$\mu$ が収束速度と定常誤差の両方に影響を与えるため、用途に応じた適切な値の選択が不可欠です。
大きな $\mu$（上限に近い値）の場合:
- 収束速度が速い: 少ないサンプル数でウィーナー解の近傍に到達
- 定常誤差が大きい: 収束後も係数が最適値の周りで大きく変動（余剰MSE、excess MSEが大きい）
- 時変環境への追従性が高い: 環境変化を素早くキャッチ
小さな $\mu$（上限から遠い値）の場合:
- 収束速度が遅い: ウィーナー解への到達に多くのサンプルが必要
- 定常誤差が小さい: 収束後の係数変動が小さく、最適解に近い性能を達成
- 時変環境への追従性が低い: 環境変化への反応が鈍い
定常状態での余剰MSEは
$$ J_{excess} = J_{ss} – J_{\min} \approx \mu \cdot J_{\min} \cdot L / 2 $$
で近似されます。ここで $J_{ss}$ は定常MSE、$J_{\min}$ はウィーナーフィルタの最小MSEです。この式から、余剰MSEは $\mu$ に比例して増加し、フィルタ長 $L$ にも比例することがわかります。
ミスアジャストメント（misadjustment） $M$ は、余剰MSEの $J_{\min}$ に対する比で定義されます。
$$ M = \frac{J_{excess}}{J_{\min}} \approx \frac{\mu L}{2} \cdot \frac{\text{tr}(\bm{R})}{\text{tr}(\bm{R})} = \mu \cdot \frac{L \sigma_x^2}{2} $$
$M = 0.1$（10%のミスアジャストメント）を目標とすると、$\mu = 0.2 / (L\sigma_x^2)$ と設定します。
このように、$\mu$ の選択は「速さ vs 精度 vs 追従性」の3要素のバランスを決める設計判断です。一般的には、定常環境では小さな $\mu$、時変環境では大きな $\mu$ が適しています。
LMSの大きな利点はその単純さですが、入力信号のパワーが変動すると収束特性が不安定になるという弱点があります。この問題を解決するのがNLMSアルゴリズムです。
NLMS（正規化LMS）アルゴリズム
LMSの問題点
LMSの更新量 $\mu \, e[n] \, \bm{x}[n]$ は入力ベクトルの大きさ $\|\bm{x}[n]\|$ に比例します。入力のパワーが時間的に変動する場合（音声信号のように無音区間と有音区間がある場合など）、修正量が不安定になります。
- 入力パワーが大きいとき: 修正量が過大になり、係数が暴れる
- 入力パワーが小さいとき: 修正量が不足し、適応が遅い
NLMSの導出
この問題を解決するために、修正量を入力ベクトルのノルムで正規化します。
時刻 $n$ での修正を「次の出力誤差をゼロにする」ように最適化する問題を考えます。$\bm{w}[n+1]$ が $\bm{w}[n]$ からなるべく動かないようにしつつ、次の時刻の出力で $d[n]$ を完全に再現する条件を課します。
$$ \min_{\bm{w}[n+1]} \|\bm{w}[n+1] – \bm{w}[n]\|^2 \quad \text{subject to} \quad \bm{w}^T[n+1]\bm{x}[n] = d[n] $$
ラグランジュの未定乗数法で解くと
$$ \bm{w}[n+1] = \bm{w}[n] + \frac{e[n]}{\|\bm{x}[n]\|^2} \bm{x}[n] $$
が得られます。ここに緩和パラメータ $\tilde{\mu}$（$0 < \tilde{\mu} < 2$）を導入して
$$ \boxed{\bm{w}[n+1] = \bm{w}[n] + \frac{\tilde{\mu}}{\|\bm{x}[n]\|^2 + \delta} e[n] \, \bm{x}[n]} $$
これがNLMSアルゴリズムです。$\delta > 0$ は $\|\bm{x}[n]\| = 0$ のときのゼロ除算を防ぐための正則化パラメータ（典型的には $10^{-8}$ 程度の小さな値）です。
NLMSの利点
NLMSにはLMSに対する明確な利点があります。
ステップサイズの設定が容易: LMSでは入力パワー $\sigma_x^2$ に依存して安定なステップサイズの範囲が変わるため、入力特性を知らないと $\mu$ を適切に設定できません。NLMSでは入力パワーで正規化されるため、$\tilde{\mu}$ は入力特性に依存せず、$0 < \tilde{\mu} < 2$ の範囲で選べば安定です。
パワー変動への頑健性: 入力パワーが大きいときは修正量が自動的に抑制され、小さいときは増幅されます。音声信号やレーダーのようにパワー変動の大きい信号に対して、LMSよりも安定した収束が得られます。
LMSとの関係: $\tilde{\mu} / \|\bm{x}[n]\|^2 = \mu$ と対応付けると、NLMSは各時刻でステップサイズを $\mu_{eff}[n] = \tilde{\mu} / \|\bm{x}[n]\|^2$ と適応的に変化させるLMSとみなせます。
ここまでで、適応フィルタの理論（ウィーナー解 → 最急降下法 → LMS → NLMS）を一通り導出しました。次に、Python で実際にLMSとNLMSを実装し、ノイズキャンセリングとチャネル等化という2つの応用例でその動作を確認しましょう。
Python実装: ノイズキャンセリング
シナリオ
ノイズキャンセリングの基本構成を実装します。目的信号 $s[n]$ にノイズ $v[n]$ が加わった観測信号 $d[n] = s[n] + v[n]$ があり、ノイズの参照信号 $x[n]$（ノイズと相関があるが目的信号とは無相関な信号）が別に得られるとします。適応フィルタはこの参照信号からノイズの推定値を生成し、観測信号から差し引いて目的信号を復元します。
import numpy as np
import matplotlib.pyplot as plt
def lms_filter(x, d, L, mu):
"""LMSアルゴリズム
x: 参照入力信号
d: 所望信号（観測信号）
L: フィルタ長
mu: ステップサイズ
"""
N = len(x)
w = np.zeros(L) # フィルタ係数
y = np.zeros(N) # フィルタ出力
e = np.zeros(N) # 誤差信号
w_history = np.zeros((N, L)) # 係数の履歴
for n in range(L, N):
# 入力ベクトル
x_vec = x[n:n-L:-1] if n >= L else np.concatenate([x[n::-1], np.zeros(L-n-1)])
x_vec = x[n-L+1:n+1][::-1]
# フィルタ出力
y[n] = np.dot(w, x_vec)
# 誤差
e[n] = d[n] - y[n]
# 係数更新
w = w + mu * e[n] * x_vec
# 履歴記録
w_history[n] = w.copy()
return y, e, w_history
def nlms_filter(x, d, L, mu_tilde, delta=1e-8):
"""NLMSアルゴリズム
x: 参照入力信号
d: 所望信号
L: フィルタ長
mu_tilde: 正規化ステップサイズ (0 < mu_tilde < 2)
delta: 正則化パラメータ
"""
N = len(x)
w = np.zeros(L)
y = np.zeros(N)
e = np.zeros(N)
w_history = np.zeros((N, L))
for n in range(L, N):
x_vec = x[n-L+1:n+1][::-1]
y[n] = np.dot(w, x_vec)
e[n] = d[n] - y[n]
# 正規化されたステップサイズ
norm_sq = np.dot(x_vec, x_vec) + delta
w = w + (mu_tilde / norm_sq) * e[n] * x_vec
w_history[n] = w.copy()
return y, e, w_history
このコードでは、LMSとNLMSの2つのアルゴリズムを関数として定義しました。どちらも入力のforループで1サンプルずつ処理する逐次処理構造になっており、リアルタイム処理の実際の動作を忠実に再現しています。w_historyに係数の履歴を記録することで、収束過程の可視化が可能です。
ノイズキャンセリングのシミュレーション
import numpy as np
import matplotlib.pyplot as plt
# 信号の生成
np.random.seed(42)
N = 5000 # サンプル数
fs = 8000 # サンプリング周波数
# 目的信号: 正弦波（100 Hz）
t = np.arange(N) / fs
s = np.sin(2 * np.pi * 100 * t)
# ノイズ源: 白色ガウスノイズ
noise_source = np.random.randn(N)
# ノイズの伝搬路（未知のシステム）: FIRフィルタ
h_true = np.array([1.0, -0.5, 0.3, -0.1, 0.05]) # 5タップ
v = np.convolve(noise_source, h_true, mode='full')[:N]
# 観測信号（目的信号 + ノイズ）
d = s + v
# 参照入力: ノイズ源そのもの
x = noise_source
# LMSフィルタの適用
L = 10 # フィルタ長
mu = 0.01 # ステップサイズ
y_lms, e_lms, w_hist_lms = lms_filter(x, d, L, mu)
# NLMSフィルタの適用
mu_tilde = 0.5
y_nlms, e_nlms, w_hist_nlms = nlms_filter(x, d, L, mu_tilde)
fig, axes = plt.subplots(4, 1, figsize=(14, 16))
# 元の信号とノイズ混入信号
ax = axes[0]
ax.plot(t, s, 'b-', linewidth=1, label='Clean signal s[n]', alpha=0.7)
ax.plot(t, d, 'r-', linewidth=0.5, alpha=0.4, label='Noisy observation d[n]')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title('Original Signal and Noisy Observation')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, N/fs)
# LMSノイズキャンセリング結果
ax = axes[1]
ax.plot(t, s, 'b-', linewidth=1, label='Clean signal', alpha=0.5)
ax.plot(t, e_lms, 'g-', linewidth=1, label='LMS output (error = cleaned)')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title(f'LMS Noise Cancelling (L={L}, mu={mu})')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, N/fs)
# NLMSノイズキャンセリング結果
ax = axes[2]
ax.plot(t, s, 'b-', linewidth=1, label='Clean signal', alpha=0.5)
ax.plot(t, e_nlms, 'm-', linewidth=1, label='NLMS output (error = cleaned)')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title(f'NLMS Noise Cancelling (L={L}, mu_tilde={mu_tilde})')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, N/fs)
# 学習曲線（MSE）
ax = axes[3]
# 移動平均でMSEを平滑化
window = 50
mse_lms = np.convolve(e_lms**2, np.ones(window)/window, mode='same')
mse_nlms = np.convolve(e_nlms**2, np.ones(window)/window, mode='same')
mse_noisy = np.convolve((d-s)**2, np.ones(window)/window, mode='same')
ax.semilogy(np.arange(N), mse_noisy, 'r-', linewidth=0.5, alpha=0.3,
label='No filtering')
ax.semilogy(np.arange(N), mse_lms, 'g-', linewidth=1.5, label=f'LMS (mu={mu})')
ax.semilogy(np.arange(N), mse_nlms, 'm-', linewidth=1.5,
label=f'NLMS (mu_tilde={mu_tilde})')
ax.set_xlabel('Sample index n')
ax.set_ylabel('MSE (smoothed)')
ax.set_title('Learning Curve: MSE vs Time')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
ax.set_xlim(0, N)
plt.tight_layout()
plt.savefig('noise_cancelling_demo.png', dpi=150, bbox_inches='tight')
plt.show()
ノイズキャンセリングのシミュレーション結果から、LMS/NLMSの動作が明確に確認できます。
-
最上段（元信号と観測信号）: 100 Hzの正弦波（青）にノイズが加わり、観測信号（赤）からは元の波形がほとんど識別できない状態です。SN比はおよそ 0 dB 程度で、厳しいノイズ環境を模擬しています。
-
2段目（LMS結果）: 初期の数百サンプルではフィルタがまだ収束していないためノイズが残っていますが、約500サンプル（0.06秒）以降は元の正弦波がきれいに復元されています。フィルタが参照ノイズからノイズの伝搬路 $h_{true}$ を学習し、その推定値を観測信号から差し引いています。
-
3段目（NLMS結果）: LMSよりも初期収束が速く、約200サンプル程度で安定した出力が得られています。NLMSの正規化効果により、入力パワーに依存しない安定した適応が実現されています。
-
最下段（学習曲線）: MSEの時間変化を対数スケールで表示しています。赤い線（フィルタなし）のMSEは一定ですが、LMS（緑）とNLMS（マゼンタ）は時間とともにMSEが減少し、定常状態に収束しています。NLMSの方が初期収束が速いことが視覚的にも確認できます。定常状態ではNLMSのMSEがやや大きいですが、これはステップサイズ $\tilde{\mu} = 0.5$ が比較的大きいためです。
ステップサイズの影響
ステップサイズ $\mu$ を変えた場合の収束速度と定常誤差のトレードオフを可視化します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
N = 5000
fs = 8000
t = np.arange(N) / fs
s = np.sin(2 * np.pi * 100 * t)
noise_source = np.random.randn(N)
h_true = np.array([1.0, -0.5, 0.3, -0.1, 0.05])
v = np.convolve(noise_source, h_true, mode='full')[:N]
d = s + v
x = noise_source
L = 10
mus = [0.001, 0.005, 0.01, 0.05, 0.1]
fig, axes = plt.subplots(2, 1, figsize=(14, 10))
# 学習曲線の比較
ax = axes[0]
window = 100
for mu in mus:
_, e, _ = lms_filter(x, d, L, mu)
mse = np.convolve(e**2, np.ones(window)/window, mode='same')
ax.semilogy(np.arange(N), mse, linewidth=1.5, label=f'mu={mu}')
ax.set_xlabel('Sample index n')
ax.set_ylabel('MSE (smoothed)')
ax.set_title('LMS Learning Curve: Effect of Step Size mu')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
ax.set_xlim(0, N)
# 定常MSE vs ステップサイズ
ax = axes[1]
steady_mses = []
for mu in mus:
_, e, _ = lms_filter(x, d, L, mu)
steady_mse = np.mean(e[N//2:]**2) # 後半のMSEを定常値とする
steady_mses.append(steady_mse)
ax.plot(mus, steady_mses, 'bo-', markersize=8, linewidth=2)
ax.set_xlabel('Step size mu')
ax.set_ylabel('Steady-state MSE')
ax.set_title('Steady-state MSE vs Step Size')
ax.set_xscale('log')
ax.set_yscale('log')
ax.grid(True, which='both', alpha=0.3)
# 理論値（近似）: J_ss ≈ J_min * (1 + mu*L*sigma_x^2/2)
sigma_x2 = np.var(x)
J_min = np.mean(s[N//2:]**2) * 0.01 # 近似的な最小MSE
mu_arr = np.logspace(-3, -0.5, 50)
J_theory = J_min * (1 + mu_arr * L * sigma_x2 / 2) + J_min
ax.plot(mu_arr, J_theory, 'r--', linewidth=1.5, alpha=0.7, label='Theoretical trend')
ax.legend()
plt.tight_layout()
plt.savefig('mu_tradeoff.png', dpi=150, bbox_inches='tight')
plt.show()
ステップサイズの比較結果から、$\mu$ の選択が性能に与える影響が定量的にわかります。
-
上図（学習曲線）: $\mu = 0.1$ は最も速く収束していますが、定常状態のMSEが大きく、収束後も大きく変動しています。$\mu = 0.001$ は収束が非常に遅く、5000サンプルでもまだ完全には収束していません。$\mu = 0.01$ は適度な収束速度と低い定常MSEを両立しており、このシナリオでは良いバランスです。
-
下図（定常MSE vs ステップサイズ）: ステップサイズと定常MSEがほぼ比例関係にあることが確認できます。$\mu$ を10倍にすると定常MSEも概ね比例して増加します。赤の破線は理論的なトレンドを示しており、実験結果とよく一致しています。
Python実装: チャネル等化
シナリオ
無線通信では、送信信号がマルチパスチャネル（反射・散乱による複数経路）を通って受信されます。チャネルの伝達特性が未知の場合、受信信号から元の送信信号を復元するために適応等化器が使われます。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)
N = 3000
# 送信信号: BPSK (+-1)
tx_symbols = 2 * (np.random.randint(0, 2, N) - 0.5) # +1 or -1
# マルチパスチャネル（未知）
h_channel = np.array([1.0, 0.5, -0.3, 0.1])
channel_len = len(h_channel)
# チャネル通過後の信号
rx_signal = np.convolve(tx_symbols, h_channel, mode='full')[:N]
# AWGN（加法性白色ガウスノイズ）
SNR_dB = 20
noise_power = np.mean(rx_signal**2) / (10**(SNR_dB/10))
noise = np.sqrt(noise_power) * np.random.randn(N)
rx_noisy = rx_signal + noise
# 適応等化器（LMS）
L_eq = 15 # 等化器のフィルタ長
delay = 7 # 判定遅延（等化器の中心タップ）
mu_eq = 0.01
# トレーニング期間: 既知のシンボルで学習
N_train = 500
w_eq = np.zeros(L_eq)
y_eq = np.zeros(N)
e_eq = np.zeros(N)
decisions = np.zeros(N)
for n in range(L_eq, N):
x_vec = rx_noisy[n-L_eq+1:n+1][::-1]
y_eq[n] = np.dot(w_eq, x_vec)
# トレーニング期間は真のシンボルを所望信号に使用
if n < N_train:
d_n = tx_symbols[n - delay] if n >= delay else 0
else:
# 判定指向モード: 判定結果を所望信号に使用
d_n = np.sign(y_eq[n])
e_eq[n] = d_n - y_eq[n]
w_eq = w_eq + mu_eq * e_eq[n] * x_vec
decisions[n] = np.sign(y_eq[n])
# ビット誤り率の計算
valid_range = slice(N_train, N)
ber_before = np.mean(np.sign(rx_noisy[valid_range]) != tx_symbols[valid_range])
delayed_tx = np.roll(tx_symbols, delay)
ber_after = np.mean(decisions[valid_range] != delayed_tx[valid_range])
fig, axes = plt.subplots(3, 2, figsize=(14, 14))
# チャネルのインパルス応答
ax = axes[0, 0]
ax.stem(np.arange(channel_len), h_channel, linefmt='b-', markerfmt='bo', basefmt='k-')
ax.set_xlabel('Tap index')
ax.set_ylabel('Amplitude')
ax.set_title('Unknown Channel Impulse Response')
ax.grid(True, alpha=0.3)
# 等化器のインパルス応答
ax = axes[0, 1]
ax.stem(np.arange(L_eq), w_eq, linefmt='r-', markerfmt='ro', basefmt='k-')
ax.set_xlabel('Tap index')
ax.set_ylabel('Amplitude')
ax.set_title('Adaptive Equalizer Coefficients (after convergence)')
ax.grid(True, alpha=0.3)
# 送信信号 vs チャネル出力（一部）
ax = axes[1, 0]
show_range = slice(100, 200)
ax.step(np.arange(100), tx_symbols[show_range], 'b-', linewidth=1.5,
label='Transmitted', where='mid')
ax.plot(np.arange(100), rx_noisy[show_range], 'r-', linewidth=0.8,
alpha=0.7, label='Received (channel + noise)')
ax.set_xlabel('Symbol index')
ax.set_ylabel('Amplitude')
ax.set_title('Transmitted vs Received Signal')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 等化後の信号
ax = axes[1, 1]
show_range2 = slice(N_train+100, N_train+200)
ax.step(np.arange(100), delayed_tx[show_range2], 'b-', linewidth=1.5,
label='Transmitted (delayed)', where='mid')
ax.plot(np.arange(100), y_eq[show_range2], 'g-', linewidth=1,
alpha=0.8, label='Equalizer output')
ax.axhline(y=0, color='gray', linestyle='--', alpha=0.3)
ax.set_xlabel('Symbol index')
ax.set_ylabel('Amplitude')
ax.set_title(f'Equalized Signal (BER: {ber_before:.4f} -> {ber_after:.4f})')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# 学習曲線
ax = axes[2, 0]
window = 50
mse_eq = np.convolve(e_eq**2, np.ones(window)/window, mode='same')
ax.semilogy(np.arange(N), mse_eq, 'g-', linewidth=1)
ax.axvline(x=N_train, color='red', linestyle='--', alpha=0.5,
label=f'Training ends (n={N_train})')
ax.set_xlabel('Symbol index n')
ax.set_ylabel('MSE (smoothed)')
ax.set_title('Equalizer Learning Curve')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
# アイパターン（等化前後）
ax = axes[2, 1]
samples_per_symbol = 1 # シンボルレート = サンプリングレート
# 等化前
for i in range(N_train, N_train+200, 2):
if i+2 < N:
ax.plot([0, 1], [rx_noisy[i], rx_noisy[i+1]], 'r-', alpha=0.1, linewidth=0.5)
# 等化後
for i in range(N_train, N_train+200, 2):
if i+2 < N:
ax.plot([2, 3], [y_eq[i], y_eq[i+1]], 'g-', alpha=0.2, linewidth=0.5)
ax.set_xlabel('Symbol period')
ax.set_ylabel('Amplitude')
ax.set_title('Eye Diagram: Before (red, left) vs After (green, right) Equalization')
ax.axhline(y=0, color='gray', linestyle='--', alpha=0.3)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('channel_equalization.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"BER before equalization: {ber_before:.4f}")
print(f"BER after equalization: {ber_after:.4f}")
チャネル等化のシミュレーション結果から、適応フィルタの通信応用における効果が明確にわかります。
-
左上（チャネル応答）: 未知のマルチパスチャネルが4タップのFIRフィルタとしてモデル化されています。メインパス（タップ0）に加えて、遅延した反射波（タップ1〜3）が存在し、符号間干渉（ISI）を引き起こします。
-
右上（等化器の係数）: 収束後の適応等化器の15タップの係数を示しています。この等化器はチャネルの逆特性を近似しており、チャネルと等化器の畳み込みがデルタ関数（単位インパルス）に近づくように学習されています。
-
左中（受信信号）: チャネルを通過した受信信号（赤）は、符号間干渉により送信信号（青の $\pm 1$）から大きく歪んでおり、単純な判定では多くのビット誤りが生じます。
-
右中（等化後）: 等化器の出力（緑）は送信信号（遅延付き）に近い $\pm 1$ の値を示しており、ビット誤り率（BER）が大幅に改善されています。等化前の高いBERから等化後は低いBERへの改善が数値で確認できます。
-
左下（学習曲線）: トレーニング期間（赤の破線の左）でMSEが急速に減少し、トレーニング終了後の判定指向モードでもMSEが安定に維持されていることがわかります。
-
右下（アイパターン）: 等化前（赤）は目が大きく閉じており判定マージンが小さいのに対し、等化後（緑）は目が開いており判定マージンが確保されています。
LMS vs NLMS の収束比較
最後に、入力パワーが変動する環境でLMSとNLMSの収束特性を比較します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(77)
N = 8000
fs = 8000
# パワーが変動する入力信号
t = np.arange(N) / fs
# 前半: 低パワー、後半: 高パワー
amplitude = np.ones(N)
amplitude[:N//2] = 0.1
amplitude[N//2:] = 1.0
# 滑らかに遷移
transition = 200
center = N // 2
for i in range(max(0, center-transition), min(N, center+transition)):
amplitude[i] = 0.1 + 0.9 * (1 + np.tanh((i - center) / (transition/4))) / 2
noise_ref = amplitude * np.random.randn(N)
# 未知システム
h_true = np.array([0.8, -0.4, 0.2, -0.1, 0.05, -0.02])
noise_filtered = np.convolve(noise_ref, h_true, mode='full')[:N]
# 目的信号
signal_clean = 0.5 * np.sin(2 * np.pi * 200 * t)
d = signal_clean + noise_filtered
x = noise_ref
L = 12
# LMS (固定ステップサイズ、低パワー時に合わせて設定)
mu_lms = 0.05
_, e_lms, w_lms = lms_filter(x, d, L, mu_lms)
# NLMS
mu_tilde = 0.5
_, e_nlms, w_nlms = nlms_filter(x, d, L, mu_tilde)
fig, axes = plt.subplots(3, 1, figsize=(14, 12))
# 入力パワーの変動
ax = axes[0]
input_power = np.convolve(x**2, np.ones(100)/100, mode='same')
ax.plot(np.arange(N)/fs, input_power, 'k-', linewidth=1.5)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Input power (smoothed)')
ax.set_title('Input Signal Power Variation')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, N/fs)
# MSE比較
ax = axes[1]
window = 100
mse_lms = np.convolve(e_lms**2, np.ones(window)/window, mode='same')
mse_nlms = np.convolve(e_nlms**2, np.ones(window)/window, mode='same')
ax.semilogy(np.arange(N)/fs, mse_lms, 'g-', linewidth=1.5, label=f'LMS (mu={mu_lms})')
ax.semilogy(np.arange(N)/fs, mse_nlms, 'm-', linewidth=1.5,
label=f'NLMS (mu_tilde={mu_tilde})')
ax.axvline(x=0.5, color='gray', linestyle='--', alpha=0.5, label='Power transition')
ax.set_xlabel('Time [s]')
ax.set_ylabel('MSE (smoothed)')
ax.set_title('LMS vs NLMS: Convergence with Varying Input Power')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
ax.set_xlim(0, N/fs)
# 出力信号の比較
ax = axes[2]
show = slice(N//2-500, N//2+500)
t_show = np.arange(show.start, show.stop) / fs
ax.plot(t_show, signal_clean[show], 'b-', linewidth=1.5, alpha=0.5, label='Clean signal')
ax.plot(t_show, e_lms[show], 'g-', linewidth=1, alpha=0.7, label='LMS output')
ax.plot(t_show, e_nlms[show], 'm-', linewidth=1, alpha=0.7, label='NLMS output')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude')
ax.set_title('Output Signal Around Power Transition')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('lms_vs_nlms_power_variation.png', dpi=150, bbox_inches='tight')
plt.show()
入力パワー変動環境でのLMS vs NLMSの比較から、NLMSの優位性が明確にわかります。
-
上段（入力パワー）: 前半（$t < 0.5$ s）ではパワーが約 0.01、後半（$t > 0.5$ s）ではパワーが約 1.0 と、100倍もの変動があります。これは音声信号の無音/有音や、通信でのフェージングに相当する状況です。
-
中段（MSE比較）: LMS（緑）は前半の低パワー区間での収束が非常に遅く、後半の高パワー区間では逆にステップサイズが実効的に大きくなってMSEが不安定になっています。NLMS（マゼンタ）は入力パワーによる正規化のおかげで、前半も後半も安定した収束を示しています。パワー遷移点（灰色の破線）前後でもNLMSのMSEはスムーズに推移しています。
-
下段（パワー遷移付近の出力）: パワー遷移点の前後500サンプルの出力を比較しています。NLMSは遷移前後で一貫して元の信号に近い出力を維持しているのに対し、LMSは遷移後に一時的に誤差が増大しています。
まとめ
本記事では、適応フィルタの理論的基盤からLMS/NLMSアルゴリズムの導出、実践的な応用までを体系的に解説しました。
- ウィーナーフィルタ: MSEを最小化する最適解 $\bm{w}_{opt} = \bm{R}^{-1}\bm{p}$ は、自己相関行列 $\bm{R}$ と相互相関ベクトル $\bm{p}$ で一意に決まる
- 最急降下法: 真の勾配 $\nabla J = -2\bm{p} + 2\bm{R}\bm{w}$ を用いた反復的最適化。収束条件は $0 < \mu < 2/\lambda_{\max}$
- LMSアルゴリズム: 瞬時勾配で真の勾配を推定する確率的勾配降下法。更新式 $\bm{w}[n+1] = \bm{w}[n] + \mu e[n]\bm{x}[n]$ は $O(L)$ の計算量で逐次処理可能
- ステップサイズ $\mu$: 収束速度と定常誤差のトレードオフを支配する最重要パラメータ。ミスアジャストメント $M \approx \mu L\sigma_x^2 / 2$
- NLMSアルゴリズム: 入力パワーで正規化することで、パワー変動に頑健な適応を実現。$\tilde{\mu} \in (0, 2)$ で安定
- 応用: ノイズキャンセリング（参照ノイズからノイズ経路を学習）、チャネル等化（受信信号からチャネル逆特性を学習）
LMSアルゴリズムは、その単純さと計算効率により、半世紀以上にわたって信号処理の現場で使われ続けています。本記事で扱った基本的なLMS/NLMSは適応フィルタの入口であり、さらに高度なアルゴリズムとしてRLS（再帰的最小二乗法）、アフィン射影法、カルマンフィルタベースの適応フィルタなどが存在します。
次のステップとして、以下の記事も参考にしてください。