適応フィルタ（LMS/RLS）の理論とPython実装

import numpy as np

def lms_filter(x, d, M, mu):
    """LMSアルゴリズムによる適応フィルタ

    Parameters
    ----------
    x : array_like
        入力信号（参照信号）
    d : array_like
        目的信号
    M : int
        フィルタ次数
    mu : float
        ステップサイズ

    Returns
    -------
    y : ndarray
        フィルタ出力
    e : ndarray
        誤差信号
    w_history : ndarray
        フィルタ係数の履歴
    """
    N = len(x)
    w = np.zeros(M)
    y = np.zeros(N)
    e = np.zeros(N)
    w_history = np.zeros((N, M))

    for n in range(M, N):
        x_vec = x[n:n-M:-1] if n >= M else np.pad(x[:n+1][::-1], (0, M-n-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

    return y, e, w_history

def nlms_filter(x, d, M, mu=0.5, eps=1e-8):
    """NLMSアルゴリズムによる適応フィルタ"""
    N = len(x)
    w = np.zeros(M)
    y = np.zeros(N)
    e = np.zeros(N)

    for n in range(M, N):
        x_vec = x[n:n-M:-1] if n >= M else np.pad(x[:n+1][::-1], (0, M-n-1))
        y[n] = np.dot(w, x_vec)
        e[n] = d[n] - y[n]
        norm = np.dot(x_vec, x_vec) + eps
        w = w + (mu / norm) * e[n] * x_vec

    return y, e

def rls_filter(x, d, M, lam=0.99, delta=1.0):
    """RLSアルゴリズムによる適応フィルタ

    Parameters
    ----------
    x : array_like
        入力信号
    d : array_like
        目的信号
    M : int
        フィルタ次数
    lam : float
        忘却係数（0 < lam <= 1）
    delta : float
        逆相関行列の初期スケーリング

    Returns
    -------
    y : ndarray
        フィルタ出力
    e : ndarray
        誤差信号
    """
    N = len(x)
    w = np.zeros(M)
    P = np.eye(M) / delta
    y = np.zeros(N)
    e = np.zeros(N)

    for n in range(M, N):
        x_vec = x[n:n-M:-1] if n >= M else np.pad(x[:n+1][::-1], (0, M-n-1))
        y[n] = np.dot(w, x_vec)
        e[n] = d[n] - y[n]

        # ゲインベクトル
        Px = P @ x_vec
        k = Px / (lam + x_vec @ Px)

        # フィルタ係数更新
        w = w + k * e[n]

        # 逆相関行列更新
        P = (P - np.outer(k, x_vec @ P)) / lam

    return y, e

import matplotlib.pyplot as plt

np.random.seed(42)

# 信号生成
N = 1000
t = np.arange(N)
signal = np.sin(2 * np.pi * 0.01 * t)  # 目的信号
noise_ref = np.random.randn(N)  # 参照ノイズ

# 未知のシステムを通ったノイズ（FIRフィルタで模擬）
h_true = np.array([0.5, -0.3, 0.2, -0.1])
noise_filtered = np.convolve(noise_ref, h_true, mode='full')[:N]

# 観測信号 = 目的信号 + フィルタされたノイズ
observed = signal + noise_filtered

# 各アルゴリズムでノイズ除去
M = 8  # フィルタ次数
y_lms, e_lms, _ = lms_filter(noise_ref, observed, M, mu=0.01)
y_nlms, e_nlms = nlms_filter(noise_ref, observed, M, mu=0.5)
y_rls, e_rls = rls_filter(noise_ref, observed, M, lam=0.99)

# プロット
fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)

axes[0].plot(t, signal, 'g', label='Original signal')
axes[0].plot(t, observed, 'r', alpha=0.5, label='Observed (signal + noise)')
axes[0].set_title('Input Signals')
axes[0].legend()

axes[1].plot(t, signal, 'g', alpha=0.5, label='Original')
axes[1].plot(t, e_lms, 'b', label='LMS output')
axes[1].set_title('LMS Noise Cancellation')
axes[1].legend()

axes[2].plot(t, signal, 'g', alpha=0.5, label='Original')
axes[2].plot(t, e_nlms, 'b', label='NLMS output')
axes[2].set_title('NLMS Noise Cancellation')
axes[2].legend()

axes[3].plot(t, signal, 'g', alpha=0.5, label='Original')
axes[3].plot(t, e_rls, 'b', label='RLS output')
axes[3].set_title('RLS Noise Cancellation')
axes[3].legend()

plt.tight_layout()
plt.show()

# 学習曲線（MSEの移動平均）
window = 50
mse_lms = np.convolve(e_lms**2, np.ones(window)/window, mode='valid')
mse_nlms = np.convolve(e_nlms**2, np.ones(window)/window, mode='valid')
mse_rls = np.convolve(e_rls**2, np.ones(window)/window, mode='valid')

plt.figure(figsize=(10, 5))
plt.semilogy(mse_lms, label='LMS')
plt.semilogy(mse_nlms, label='NLMS')
plt.semilogy(mse_rls, label='RLS')
plt.xlabel('Sample')
plt.ylabel('MSE')
plt.title('Learning Curves')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()