Question

本章将简要介绍自适应滤波器的原理以及其最常用的算法 NLMS，并给 NLMS 算法的两种实现方法：用纯 Python 编写，和用 ctypes 调用 C 语言编写。最后将对 NLMS 算法进行一些的实验。

自适应滤波器简介

近年来，随着数字信号处理器的功能的不断增强，自适应信号处理 (adaptive signal process) 活跃在噪声消除、回声控制、信号预测、声音定位等众多信号处理领域。

尽管其应用领域十分广泛，但基本的系统构造大致只有如下几种分类。

系统辨识

所谓系统辨识(system identification)，就是通过对未知系统的输入输出进行观测，构造一个滤波器使得它在同样的输入的情况下，输出信号和未知系统相同。简而言之，就是通过观测未知系统对输入的反应，探知其内部情况。为了探知内情而使用的输入信号我们称之为参照信号。

系统识别(System Identification) 的框图

如上图所示参照信号 x(j) 同时输入到未知系统和自适应滤波器 H 中，未知系统的输出为 y(j)， 自适应滤波器的输出为 u(j)，由于观测误差或者外部噪声的干扰，实际观测到的未知系统的输出为 d(j)=y(j)+n(j)，n(j) 被称为外部干扰。通过求的 d(j) 和 u(j) 之间的误差 e(j)=d(j)-u(j)，我们可以知道自适应滤波器 H 和未知系统还有多少差别，通过这个误差我们更新 H 的内部参数，使得它更加靠近未知系统。

上面各个公式中的 j 表示某一时刻，因为我们讨论的是数字信号处理，已经对所有的信号进行取样，因此可以把 j 简单的看作取样点的下标。

信号预测

所谓信号预测就是通过信号过去的值预测（计算）现在的值，下面是信号预测的系统框图。

信号预测(Predication) 框图

x(j) 是待预测的信号，假设我们无法完美地观测此信号，因此导入一个外部干扰 n(j)，这样 d(j)=x(j)+n(j) 就是我们观测到的待预测信号。

通过延迟器将 d(j) 进行延时得到 d(j-D)，并把 d(j-D) 输入到自适应滤波器 H 中，得到其输出为 u(j)，u(j) 就是自适应滤波器通过待预测信号过去的值预测出的现在的值，计算观测值 d(j) 和预测值 u(j) 之间的误差 e(j)=d(j)-u(j)，通过 e(j) 更新自适应滤波器 H 的内部系数使得其输出更加接近 d(j)。

如果 x(j) 存在白色噪声的成分和周期信号的成分，由于白色噪声是完全不自相关，无法预测的信号，因此通过过去的值 x(j-D) 所能预测的只能是其中的周期信号的成分。这样自适应滤波器 H 的输出信号 u(j) 就会与周期信号成分渐渐逼近，而 e(j) 则是剩下的不可预测的白色噪声的成分。因此自适应滤波器也可以运用于噪声消除。

信号均衡

信号均衡(Equalization) 框图

当信号 x(j) 通过未知系统之后变成 y(j)，未知系统对信号 x(j) 进行了某种改变，使得其波形产生歪曲。我们希望均衡器矫正这种歪曲，也就是通过 y(j) 重建原始信号 x(j)，由于因果律还原原始信号 x(j) 是不可能的，我们只能还原其延时了的信号 x(j-D)。x(j) 和 x(j-D) 除了时间上的延迟之外，其它特性完全相同。

这里我们将观测到的未知系统的输出 y(j)+n(j) 输入到自适应滤波器 H 中，通过 H 的系数更新使得其输出 u(j) 逐渐逼近原始信号的延时 x(j-D)。这样我们就构建了一个滤波器 H 使得它与未知系统的卷积正好等于一个脉冲传递函数。也就是说 H 的频域特性恰好能抵消未知系统的所带来的改变。

NLMS 计算公式

自适应滤波器中最重要的一个环节就是其系数的更新算法，如果不对自适应滤波器的系数更新的话，那么它就只是一个普通的 FIR 滤波器了。系数更新算法有很多种类，最基本、常用、简单的一种方法叫做 NLMS(归一化最小均方)，让我们先来看看它的数学公式表达：

设置自适应滤波器系数 的所有初始值为 0, 的长度为 I。

对每个取样值进行如下计算，其中 n=0, 1, 2, ...

自适应滤波器系数 是一个长度为 I 的矢量，也就是一个长度为 I 的 FIR 滤波器。在时刻 n，滤波器的每个系数对应的输入信号为 ，它也是一个长度为 I 的矢量。这两个矢量的点乘即为滤波器的输出和目标信号 d(n) 之间的差为 e(n)，然后根据 e(n) 和 ， 更新滤波器的系数。

数学公式总是令人难以理解的，下面我们以图示为例进行说明：

NLMS 算法示意图

图中假设自适应滤波器 h 的长度为 4，在时刻 7 滤波器的输出为:

u[7] = h[0]*x[7] + h[1]*x[6] + h[2]*x[5] + h[3]*x[4]

滤波器的输入信号的平方和 powerX 为:

powerX = x[4]*x[4] + x[5]*x[5] + x[6]*x[6] + x[7]*x[7]

未知系统的输出 d[7]和滤波器的输出 u[7]之间的差为:

e[7] = d[7] - u[7]

使用 u[7]和 x[4]..x[7]对滤波器的系数更新:

h[4] = h[4] + u * e[7]*x[4]/powerX
h[4] = h[5] + u * e[7]*x[5]/powerX
h[4] = h[6] + u * e[7]*x[6]/powerX
h[4] = h[7] + u * e[7]*x[7]/powerX

其中参数 u 成为更新系数，为 0 到 1 之间的一个实数，此值越大系数更新的速度越快。对于每个时刻 i 都需要进行上述的计算，因此滤波器的系数对于每个参照信号 x 的取样都更新一次。

NumPy 实现

按照上面介绍的 NLMS 算法，我们很容易写出用 NumPy 实现的 NLMS 计算程序：

# -*- coding: utf-8 -*-
# filename: nlms_numpy.py

import numpy as np

# 用 Numpy 实现的 NLMS 算法
# x 为参照信号，d 为目标信号，h 为自适应滤波器的初值
# step_size 为更新系数
def nlms(x, d, h, step_size=0.5):
   i = len(h)
   size = len(x)
   # 计算输入到 h 中的参照信号的乘方 he
   power = np.sum( x[i:i-len(h):-1] * x[i:i-len(h):-1] )
   u = np.zeros(size, dtype=np.float64)

   while True:
     x_input = x[i:i-len(h):-1]
     u[i] = np.dot(x_input , h)
     e = d[i] - u[i]
     h += step_size * e / power * x_input

     power -= x_input[-1] * x_input[-1] # 减去最早的取样
     i+=1
     if i >= size: return u
     power += x[i] * x[i] # 增加最新的取样

为了节省计算时间，我们用一个临时变量 power 保存输入到滤波器 h 中的参照信号 x 的能量。在对于 x 中的每个取样的循环中，power 减去 x 中最早的一个取样值的乘方，增加最新的一个取样值的乘方。这样为了计算参照信号的能量，每次循环只需要计算两次乘法和两次加法即可。

nlms 函数的输入为参照信号 x、目标信号 d 和自适应滤波器的系数 h。因为在后面的模拟计算中，d 是 x 和未知系统的脉冲响应的卷积而计算的来，它的长度会大于 x 的参数，因此循环体的循环次数以参照信号的长度为基准。

为了对自适应滤波器的各种应用进行模拟，我们还需要如下的几个辅助函数，完整的程序请参考 NLMS 算法的模拟测试 。

def make_path(delay, length):
   path_length = length - delay
   h = np.zeros(length, np.float64)
   h[delay:] = np.random.standard_normal(path_length) * np.exp( np.linspace(0, -4, path_length) )
   h /= np.sqrt(np.sum(h*h))
   return h

make_path 产生一个长度为 length，最小延时为 delay 的指数衰减的波形。这种波形和封闭空间的声音的传递函数有些类似之处，因此在计算机上进行声音的算法模拟时经常用这种波形作为系统的传递函数。:

def plot_converge(y, u, label=""):
  size = len(u)
  avg_number = 200
  e = np.power(y[:size] - u, 2)
  tmp = e[:int(size/avg_number)*avg_number]
  tmp.shape = -1, avg_number
  avg = np.average( tmp, axis=1 )
  pl.plot(np.linspace(0, size, len(avg)), 10*np.log10(avg), linewidth=2.0, label=label)

def diff_db(h0, h):
   return 10*np.log10(np.sum((h0-h)*(h0-h)) / np.sum(h0*h0))

plot_converge 绘制信号 y 和信号 u 之间的误差，每 avg_number 个取样点就上一次误差的乘方的平均值。我们将用 plot_converge 函数绘制未知系统的输出 y 和自适应滤波器的输出 u 之间的误差。观察自适应滤波器是如何收敛的，以评价自适应滤波器的收敛特性。diff_db 函数同样是用来评价自适应滤波器的收敛特性，不过他是直接计算未知系统的传递函数 h0 和自适应滤波器的传递函数 h 之间的误差。下面我们会看到这两个函数得到的收敛值是相同的。

系统辨识模拟

我们用下面的函数调用 nlms 算法对系统辨识应用进行模拟：

def sim_system_identify(nlms, x, h0, step_size, noise_scale):
     y = np.convolve(x, h0)
     d = y + np.random.standard_normal(len(y)) * noise_scale # 添加白色噪声的外部干扰
     h = np.zeros(len(h0), np.float64) # 自适应滤波器的长度和未知系统长度相同，初始值为 0
     u = nlms( x, d, h, step_size )
     return y, u, h

系统识别(System Identification) 的框图

此函数的参数分别为：

• nlms : nlms 算法的实现函数
• x : 参照信号
• h0 : 未知系统的传递函数，虽然是未知系统，但是计算机模拟时它是已知的
• step_size : nlms 算法的更新系数
• noise_scale : 外部干扰的系数，此系数决定外部干扰的大小，0 表示没有外部干扰

函数的返回值分别为：

• y : 未知系统的输出，不包括外部干扰
• u : 自适应滤波器的输出
• h : 自适应滤波器的最终的系数

最后我们用下面的函数创建未知系统 h0， 参照信号 x，然后调用 sim_system_identify 函数得到结果并且绘图：

def system_identify_test1():
  h0 = make_path(32, 256) # 随机产生一个未知系统的传递函数
  x = np.random.standard_normal(10000)  # 参照信号为白噪声
  y, u, h = sim_system_identify(nlms_numpy.nlms, x, h0, 0.5, 0.1)
  print diff_db(h0, h)
  pl.figure( figsize=(8, 6) )
  pl.subplot(211)
  pl.subplots_adjust(hspace=0.4)
  pl.plot(h0, c="r")
  pl.plot(h, c="b")
  pl.title(u"未知系统和收敛后的滤波器的系数比较")
  pl.subplot(212)
  plot_converge(y, u)
  pl.title(u"自适应滤波器收敛特性")
  pl.xlabel("Iterations (samples)")
  pl.ylabel("Converge Level (dB)")
  pl.show()

自适应滤波器收敛之后的系数和收敛速度

上部的图显示的是未知系统（红色) 和自适应滤波器（蓝色) 的传递函数的系数，我们看到自适应滤波器已经十分接近未知系统了。diff_db(h0, h) 的输出为-25.35dB。下部的图通过绘制 y 和 u 之间的误差，显示了自适应滤波器的收敛过程。我们看到经过约 3000 点的计算之后，收敛过程已经饱和，最终的误差为-25dB 左右，和 diff_db 计算的结果一致。

从图中可以看到收敛过程的两个重要特性：收敛时间和收敛精度。参照信号的特性、外部干扰的大小和更新系数都会影响这两个特性。下面让我们看看参照信号为白色噪声、外部干扰的能量固定时，更新系数对它们影响：

def system_identify_test2():
  h0 = make_path(32, 256) # 随机产生一个未知系统的传递函数
  x = np.random.standard_normal(20000)  # 参照信号为白噪声
  pl.figure(figsize=(8,4))
  for step_size in np.arange(0.1, 1.0, 0.2):
    y, u, h = sim_system_identify(nlms_numpy.nlms, x, h0, step_size, 0.1)
    plot_converge(y, u, label=u"μ=%s" % step_size)
  pl.title(u"更新系数和收敛特性的关系")
  pl.xlabel("Iterations (samples)")
  pl.ylabel("Converge Level (dB)")
  pl.legend()
  pl.show()

更新系数和收敛速度的关系

下面是更新系数固定，外部干扰能量变化时的收敛特性：

def system_identify_test3():
  h0 = make_path(32, 256) # 随机产生一个未知系统的传递函数
  x = np.random.standard_normal(20000)  # 参照信号为白噪声
  pl.figure(figsize=(8,4))
  for noise_scale in [0.05, 0.1, 0.2, 0.4, 0.8]:
    y, u, h = sim_system_identify(nlms_numpy.nlms, x, h0, 0.5, noise_scale)
    plot_converge(y, u, label=u"noise=%s" % noise_scale)
  pl.title(u"外部干扰和收敛特性的关系")
  pl.xlabel("Iterations (samples)")
  pl.ylabel("Converge Level (dB)")
  pl.legend()
  pl.show()

外部干扰噪声和收敛速度的关系

从上面的图可以看出，当外部干扰的振幅增加一倍、能能量增加 6dB 时，收敛精度降低 6dB。而由于更新系数相同，所以收敛过程中的收敛速度都是一样的。

信号均衡模拟

对于信号均衡的应用我们用如下的程序进行模拟：

def sim_signal_equation(nlms, x, h0, D, step_size, noise_scale):
  d = x[:-D]
  x = x[D:]
  y = np.convolve(x, h0)[:len(x)]
  h = np.zeros(2*len(h0)+2*D, np.float64)
  y += np.random.standard_normal(len(y)) * noise_scale
  u = nlms(y, d, h, step_size)
  return h

信号均衡(Equalization) 框图

sim_signal_equation 函数的参数：

• nlms : nlms 算法的实现函数
• x : 未知系统的输入信号
• h0 : 未知系统的传递函数
• D : 延迟器的延时参数
• step_size : nlms 算法的更新系数
• noise_scale : 外部干扰的系数，此系数决定外部干扰的大小，0 表示没有外部干扰

在函数中的各个局部变量：

• d : 输入信号经过延迟器之后的信号
• y : 未知系统的输出
• h : 自适应滤波器的系数，它的长度要足够长，程序中使用 2 倍延时 + 2 倍未知系统的传递函数的长度

函数的返回值为自适应滤波器收敛后的系数，它能够均衡 h0 对输入信号所造成的影响。我们通过下面的函数产生数据、调用模拟函数以及绘制结果：

def signal_equation_test1():
  h0 = make_path(5, 64)
  D = 128
  length = 20000
  data = np.random.standard_normal(length+D)
  h = sim_signal_equation(nlms_numpy.nlms, data, h0, D, 0.5, 0.1)
  pl.figure(figsize=(8,4))
  pl.plot(h0, label=u"未知系统")
  pl.plot(h, label=u"自适应滤波器")
  pl.plot(np.convolve(h0, h), label=u"二者卷积")
  pl.title(u"信号均衡演示")
  pl.legend()
  w0, H0 = scipy.signal.freqz(h0, worN = 1000)
  w, H = scipy.signal.freqz(h, worN = 1000)
  pl.figure(figsize=(8,4))
  pl.plot(w0, 20*np.log10(np.abs(H0)), w, 20*np.log10(np.abs(H)))
  pl.title(u"未知系统和自适应滤波器的振幅特性")
  pl.xlabel(u"圆频率")
  pl.ylabel(u"振幅(dB)")
  pl.show()

如果延迟器的延时 D 不够的话，会由于因果律使得自适应滤波器无法收敛。因此这里我们采用的 D 的长度为 h0 的长度的 2 倍。下图显示 h0, h 和它们的卷积。我们看到 h0 和 h 的卷积正好是一个脉冲，其延时为正好等于 D(128)。

未知系统和自适应滤波器的级联（卷积) 近似为标准延迟

下图显示未知系统的频率响应（蓝色) 和自适应滤波器的频率响应（绿色)，我们看到二者正好相反，也就是说自适应滤波器均衡了未知系统对信号的影响。

未知系统和自适应滤波器的频率响应正好相反

卷积逆运算

虽然卷积运算最终能归结为简单的加法和乘法运算，然而卷积的逆运算就不是很容易计算了。我们知道两个线性系统 h1 和 h2 的级联 h3 可以用它们的脉冲响应的卷积计算求得，而所谓卷积的逆运算可以想象为已知 h3 和 h1，求一个 h2 使它和 h1 级联之后正好等于 h3。

根据卷积的计算公式可知，如果 h1 的长度为 100，h3 的长度为 199，那么 h2 的长度则为 100，因为 h2 的每个系数都是未知的，于是就有 100 个未知数，而这 100 个未知数需要满足 199 个线性方程：h3 中的每个系数都有一个方程与之对应。由于方程数大于未知数的个数，显然对于任意的 h1 和 h3 并不能保证有一个 h2 使得它和 h1 的卷积正好等于 h3。

既然不能精确求解，那么卷积的逆运算就变成了一个误差最小化的优化问题。用自适应滤波器计算卷积的逆运算和计算信号均衡类似，将白色噪声 x 输入到 h1 中得到信号 u，将 x 输入到 h3 中得到信号 d，然后使用 u 作为参照信号，d 作为目标信号进行 NLMS 计算，最终收敛后的自适应滤波器的系数就是 h2。

下面的程序模拟这一过程：

# -*- coding: utf-8 -*-
import numpy as np
import pylab as pl
from nlms_numpy import nlms
import scipy.signal as signal

def inv_convolve(h1, h3, length):
  x = np.random.standard_normal(10000)
  u = signal.lfilter(h1, 1, x)
  d = signal.lfilter(h3, 1, x)
  h = np.zeros(length, np.float64)
  nlms(u, d, h, 0.1)
  return h

h1 = np.fromfile("h1.txt", sep="\n")
h1 /= np.max(h1)
h3 = np.fromfile("h3.txt", sep="\n")
h3 /= np.max(h3)

pl.rc('legend', fontsize=10)
pl.subplot(411)
pl.plot(h3, label="h3")
pl.plot(h1, label="h1")
pl.legend()
pl.gca().set_yticklabels([]) 
for idx, length in enumerate([128, 256, 512]):
  pl.subplot(412+idx)
  h2 = inv_convolve(h1, h3, length)
  pl.plot(np.convolve(h1, h2)[:len(h3)], label="h1*h2(%s)" % length)
  pl.legend()
  pl.gca().set_yticklabels([]) 
  pl.gca().set_xticklabels([]) 

pl.show()

下面是程序的计算结果：

卷积逆运算演示

程序中的 h1 和 h3 从文本文件中读取而得，它们是 ANC(能动噪声控制) 系统中实际测量的脉冲响应。如果能找到一个 h2 满足卷积条件的话，就能够有效的进行噪声控制。

程序计算出 h2 的长度分别为 128, 256, 512 时的结果，可以看出 h2 越长结果越精确。

DLL 函数的编写

ctypes 的 python 接口