自适应滤波器实验代码（LMS/NLMS/APA/APSA/VSS_APSA/MVSS_APSA）

效果图

 

 

完整代码

import numpy as np
import matplotlib.pyplot as plt


# Define the LMS algorithm
def LMS(x, d, mu, order):
    w = np.zeros(order)  # Initialize weights to zeros
    y = np.zeros_like(d)  # Initialize output to zeros
    e = np.zeros_like(d)
    for n in range(order, len(x)):
        x_n = x[n - order:n][::-1]  # Select current input sample
        y[n] = np.dot(w, x_n)  # Compute output using current weights
        e[n] = d[n] - y[n]  # Compute error between desired and actual output
        w += mu * e[n] * x_n  # Update weights using error and learning rate

    return y, e


def NLMS(x, d, m, order):
    w = np.zeros(order)
    y = np.zeros_like(d)
    e = np.zeros_like(d)
    for n in range(order, len(x)):
        x_n = x[n - order:n][::-1]
        y[n] = np.dot(w, x_n)
        e[n] = d[n] - y[n]
        w += m * x_n * e[n] / np.sum(x_n ** 2)
    return y, e


def APA(x: list, d: list, mu: float, M: int, L: int, delta: float) -> tuple:
    w = np.zeros(order)
    y = np.zeros_like(d)
    e = np.zeros_like(d)
    for n in range(M + L, len(x)):
        x_n = x[n - order:n][::-1]
        y[n] = np.dot(w, x_n)
        e[n] = d[n] - y[n]
        U = np.zeros((M, L))
        D = np.zeros(M)
        for m in range(M):
            U[m] = x[n - m - L:n - m][::-1]
            D[m] = d[n - m]
        U = U.T
        E = D - np.dot(U.T, w)
        w += mu * np.dot(np.dot(U, np.linalg.inv(np.dot(U.T, U) + delta * np.eye(M))), E)
    return y, e


def APSA(x: list, d: list, mu: float, M: int, L: int) -> tuple:
    w = np.zeros(order)
    y = np.zeros_like(d)
    e = np.zeros_like(d)
    for n in range(M + L, len(x)):
        x_n = x[n - order:n][::-1]
        y[n] = np.dot(w, x_n)
        e[n] = d[n] - y[n]
        U = np.zeros((M, L))
        D = np.zeros(M)
        for m in range(M):
            U[m] = x[n - m - L:n - m][::-1]
            D[m] = d[n - m]
        U = U.T
        E = D - np.dot(U.T, w)
        w += mu * np.dot(U, np.sign(E)) / np.sqrt(np.dot(np.dot(np.dot(np.sign(E.T), U.T), U), np.sign(E)))
    return y, e


def VSS_APSA(x: list, d: list, M: int, L: int, mu: float = 0.1, alpha: float = 0.96) -> tuple:
    w = np.zeros(order)
    y = np.zeros_like(d)
    e = np.zeros_like(d)
    for n in range(M + L, len(x)):
        x_n = x[n - order:n][::-1]
        y[n] = np.dot(w, x_n)
        e[n] = d[n] - y[n]
        U = np.zeros((M, L))
        D = np.zeros(M)
        for m in range(M):
            U[m] = x[n - m - L:n - m][::-1]
            D[m] = d[n - m]
        U = U.T
        E = D - np.dot(U.T, w)
        mu = alpha * mu + (1 - alpha) * np.min(
            [np.sum(np.abs(E)) / np.sqrt(np.dot(np.dot(np.dot(np.sign(E.T), U.T), U), np.sign(E))), mu])
        w += mu * np.dot(U, np.sign(E)) / np.sqrt(np.dot(np.dot(np.dot(np.sign(E.T), U.T), U), np.sign(E)))
    return y, e


def MVSS_APSA(x: list, d: list, M: int, L: int, mu: float = 0.1, alpha: float = 0.96, Nw: int = 8) -> tuple:
    w = np.zeros(order)
    y = np.zeros_like(d)
    e = np.zeros_like(d)
    R = np.zeros(L)
    C = 1.483 * (1 + 5 / (Nw - 1))
    DELTA_U_M = 0
    DELTA_E_M = 0
    for n in range(M + L, len(x)):
        x_n = x[n - order:n][::-1]
        y[n] = np.dot(w, x_n)
        e[n] = d[n] - y[n]
        U = np.zeros((M, L))
        D = np.zeros(M)
        for m in range(M):
            U[m] = x[n - m - L:n - m][::-1]
            D[m] = d[n - m]
        U = U.T
        E = D - np.dot(U.T, w)
        if n - M - L % Nw == 0 and n - M - L >= Nw:
            EA = y[n - Nw:n][::-1] - d[n - Nw:n][::-1]
            EB = (y[n - Nw:n][::-1] - d[n - Nw:n][::-1]) ** 2
            R = alpha * R + C * (1 - alpha) * U[:, 0] * np.median(EA)
            DELTA_U_M = alpha * DELTA_U_M + (1 - alpha) * x[n] ** 2
            DELTA_E_M = alpha * DELTA_E_M + C * (1 - alpha) * np.median(EB)
            DELTA_V_M = DELTA_E_M - np.dot(R.T, R) / DELTA_U_M
            mu = alpha * mu + (1 - alpha) * np.min([(np.sum(np.abs(E)) - M * np.sqrt(2 / np.pi) * np.sqrt(
                DELTA_V_M)) / np.sqrt(np.dot(np.dot(np.dot(np.sign(E.T), U.T), U), np.sign(E))), mu])
        w += mu * np.dot(U, np.sign(E)) / np.sqrt(np.dot(np.dot(np.dot(np.sign(E.T), U.T), U), np.sign(E)))
    return y, e


n = 2000
t = np.arange(n)
original = np.sin(0.1 * t)
noisy = np.random.randn(n)
x = original + 0.2 * noisy
d = np.sin(0.1 * t)

# Define the filter coefficients
order = 10

# Apply the LMS filter
y1, e1 = LMS(x, d, 0.1, order)
y2, e2 = NLMS(x, d, 0.1, order)
y3, e3 = APA(x, d, 0.1, 5, 10, 0.001)
y4, e4 = APSA(x, d, 0.1, 5, 10)
y5, e5 = VSS_APSA(x, d, 5, 10)
y6, e6 = MVSS_APSA(x, d, 5, 10)
# Plot the original and filtered signals


fig, ax = plt.subplots(9, 1, figsize=(8, 6))
ax[0].plot(t, original)
ax[0].set_title('Original Signal')
ax[1].plot(t, noisy)
ax[1].set_title('Noisy Signal')
ax[2].plot(t, x)
ax[2].set_title('Input Signal')
ax[3].plot(t, y1)
ax[3].set_title('LMS Algorithm')
ax[4].plot(t, y2)
ax[4].set_title('NLMS Algorithm')
ax[5].plot(t, y3)
ax[5].set_title('APA Algorithm')
ax[6].plot(t, y4)
ax[6].set_title('APSA Algorithm')
ax[7].plot(t, y5)
ax[7].set_title('VSS_APSA Algorithm')
ax[8].plot(t, y6)
ax[8].set_title('MVSS_APSA Algorithm')

plt.figure(2)
plt.plot([np.sum(e1[i - 20:i] ** 2) / 20 for i in range(20, len(e1), 20)], label="LMS")
plt.plot([np.sum(e2[i - 20:i] ** 2) / 20 for i in range(20, len(e2), 20)], label="NLMS")
plt.plot([np.sum(e3[i - 20:i] ** 2) / 20 for i in range(20, len(e3), 20)], label="APA")
plt.plot([np.sum(e4[i - 20:i] ** 2) / 20 for i in range(20, len(e4), 20)], label="APSA")
plt.plot([np.sum(e5[i - 20:i] ** 2) / 20 for i in range(20, len(e5), 20)], label="VSS_APSA")
plt.plot([np.sum(e6[i - 20:i] ** 2) / 20 for i in range(20, len(e6), 20)], label="MVSS_APSA")
plt.legend()
plt.show()

参考文献：

【1】变步长仿射投影符号-李雪蕊

【2】Variable Step-Size Affine Projection Sign Algorithm，Variable Step-Size Affine Projection Sign Algorithm

【3】A Recursive Least M-Estimate (RLM) Adaptive Filter for Robust Filtering in Impulse Noise Y. Zou, S. C. Chan, and T. S. Ng