適応フィルタは、入力信号の統計的性質が時間とともに変化する環境や、フィルタの最適係数が事前に分からない状況において、フィルタ係数を自動的に更新する仕組みです。固定係数のFIRフィルタやIIRフィルタとは異なり、適応フィルタは動作中に逐次的に学習し、最適な係数に収束します。
適応フィルタは、ノイズキャンセリング（能動騒音制御）、エコーキャンセリング（通信）、チャネル等化（ディジタル通信）、システム同定（制御工学）など、非常に幅広い応用を持つ重要な信号処理技術です。
本記事の内容
- 適応フィルタの必要性（時変環境・未知システム）
- ウィーナーフィルタ（最適解）の導出
- MSEの2次形式と誤差曲面
- 最急降下法によるウィーナー解への収束
- LMS（Least Mean Squares）アルゴリズムの導出
- 収束条件 $0 < \mu < 2/\lambda_{\max}$ の導出
- ミスアジャストメントの定義
- NLMS（正規化LMS）アルゴリズム
- Pythonによるノイズキャンセリング・システム同定・エコーキャンセリングの実装
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
適応フィルタの必要性
従来のフィルタ設計（窓関数法、バタワース設計など）では、信号とノイズの周波数特性が既知であることが前提です。しかし現実の多くの場面では、
- ノイズの統計的性質が時間とともに変化する（非定常環境）
- システムの伝達関数が未知（未知システムの同定）
- 信号とノイズの分離が周波数領域だけでは困難（同一帯域の干渉）
といった問題が生じます。適応フィルタは、誤差信号を使ってフィルタ係数をリアルタイムに更新することで、これらの問題に対処します。
適応フィルタの基本構成
適応フィルタの基本構成は以下の通りです。
- 入力信号: $\bm{x}[n] = [x[n], x[n-1], \ldots, x[n-L+1]]^T$（$L$ タップの入力ベクトル）
- フィルタ係数: $\bm{w}[n] = [w_0[n], w_1[n], \ldots, w_{L-1}[n]]^T$
- フィルタ出力: $y[n] = \bm{w}^T[n] \bm{x}[n]$
- 所望信号: $d[n]$
- 誤差信号: $e[n] = d[n] – y[n] = d[n] – \bm{w}^T[n] \bm{x}[n]$
適応アルゴリズムの目的は、誤差信号 $e[n]$ を何らかの意味で最小化するように $\bm{w}[n]$ を更新することです。
ウィーナーフィルタ（最適解）
コスト関数の定義
最も基本的なコスト関数は 平均二乗誤差（MSE: Mean Squared Error）です。
$$ \begin{equation} J(\bm{w}) = E[|e[n]|^2] = E[|d[n] – \bm{w}^T\bm{x}[n]|^2] \end{equation} $$
ここで $E[\cdot]$ は期待値演算子です。
MSEの展開
$J(\bm{w})$ を展開します。
$$ \begin{align} J(\bm{w}) &= E[(d[n] – \bm{w}^T\bm{x}[n])^2] \\ &= E[d^2[n] – 2d[n]\bm{w}^T\bm{x}[n] + \bm{w}^T\bm{x}[n]\bm{x}^T[n]\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} \end{align} $$
ここで以下の量を定義します。
- 入力の自己相関行列: $\bm{R} = E[\bm{x}[n]\bm{x}^T[n]]$（$L \times L$ の正定値対称行列）
- 相互相関ベクトル: $\bm{p} = E[d[n]\bm{x}[n]]$（$L \times 1$ のベクトル）
- 所望信号の電力: $\sigma_d^2 = E[d^2[n]]$
これらを使うと、MSEは次のように書けます。
$$ \begin{equation} J(\bm{w}) = \sigma_d^2 – 2\bm{w}^T\bm{p} + \bm{w}^T\bm{R}\bm{w} \end{equation} $$
これは $\bm{w}$ に関する 2次形式（下に凸の放物面）です。$\bm{R}$ が正定値であるため、$J(\bm{w})$ には唯一の最小値が存在します。
ウィーナー-ホップ方程式の導出
$J(\bm{w})$ を最小にする $\bm{w}$ を求めるため、勾配をゼロとおきます。
$$ \nabla_{\bm{w}} J(\bm{w}) = -2\bm{p} + 2\bm{R}\bm{w} = \bm{0} $$
この勾配の計算を詳しく見ましょう。$J(\bm{w}) = \sigma_d^2 – 2\bm{w}^T\bm{p} + \bm{w}^T\bm{R}\bm{w}$ の各項について、
$$ \begin{align} \frac{\partial}{\partial \bm{w}}(\sigma_d^2) &= \bm{0} \\ \frac{\partial}{\partial \bm{w}}(-2\bm{w}^T\bm{p}) &= -2\bm{p} \\ \frac{\partial}{\partial \bm{w}}(\bm{w}^T\bm{R}\bm{w}) &= 2\bm{R}\bm{w} \quad (\because \bm{R} = \bm{R}^T) \end{align} $$
$\nabla_{\bm{w}} J = \bm{0}$ より、
$$ \bm{R}\bm{w} = \bm{p} $$
これが ウィーナー-ホップ方程式 です。解は
$$ \begin{equation} \bm{w}_{\text{opt}} = \bm{R}^{-1}\bm{p} \end{equation} $$
となります。これがウィーナーフィルタの最適係数であり、MSEを最小にする解です。
最小MSE
最適係数 $\bm{w}_{\text{opt}}$ を $J(\bm{w})$ に代入すると、最小MSEが得られます。
$$ \begin{align} J_{\min} &= \sigma_d^2 – 2\bm{w}_{\text{opt}}^T\bm{p} + \bm{w}_{\text{opt}}^T\bm{R}\bm{w}_{\text{opt}} \\ &= \sigma_d^2 – 2\bm{p}^T\bm{R}^{-1}\bm{p} + \bm{p}^T\bm{R}^{-1}\bm{R}\bm{R}^{-1}\bm{p} \\ &= \sigma_d^2 – \bm{p}^T\bm{R}^{-1}\bm{p} \end{align} $$
$$ \begin{equation} J_{\min} = \sigma_d^2 – \bm{p}^T\bm{w}_{\text{opt}} \end{equation} $$
誤差曲面と固有値分解
誤差曲面の構造
MSE $J(\bm{w})$ を $\bm{w}_{\text{opt}}$ 周りで展開してみましょう。$\bm{v} = \bm{w} – \bm{w}_{\text{opt}}$ とおくと、
$$ \begin{align} J(\bm{w}) &= J(\bm{w}_{\text{opt}} + \bm{v}) \\ &= \sigma_d^2 – 2(\bm{w}_{\text{opt}} + \bm{v})^T\bm{p} + (\bm{w}_{\text{opt}} + \bm{v})^T\bm{R}(\bm{w}_{\text{opt}} + \bm{v}) \\ &= J_{\min} + \bm{v}^T\bm{R}\bm{v} \end{align} $$
最後の等号は、$\bm{R}\bm{w}_{\text{opt}} = \bm{p}$ を使って展開・整理すると得られます。
$\bm{R}$ を固有値分解すると $\bm{R} = \bm{Q}\bm{\Lambda}\bm{Q}^T$（$\bm{\Lambda} = \text{diag}(\lambda_0, \lambda_1, \ldots, \lambda_{L-1})$）です。$\bm{v}’ = \bm{Q}^T\bm{v}$ と座標変換すると、
$$ J = J_{\min} + \sum_{i=0}^{L-1} \lambda_i v_i’^2 $$
これは各固有値方向に独立な放物線の和です。固有値 $\lambda_i$ が大きい方向は曲率が急で、小さい方向は緩やかです。
固有値の広がり（eigenvalue spread）$\chi = \lambda_{\max}/\lambda_{\min}$ が大きいほど、誤差曲面は細長い楕円形になり、収束が遅くなります。
最急降下法
アルゴリズム
ウィーナーフィルタの最適解 $\bm{w}_{\text{opt}} = \bm{R}^{-1}\bm{p}$ を直接計算するには $\bm{R}$ と $\bm{p}$ の事前知識が必要です。適応フィルタでは、代わりに反復法を使います。
最急降下法（Steepest Descent）は、勾配の逆方向にフィルタ係数を更新する方法です。
$$ \begin{equation} \bm{w}[n+1] = \bm{w}[n] – \frac{\mu}{2}\nabla_{\bm{w}} J(\bm{w}[n]) \end{equation} $$
ここで $\mu$ はステップサイズ（学習率）です。係数 $1/2$ は後の式を簡潔にするための慣習です。
勾配 $\nabla_{\bm{w}} J = -2\bm{p} + 2\bm{R}\bm{w}$ を代入すると、
$$ \begin{equation} \bm{w}[n+1] = \bm{w}[n] + \mu(\bm{p} – \bm{R}\bm{w}[n]) \end{equation} $$
収束条件の導出
収束を解析するため、$\bm{v}[n] = \bm{w}[n] – \bm{w}_{\text{opt}}$ を使います。
$$ \begin{align} \bm{v}[n+1] &= \bm{w}[n+1] – \bm{w}_{\text{opt}} \\ &= \bm{w}[n] + \mu(\bm{p} – \bm{R}\bm{w}[n]) – \bm{w}_{\text{opt}} \\ &= \bm{v}[n] + \mu(\bm{p} – \bm{R}(\bm{v}[n] + \bm{w}_{\text{opt}})) \\ &= \bm{v}[n] + \mu\bm{p} – \mu\bm{R}\bm{v}[n] – \mu\bm{R}\bm{w}_{\text{opt}} \\ &= \bm{v}[n] – \mu\bm{R}\bm{v}[n] \quad (\because \bm{R}\bm{w}_{\text{opt}} = \bm{p}) \\ &= (\bm{I} – \mu\bm{R})\bm{v}[n] \end{align} $$
したがって、
$$ \bm{v}[n] = (\bm{I} – \mu\bm{R})^n \bm{v}[0] $$
$\bm{v}[n] \to \bm{0}$（$n \to \infty$）となるためには、$(\bm{I} – \mu\bm{R})$ のすべての固有値の絶対値が1未満でなければなりません。$\bm{R}$ の固有値を $\lambda_i$ とすると、$(\bm{I} – \mu\bm{R})$ の固有値は $1 – \mu\lambda_i$ です。
$$ |1 – \mu\lambda_i| < 1 \quad \Leftrightarrow \quad 0 < \mu\lambda_i < 2 $$
すべての $i$ についてこれが成り立つ条件は、最大固有値 $\lambda_{\max}$ に対して
$$ \begin{equation} 0 < \mu < \frac{2}{\lambda_{\max}} \end{equation} $$
これが最急降下法の 収束条件 です。
LMSアルゴリズムの導出
瞬時勾配推定
最急降下法の更新則 $\bm{w}[n+1] = \bm{w}[n] + \mu(\bm{p} – \bm{R}\bm{w}[n])$ を使うには、$\bm{R}$ と $\bm{p}$ の知識が必要です。しかし実際にはこれらは未知です。
Widrow と Hoff（1960年）は、期待値を瞬時値で近似するという大胆なアイデアを提案しました。
$$ \begin{align} \bm{R} = E[\bm{x}[n]\bm{x}^T[n]] &\approx \bm{x}[n]\bm{x}^T[n] \\ \bm{p} = E[d[n]\bm{x}[n]] &\approx d[n]\bm{x}[n] \end{align} $$
この近似を最急降下法に代入すると、
$$ \begin{align} \bm{w}[n+1] &= \bm{w}[n] + \mu(d[n]\bm{x}[n] – \bm{x}[n]\bm{x}^T[n]\bm{w}[n]) \\ &= \bm{w}[n] + \mu\bm{x}[n](d[n] – \bm{x}^T[n]\bm{w}[n]) \\ &= \bm{w}[n] + \mu\bm{x}[n]e[n] \end{align} $$
$$ \begin{equation} \bm{w}[n+1] = \bm{w}[n] + \mu e[n] \bm{x}[n] \end{equation} $$
これが LMS（Least Mean Squares）アルゴリズム です。
LMSアルゴリズムの手順
- 初期化: $\bm{w}[0] = \bm{0}$
- 各時刻 $n$ で: – フィルタ出力: $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]$
計算量は1反復あたり $O(L)$ の乗算と加算であり、非常に軽量です。
LMSの収束条件
LMSの収束条件は、最急降下法と同様に
$$ 0 < \mu < \frac{2}{\lambda_{\max}} $$
ですが、$\lambda_{\max}$ が未知の場合は、$\text{tr}(\bm{R}) = \sum_{i=0}^{L-1}\lambda_i \ge \lambda_{\max}$ を使って
$$ 0 < \mu < \frac{2}{\text{tr}(\bm{R})} = \frac{2}{L \cdot \sigma_x^2} $$
（$\sigma_x^2 = E[x^2[n]]$ は入力信号の電力）として安全な上限を設定できます。
実用的には、$\mu$ は上限値の 1/10 程度に設定されることが多いです。
ミスアジャストメント
LMSは瞬時勾配を使うため、勾配推定にノイズが含まれます。そのため、収束後も $\bm{w}[n]$ は $\bm{w}_{\text{opt}}$ の周りで揺動し、MSEは $J_{\min}$ より大きな値に落ち着きます。
この定常状態での余剰MSEを excess MSE と呼び、
$$ J_{\text{excess}} = J_{\infty} – J_{\min} $$
で定義されます。ミスアジャストメント $\mathcal{M}$ は、excess MSE と $J_{\min}$ の比です。
$$ \begin{equation} \mathcal{M} = \frac{J_{\text{excess}}}{J_{\min}} \approx \mu \cdot \text{tr}(\bm{R}) = \mu L \sigma_x^2 \end{equation} $$
この近似は $\mu$ が十分小さい場合に成り立ちます。ミスアジャストメントは $\mu$ に比例するので、$\mu$ を小さくすれば定常MSEは改善しますが、収束速度が遅くなるトレードオフがあります。
NLMS（正規化LMS）アルゴリズム
LMSの問題点は、ステップサイズ $\mu$ が固定であるため、入力信号の電力が変動すると収束特性が変わることです。NLMS（Normalized LMS）は入力信号の電力で正規化することでこの問題を解決します。
$$ \begin{equation} \bm{w}[n+1] = \bm{w}[n] + \frac{\tilde{\mu}}{\|\bm{x}[n]\|^2 + \delta} e[n] \bm{x}[n] \end{equation} $$
ここで $\tilde{\mu}$（$0 < \tilde{\mu} < 2$）は正規化ステップサイズ、$\delta$ はゼロ除算防止の小さな正の定数です。
NLMSの導出は、各時刻で以下の制約付き最適化問題を解くことで得られます。
$$ \min_{\bm{w}[n+1]} \|\bm{w}[n+1] – \bm{w}[n]\|^2 \quad \text{s.t.} \quad d[n] = \bm{w}^T[n+1]\bm{x}[n] $$
「前の係数からの変化を最小にしつつ、現在の誤差をゼロにする」という意味です。ラグランジュの未定乗数法で解くと、
$$ \bm{w}[n+1] = \bm{w}[n] + \frac{e[n]}{\|\bm{x}[n]\|^2}\bm{x}[n] $$
これに $\tilde{\mu}$ を掛けて正規化ステップサイズとしたものがNLMSです。
Pythonでの実装
ノイズキャンセリング
ノイズキャンセリングは適応フィルタの最も代表的な応用です。所望信号にノイズが重畳した観測信号から、ノイズの参照信号を使ってノイズ成分を推定・除去します。
import numpy as np
import matplotlib.pyplot as plt
def lms_filter(x, d, L, mu):
"""LMSアルゴリズム
Parameters
----------
x : ndarray
入力信号（参照ノイズ）
d : ndarray
所望信号（ノイズ混入信号）
L : int
フィルタ長
mu : float
ステップサイズ
Returns
-------
y : ndarray
フィルタ出力（推定ノイズ）
e : ndarray
誤差信号（ノイズ除去後の信号）
w_history : ndarray
フィルタ係数の履歴
"""
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 L > 1 else np.array([x[n]])
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
# --- ノイズキャンセリングのシミュレーション ---
np.random.seed(42)
N = 2000 # サンプル数
# 所望信号（正弦波）
n = np.arange(N)
s = np.sin(2 * np.pi * 0.02 * n) + 0.5 * np.sin(2 * np.pi * 0.05 * n)
# ノイズ源（ランダムノイズ）
noise_ref = np.random.randn(N)
# ノイズがパスを通って信号に重畳（未知パスを模擬）
# 未知パス: 3タップFIRフィルタ
h_path = np.array([1.0, -0.5, 0.3])
noise_filtered = np.convolve(noise_ref, h_path, mode='full')[:N]
# 観測信号 = 所望信号 + フィルタ済みノイズ
d_obs = s + noise_filtered
# LMSによるノイズキャンセリング
L = 8 # フィルタ長
mu = 0.01 # ステップサイズ
y_est, e_clean, w_hist = lms_filter(noise_ref, d_obs, L, mu)
# 可視化
fig, axes = plt.subplots(4, 1, figsize=(12, 14))
axes[0].plot(n, s, 'g', linewidth=1.5, label='Original signal s[n]')
axes[0].set_xlabel('Sample')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Original Signal')
axes[0].legend()
axes[0].grid(True)
axes[1].plot(n, d_obs, 'b', alpha=0.7, label='Noisy observation d[n]')
axes[1].set_xlabel('Sample')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Noisy Observation (Signal + Noise)')
axes[1].legend()
axes[1].grid(True)
axes[2].plot(n, e_clean, 'r', linewidth=1.5, label='LMS output e[n]')
axes[2].plot(n, s, 'g--', alpha=0.5, label='Original signal')
axes[2].set_xlabel('Sample')
axes[2].set_ylabel('Amplitude')
axes[2].set_title('Noise-Cancelled Signal (LMS Error Output)')
axes[2].legend()
axes[2].grid(True)
# 学習曲線
mse = e_clean**2
# 移動平均で平滑化
window_size = 50
mse_smooth = np.convolve(mse, np.ones(window_size)/window_size, mode='valid')
axes[3].plot(mse_smooth, 'b')
axes[3].set_xlabel('Sample')
axes[3].set_ylabel('MSE')
axes[3].set_title('Learning Curve (Smoothed MSE)')
axes[3].set_yscale('log')
axes[3].grid(True)
plt.tight_layout()
plt.show()
ノイズキャンセリングでは、適応フィルタが参照ノイズから観測信号中のノイズ成分を推定します。誤差信号 $e[n]$ がノイズ除去後の信号になるという点が直感的でないかもしれませんが、$e[n] = d[n] – y[n] = (s[n] + \text{noise}) – \hat{\text{noise}} \approx s[n]$ と理解できます。
システム同定
未知のシステム（FIRフィルタ）を適応フィルタで同定する実験です。
import numpy as np
import matplotlib.pyplot as plt
def nlms_filter(x, d, L, mu_tilde, delta=1e-6):
"""NLMSアルゴリズム
Parameters
----------
x : ndarray
入力信号
d : ndarray
所望信号
L : int
フィルタ長
mu_tilde : float
正規化ステップサイズ (0 < mu_tilde < 2)
delta : float
正則化パラメータ
Returns
-------
y, e, w_history
"""
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]
# NLMS更新（入力電力で正規化）
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
# --- システム同定 ---
np.random.seed(123)
N = 3000
# 未知システム（同定対象）
h_unknown = np.array([0.5, 1.2, -0.8, 0.3, -0.1, 0.05])
L_true = len(h_unknown)
# 入力信号（白色ガウスノイズ）
x = np.random.randn(N)
# 未知システムの出力
d = np.convolve(x, h_unknown, mode='full')[:N]
# 観測ノイズを追加
d += 0.01 * np.random.randn(N)
# LMSとNLMSの比較
L = 8 # フィルタ長（真の長さより少し大きめ）
# LMS
mu_lms = 0.02
_, e_lms, w_hist_lms = lms_filter(x, d, L, mu_lms)
# NLMS
mu_nlms = 0.5
_, e_nlms, w_hist_nlms = nlms_filter(x, d, L, mu_nlms)
fig, axes = plt.subplots(3, 1, figsize=(12, 12))
# 学習曲線の比較
mse_lms = np.convolve(e_lms**2, np.ones(100)/100, mode='valid')
mse_nlms = np.convolve(e_nlms**2, np.ones(100)/100, mode='valid')
axes[0].plot(mse_lms, 'b', label=f'LMS (μ={mu_lms})')
axes[0].plot(mse_nlms, 'r', label=f'NLMS (μ̃={mu_nlms})')
axes[0].set_xlabel('Sample')
axes[0].set_ylabel('MSE')
axes[0].set_title('Learning Curve: LMS vs NLMS')
axes[0].set_yscale('log')
axes[0].legend()
axes[0].grid(True)
# 収束後のフィルタ係数
h_lms = w_hist_lms[-1, :L_true]
h_nlms = w_hist_nlms[-1, :L_true]
x_pos = np.arange(L_true)
width = 0.25
axes[1].bar(x_pos - width, h_unknown, width, label='True', color='green', alpha=0.8)
axes[1].bar(x_pos, h_lms, width, label='LMS', color='blue', alpha=0.8)
axes[1].bar(x_pos + width, h_nlms, width, label='NLMS', color='red', alpha=0.8)
axes[1].set_xlabel('Tap Index')
axes[1].set_ylabel('Coefficient Value')
axes[1].set_title('Identified Filter Coefficients')
axes[1].legend()
axes[1].grid(True)
# 係数の収束過程
for i in range(L_true):
axes[2].plot(w_hist_nlms[:, i], label=f'w[{i}]')
axes[2].axhline(y=h_unknown[i], color='gray', linestyle='--', alpha=0.3)
axes[2].set_xlabel('Sample')
axes[2].set_ylabel('Coefficient Value')
axes[2].set_title('NLMS: Coefficient Convergence')
axes[2].legend(fontsize=8)
axes[2].grid(True)
plt.tight_layout()
plt.show()
# 同定精度の表示
print('=== System Identification Results ===')
print(f'True coefficients: {h_unknown}')
print(f'LMS identified: {np.round(h_lms, 4)}')
print(f'NLMS identified: {np.round(h_nlms, 4)}')
print(f'LMS error norm: {np.linalg.norm(h_lms - h_unknown):.6f}')
print(f'NLMS error norm: {np.linalg.norm(h_nlms - h_unknown):.6f}')
エコーキャンセリング
通信システムにおけるエコーキャンセリングのシミュレーションです。スピーカーから出た音がマイクに回り込む「音響エコー」を適応フィルタで除去します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as sig
def lms_echo_cancel(far_end, mic, L, mu):
"""エコーキャンセリング用LMS
Parameters
----------
far_end : ndarray
遠端信号（スピーカーに送る信号）
mic : ndarray
マイク信号（エコー + 近端音声）
L : int
適応フィルタ長
mu : float
ステップサイズ
Returns
-------
echo_est : ndarray
推定エコー
error : ndarray
エコー除去後の信号
erle : ndarray
ERLE (Echo Return Loss Enhancement)
"""
N = len(far_end)
w = np.zeros(L)
echo_est = np.zeros(N)
error = np.zeros(N)
for n in range(L, N):
x_vec = far_end[n-L+1:n+1][::-1]
echo_est[n] = np.dot(w, x_vec)
error[n] = mic[n] - echo_est[n]
# NLMS更新
norm_sq = np.dot(x_vec, x_vec) + 1e-6
w = w + (mu / norm_sq) * error[n] * x_vec
# ERLE計算（移動平均）
win = 256
mic_power = np.convolve(mic**2, np.ones(win)/win, mode='same')
err_power = np.convolve(error**2, np.ones(win)/win, mode='same')
erle = 10 * np.log10(mic_power / (err_power + 1e-12) + 1e-12)
return echo_est, error, erle
# --- エコーキャンセリングのシミュレーション ---
np.random.seed(456)
fs = 8000
T = 3.0
N = int(fs * T)
t = np.arange(N) / fs
# 遠端信号（音声を模擬: チャープ信号 + 正弦波）
far_end = np.zeros(N)
# 区間ごとに異なる信号
far_end[:N//3] = sig.chirp(t[:N//3], f0=200, f1=1000, t1=T/3)
far_end[N//3:2*N//3] = np.sin(2*np.pi*300*t[N//3:2*N//3])
far_end[2*N//3:] = sig.chirp(t[2*N//3:], f0=500, f1=200, t1=T/3)
far_end *= 0.8
# エコーパス（室内音響を模擬）
echo_path_len = 128
echo_path = np.zeros(echo_path_len)
# 直接波 + 早期反射 + 後期残響
echo_path[5] = 0.8 # 直接波（遅延5サンプル）
echo_path[15] = -0.4 # 第1反射
echo_path[30] = 0.2 # 第2反射
echo_path[50] = -0.1 # 第3反射
# 指数減衰する残響
np.random.seed(789)
reverb = np.random.randn(echo_path_len) * np.exp(-np.arange(echo_path_len) / 30) * 0.05
echo_path += reverb
# エコー信号
echo = np.convolve(far_end, echo_path, mode='full')[:N]
# 近端音声（存在しない区間と存在する区間を設定）
near_end = np.zeros(N)
# 1.5秒〜2.0秒に近端音声あり
near_start = int(1.5 * fs)
near_end_idx = int(2.0 * fs)
near_end[near_start:near_end_idx] = 0.3 * np.sin(2*np.pi*500*t[near_start:near_end_idx])
# マイク信号 = エコー + 近端音声 + 背景雑音
background_noise = 0.01 * np.random.randn(N)
mic = echo + near_end + background_noise
# エコーキャンセリング
L_filter = 160 # フィルタ長
mu_ec = 0.5 # NLMS ステップサイズ
echo_est, error, erle = lms_echo_cancel(far_end, mic, L_filter, mu_ec)
# 可視化
fig, axes = plt.subplots(5, 1, figsize=(14, 16))
axes[0].plot(t, far_end, 'b')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Far-end Signal (Speaker)')
axes[0].grid(True)
axes[1].plot(t, echo, 'r', alpha=0.7, label='Echo')
axes[1].plot(t, near_end, 'g', alpha=0.7, label='Near-end')
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Echo and Near-end Signal')
axes[1].legend()
axes[1].grid(True)
axes[2].plot(t, mic, 'b', alpha=0.7)
axes[2].set_xlabel('Time (s)')
axes[2].set_ylabel('Amplitude')
axes[2].set_title('Microphone Signal (Echo + Near-end + Noise)')
axes[2].grid(True)
axes[3].plot(t, error, 'r', alpha=0.7, label='After echo cancel')
axes[3].plot(t, near_end, 'g--', alpha=0.7, label='Near-end (reference)')
axes[3].set_xlabel('Time (s)')
axes[3].set_ylabel('Amplitude')
axes[3].set_title('Echo-Cancelled Signal')
axes[3].legend()
axes[3].grid(True)
axes[4].plot(t, erle, 'b')
axes[4].set_xlabel('Time (s)')
axes[4].set_ylabel('ERLE (dB)')
axes[4].set_title('Echo Return Loss Enhancement')
axes[4].set_ylim([0, 40])
axes[4].grid(True)
plt.tight_layout()
plt.show()
エコーキャンセリングでは、適応フィルタが遠端信号から音響エコーパスのインパルス応答を推定し、推定エコーをマイク信号から差し引きます。ERLE（Echo Return Loss Enhancement）はエコー除去の効果を dB で表す指標で、値が大きいほどエコー除去性能が良いことを示します。
まとめ
本記事では、適応フィルタの基礎理論とLMSアルゴリズムについて解説しました。
- ウィーナーフィルタはMSEを最小化する最適解 $\bm{w}_{\text{opt}} = \bm{R}^{-1}\bm{p}$ であり、ウィーナー-ホップ方程式から導かれる
- MSEは $\bm{w}$ の2次形式であり、誤差曲面は楕円形の放物面。固有値の広がりが収束速度に影響する
- LMSアルゴリズムは瞬時勾配推定に基づき、更新則 $\bm{w}[n+1] = \bm{w}[n] + \mu e[n]\bm{x}[n]$ で表される
- 収束条件は $0 < \mu < 2/\lambda_{\max}$ であり、ミスアジャストメントは $\mu$ に比例する
- NLMSは入力電力で正規化することで、入力信号の変動に対するロバスト性を改善する
- ノイズキャンセリング・システム同定・エコーキャンセリングは適応フィルタの代表的な応用
次のステップとして、以下の記事も参考にしてください。