自适应滤波方法—

自适应滤波方法——LMS算法

news2024/11/22 21:05:51

自适应滤波器

自适应滤波器：一种能够根据输入信号自动调整自身参数的数字滤波器
非自适应滤波器：具有静态滤波器系统的数字滤波器，静态系数构成了滤波器的传递函数
对于一些应用（如系统辨识、预测、去噪等），无法事先知道需要操作的参数，必须使用自适应的系数进行处理，这种情况下通常使用自适应滤波器
自适应滤波器处理语音信号时，不需要实现知道输入信号和噪声的统计特性，滤波器自身能够在工作过程中学习或估计信号的统计特性，并以此为依据调整自身参数，以达到某种准则/代价函数下的最优滤波效果
自适应滤波器是处理非平稳信号的一种有效手段
下面从一个N阶的滤波器出发，逐步过渡到自适应滤波器

N阶滤波器

一个N阶滤波器，参数为 $\text{w}(n)$ ，包含N个独立的滤波器系数
滤波器的输出为输入信号和滤波器系数的线性卷积：
$y(n)=\sum_{i=0}^{N-1}w_i(n)x(n-i)=\text{w}^T(n)\text{x}(n)=\text{x}^T(n)\text{w}(n)$
其中，
- $\text{x}(n)=[x(n),x(n-1),...,x(n-N+1)]^T$
- $\text{w}(n)=[w_0(n),w_1(n),...,w_{N-1}(n)]^T$
- 括号里的n表示不同的时刻
定义期望输出为d(n)，误差序列为： $e(n)=d(n)-y(n)=d(n)-\text{w}^T(n)\text{x}(n)$
由于n时刻的误差是根据n时刻的滤波器系数进行定义的，因此该误差被称为“后验误差”

标准LMS算法

根据最小化均方误差（MMSE）准则，最小化目标函数为
$J(\text{w})=E[|e(n)|^2]=E[|d(n)-\text{w}^T(n)\text{x}(n)|^2]$
为了最小化均方误差函数，需计算 $J(\text{w})$ 对w的导数，令导数为零：
$E[\text{x}(n)\text{x}^T(n)]\text{w}(n)-E[\text{x}(n)d(n)]=R\text{w}-r=0$
可得到： $\text{w}_{opt}=R^{-1}r$
称 $\text{w}_{opt}$ 定义的滤波器为Wiener滤波器，Wiener滤波器是均方误差最小意义上的统计最优滤波器
Wiener滤波器的缺陷：
- 不适用于非平稳过程，而语音就是非平稳信号
- 数学期望的计算需要利用全部的历史数据，来估计当前的滤波器系数，这不是合理的选择
由于涉及数学期望，在实际应用中，无法获得信号真实的自相关矩阵R，以及输入信号和期望输出之前的互相关向量r。因此，在实际应用中，利用瞬时误差的平方代替目标函数中的数学期望，再求取瞬时梯度： $\triangledown (n)=-2e(n)\text{x}(n)$ ，然后沿着负梯度方向更新滤波器系数：
$\text{w}(n+1)=\text{w}(n)+2\mu \text{x}(n)e(n)$
上式为标准时序LMS算法的更新公式，为逐点更新，每当给定一个新的x(n)和d(n)，滤波器系数就更新一次
标准LMS算法的执行流程：
- 初始化w(0)、x(0)
- 对每一个新的输入值x(n)，计算输出值y(n)
- 利用期望输出d(n)，计算误差值e(n)，从而得到梯度
- 更新滤波器系数，得到w(n+1)
- 返回步骤2，直至结束，可得到输出序列和误差序列
标准LMS算法
- 基本思想是梯度下降
- 优点是算法简单，易实现
- 缺点是收敛速度慢，因为存在梯度噪声，所以跟踪性能比较差

block（分块）LMS算法

用更多采样点，更新一次滤波器系数：
$\begin{aligned} \text{w}(n+1)&=\text{w}(n)+2\mu \text{x}(n)e(n) \\ \text{w}(n+2)&=\text{w}(n+1)+2\mu \text{x}(n+1)e(n+1) \\ ... \\ \text{w}(n+L)&=\text{w}(n+L-1)+2\mu \text{x}(n+L-1)e(n+L-1) \end{aligned}$
得： $\text{w}(n+L)=\text{w}(n)+2\mu \sum_{m=0}^{L-1} \text{x}(n+m)e(n+m)$
分块更新，每L点更新一次滤波器系数
注意，误差序列是通过同一个滤波器得到的： $e(n+m)=d(n+m)-\text{x}^T(n+m)\text{w}(n)$
为了后续推到方便，令L=N，也就是说收集到采样点个数等于滤波器阶数时，才进行滤波器系数更新
与标准LMS相比，block LMS的更新频率降低了L倍，因此定义一个新的时间索引k来代替n， $k=\frac{n}{L}$
从而block LMS更新公式为： $\text{w}(k+1)=\text{w}(k)+2\mu \sum_{m=0}^{L-1} \text{x}(kL+m)e(kL+m)$
block LMS中，瞬时梯度为： $\triangledown (k)=-2 \sum_{m=0}^{L-1} \text{x}(kL+m) e(kL+m)$
block LMS的核心运算：
- 输出序列为输入信号与滤波器系数的线性卷积： $y(n+m)=\text{x}^T(n+m)\text{w}(n)，m \in [0,N-1]$
- 瞬时梯度为输入信号与误差序列的线性相关： $\triangledown (k)=-2 \sum_{m=0}^{L-1} \text{x}(kL+m) e(kL+m)$
利用FFT计算线性卷积的方法有Overlap-save method和Overlap-add method

Overlap-save method

线性卷积与圆周卷积的关系：一般地，如果两个有限长序列的长度为 $N_1、N_2$ ，且满足 $N_1 \ge N_2$ ，则圆周卷积的后 $N_1-N_2+1$ 个点，与线性卷积的结果一致
线性相关与圆周相关的关系：一般地，如果两个有限长序列的长度为 $N_1、N_2$ ，且满足 $N_1 \ge N_2$ ，则圆周相关的前 $N_1-N_2+1$ 个点，与线性相关的结果一致
希望理解上述两个关系，可参考数字信号的基本运算——线性卷积（相关）和圆周卷积（相关）
为了得到N个点的输出序列： $y(n+m)=\text{x}^T(n+m)\text{w}(n)，m \in [0,N-1]$ ，要保证至少有N个点的线性卷积和圆周卷积的结果相同，即 $N_1-N_2+1 \ge N$
由于 $N_1 \ge N_2$ ，且滤波器阶数为N，因此 $N_2=N$ ，所以条件变为： $N_1 \ge 2N-1$
为了方便FFT的计算，令 $N_1=2N$
根据下图，构造x(n)和w(n)
分别计算输入信号向量和滤波器向量的傅里叶变换
$\begin{aligned} X(k)&=diag\{ F[x(kN-N), ..., x(kN-1), x(kN), ..., x(kN+N-1)] \} \\ W(k)&=F[\text{w}^T(k), 0, ..., 0]^T \\ Y(k)&=X(k)W(k) \end{aligned}$
每个Y(k)的长度为2N，输出序列y(k)的N点线性卷积等于Y(k)经过傅里叶逆变换的后N个点
$\begin{aligned} y(k)&=[y(kN), y(kN+1), y(kN), ..., y(kN+N-1)] \\ &=\text{last N samples of }F^{-1}Y(k) \end{aligned}$
y(k)的计算流程如下：
接下来计算N个点的瞬时梯度之和： $\triangledown (k)=-2 \sum_{m=0}^{L-1} \text{x}(kL+m) e(kL+m)$
对于线性相关运算，仍然是变换到频域，计算输入信号的共轭谱与误差序列的乘积
X(k)在上面已求得，接下来将误差序列 $e (k) = [e (k N), e (k N + 1), ..., e (k N + N - 1)]$ 也扩展到2N的长度
注意：计算卷积时，在w(n)的后面补N个0，计算相关时，在e(n)的前面补N个0
误差序列的傅里叶变换 $E(k)=F[0, 0, ..., 0, e^T(k)]^T$
瞬时梯度的傅里叶变换 $\bigtriangledown (k)=X^H(k)E(k)$
每个 $\bigtriangledown (k)$ 的长度为2N，瞬时梯度之和的N点线性相关等于 $\bigtriangledown (k)$ 经过傅里叶逆变换的前N个点
$\begin{aligned} \vec{\triangledown}(k)&=[\triangledown(kN), \triangledown(kN+1), \triangledown(kN), ..., \triangledown(kN+N-1)] \\ &=\text{first N samples of }F^{-1}\bigtriangledown(k) \end{aligned}$