维纳滤波(Wiener Filter)

Wiener Filter

因为最近看文章接触了维纳滤波，所以这里写一下Weiner Filter的一些简单理解和推导。

基本定义

维纳滤波是一种在含噪声的时序信号把信号提取出来的滤波器，其基本框图如下：

简单的维纳滤波其实就是通过一个FIR滤波器，去除噪声的过程。在这里，$h$的作用也可以理解为： 通过训练集的数据对信号和噪声的建模，然后通过前几个点的信息，预测当前时刻的噪声信号所占的比例，然后去除掉，剩下的就是预测的时序信号了。维纳滤波作为一种使用很广泛的滤波器，其变化的形式也有很多种，可以是单输入输出的，也可以是多输入输出的。$h$所表示的变换也可以写成非线性；$h$可以是有限长的FIR滤波器，也可以是无限长的IIR滤波器。要取决于当前你所解决的问题。但是维纳滤波的基本思想还是一致的。通过滤波（矩阵或者其他模型的形式）来从信号和噪声的混合中提取信号。所以维纳滤波的核心，就是计算这个滤波器（矩阵$h$或者模型的参数）。也就是解Wiener-Hopf方程。

本文用比较简单的单输入输出，且只考虑有限长滤波（即认为当前时刻的信号只和前有限个时间点的信号相关）。

公式推导

首先，对于图1中的滤波器：

$$y(n)=x(n)\ast h(n)=(s(n)+v(n))\ast h(n)\text{\hspace{1em}(1)}$$

其中$\ast$表示卷积，$x(n)$表示输入信号， $y(n)$表示输出信号， $s(n)$表示输入信号中，有用的信号部分；$v(n)$表示输入信号中的噪声部分。

维纳滤波的目标是，保证输出$y(n)$和真实信号$s(n)$的差别最小，由于$y(n)$和$s(n)$是时序信号，所以要保证两者的均方误差最小，所以有：

E{e2(n)}=E{(y(n)−s(n))2}=E{(x(n)∗h(n)−s(n))2}(2)E\{e^2(n)\} = E\{(y(n)-s(n))^2\} = E\{(x(n)*h(n)-s(n))^2\} \tag{2} E{e2(n)}=E{(y(n)−s(n))2}=E{(x(n)∗h(n)−s(n))2}(2)

即求使得Eq(2)最小的$h$。所以E{e2}E\{e^2\}E{e2}对$h$求偏导。有：

∂E{e2(n)}∂h=2E{e(n)∗∂e(n)∂h}=0(3)\frac{\partial{E\{e^2(n)\}}}{\partial{h}} = 2E\{e(n) * \frac{\partial{e(n)}}{\partial{h}}\} = 0 \tag{3} ∂h∂E{e2(n)}​=2E{e(n)∗∂h∂e(n)​}=0(3)

∂E{e2(n)}∂h=2E{[∑m=0N−1h(m)x(n−m)−s(n)]x(n−j)},j=0,1,...,N−1(4)\frac{\partial{E\{e^2(n)\}}}{\partial{h}} = 2E\{[\sum_{m=0}^{N-1}{h(m)x(n-m) - s(n)}]x(n-j)\}, j = 0, 1, ... , N-1 \tag{4} ∂h∂E{e2(n)}​=2E{[m=0∑N−1​h(m)x(n−m)−s(n)]x(n−j)},j=0,1,...,N−1(4)

∂E{e2(n)}∂h=2∑m=1N−1h(m)E{x(n−j)x(n−m)}−2E{s(n)x(n−j)}=0,j=0,1,...,N−1(5)\frac{\partial{E\{e^2(n)\}}}{\partial{h}} = 2\sum_{m=1}^{N-1}{h(m)}E\{x(n-j)x(n-m)\} - 2E\{s(n)x(n-j)\} = 0, j = 0, 1, ..., N-1 \tag{5} ∂h∂E{e2(n)}​=2m=1∑N−1​h(m)E{x(n−j)x(n−m)}−2E{s(n)x(n−j)}=0,j=0,1,...,N−1(5)

我们设$x$和$s$的相关系数为$R_{xs}$，则有：

$$R_{xs}(j)=\sum _{m=0}^{N-1}h(m)R_{xx}(j-m),j=0,1,...,N-1\text{\hspace{1em}(6)}$$

其中，$R_{xx}(j-m)$表示$x(n-j)$和$x(n-m)$的相关系数，这里$m$是固定的，$j$是变化的。且$m>=0$，$R_{xs}(j)$表示$x(n-j)$和$s(n)$的相关系数。上述公式中，$n$表示的是时序信号中的时间点。

然后，根据Eq(6)，可以得到$N$个线性方程：
$$\{\begin{array}{l}{R}_{xs}(0)=h(0)R_{xx}(0)+h(1)R_{xx}(1)+...+h(N-1)R_{xx}(N-1)\\ R_{xs}(1)=h(1)R_{xx}(1)+h(0)R_{xx}(0)+...+h(N-1)R_{xx}(N-2)\\ ...\\ R_{xs}(N-1)=h(N-1)R_{xx}(N-1)+h(N-2)R_{xx}(N-2)+...+h(0)R_{xx}(0)\end{array}\text{\hspace{1em}(7)}$$
写成矩阵形式，有：

$${\mathbf{R}}_{\mathbf{xx}}\mathbf{H}={\mathbf{R}}_{\mathbf{xs}}\text{\hspace{1em}(8)}$$

其中， $\mathbf{H}=[h(0),h(1),...,h(N-1){]}^T$是需要求的滤波器参数

$${\mathbf{R}}_{\mathbf{xs}}=[R_{xs}(0),R_{xs}(1),...,R_{xs}(N-1){]}^T$$是$x$和$s$的相关系数
$${\mathbf{R}}_{\mathbf{xx}}=\left[\begin{array}{ccccc}{R}_{xx}(0)& R_{xx}(1)& ...& R_{xx}(N-1)& \\ R_{xx}(1)& R_{xx}(0)& ...& R_{xx}(N-2)& \\ ⋮& ⋮& ...& ⋮& \\ R_{xx}(N-1)& R_{xx}(N-2)& ...& R_{xx}(0)& \end{array}\right]\text{\hspace{1em}(9)}$$

所以根据Eq(8)可以求得：

$$\mathbf{H}={\mathbf{R}}_{\mathbf{xx}}^{-1}{\mathbf{R}}_{\mathbf{xs}}\text{\hspace{1em}(10)}$$

此时，信号的均方误差最小，根据Eq(2)，可得：

E{e2(n)}=E{(s(n)−∑m=0N−1h(m)x(n−m))2}(11)E\{e^2(n)\} = E\{(s(n)-\sum_{m=0}^{N-1}h(m)x(n-m))^2\} \tag{11}E{e2(n)}=E{(s(n)−m=0∑N−1​h(m)x(n−m))2}(11)

E{e2(n)}=E{s2(n)−2s(n)∑m=0N−1h(m)x(n−m)+∑m=0N−1∑r=0N−1h(m)x(n−m)h(r)x(n−r)}E\{e^2(n)\} = E\{s^2(n) - 2s(n)\sum_{m=0}^{N-1}h(m)x(n-m)+\sum_{m=0}^{N-1}\sum_{r=0}^{N-1}{h(m)x(n-m)h(r)x(n-r)}\}E{e2(n)}=E{s2(n)−2s(n)m=0∑N−1​h(m)x(n−m)+m=0∑N−1​r=0∑N−1​h(m)x(n−m)h(r)x(n−r)}

E{e2(n)}=Rss(0)−2∑m=0N−1h(m)Rxs(m)+∑m=0N−1h(m)∑r=0N−1h(r)Rxx(m−r)E\{e^2(n)\}=R_{ss}(0)-2\sum_{m=0}^{N-1}{h(m)R_{xs}(m)+\sum_{m=0}^{N-1}{h(m)}\sum_{r=0}^{N-1}{h(r)R_{xx}(m-r)}}E{e2(n)}=Rss​(0)−2m=0∑N−1​h(m)Rxs​(m)+m=0∑N−1​h(m)r=0∑N−1​h(r)Rxx​(m−r)

根据Eq(5)，可得：

E{e2(n)}=Rss(0)−∑m=0N−1h(m)Rxs(m)(12)E\{e^2(n)\} = R_{ss}(0) - \sum_{m=0}^{N-1}{h(m)R_{xs}(m)} \tag{12}E{e2(n)}=Rss​(0)−m=0∑N−1​h(m)Rxs​(m)(12)

假设信号$s(n)$和噪声$v(n)$互相独立，那么有：

$$R_{sv}=R_{vs}=0$$

$$R_{xs}=R_{ss}+R_{vs}=R_{ss}$$

$$R_{xx}=R_{ss}+R_{sv}+R_{vs}+R_{vv}=R_{ss}+R_{vv}$$

则，根据Eq(12)，有：

E{e2(n)}=Rss(0)−∑m=0N−1h(m)Rss(m)(14)E\{e^2(n)\} = R_{ss}(0) - \sum_{m=0}^{N-1}{h(m)R_{ss}(m)} \tag{14}E{e2(n)}=Rss​(0)−m=0∑N−1​h(m)Rss​(m)(14)

至此，最简单的维纳滤波的基本公式推导完成，如果涉及到多输入多输出的维纳滤波，会更加复杂，这里不做推导。