维纳滤波（1）

1. 概述

当系统中的有效信号和噪声都是随机过程，信号和噪声的频谱还可能重叠（比如有效信号是高斯-马尔可夫过程，噪声是白噪声），根据频域参数设计滤波器的方法就不再适用。

维纳滤波器可以在一些场合解决上述为题，其设计原则是均方误差（的期望）最小。我们从相对简单的单参数滤波器开始。

2. 参数滤波器

设输入信号为$x(t)+n(t)$，其中$n(t)$为噪声，系统冲击响应为$g(t)$，输出为$y(t)$。

将输入输出写成Laplace形式：

$Y(s)=G(s)[X(s)+N(s)]$         (2.1)

误差为有效信号与滤波器输出的差：

$e(t)=x(t)-y(t)$                     (2.2)

$E(s)=X(s)-Y(s)$                 (2.3)

(2.1)式代入(2.3)式：

$E(s)=X(s)-G(s)[X(s)+N(s)]=[1-G(s)]X(s)-G(s)N(s)$                  (2.4)

从上式可以看出，误差有两个来源：一是输入信号被传递函数“编辑”后与原始信号的差；二是系统处理后的噪声。

更进一步的说，误差第一项可看做$X(s)$通过系统$1-G(s)$，第二项可看做$N(s)$通过系统$G(s)$。

如果信号与噪声是不相关的，误差的功率谱就为$[1-G(s)][1-G(-s)]S_x(s)+G(s)G(-s)S_n(s)$。

将该功率谱作逆变换得到自相关，并对自相关取$ au=0$，均方误差就可以写为如下形式：

$E[e^2 ]=frac{1}{2pi j}int_{-jinfty}^{+jinfty}[1-G(s)][1-G(-s)]S_x(s)ds+frac{1}{2pi j}int_{-jinfty}^{+jinfty}G(s)G(-s)S_n(s)ds$        (2.5)

其中，$S_x(s)$是有效信号的功率谱，$S_n(s)$是噪声的功率谱。

在设计该种滤波器时，一般方法是使用带参的传递函数，于是(2.5)式也是一个带参的式子。针对具体问题，将均方误差对该参数求导，就可以得出满足最小均方误差条件的参数值。

上面是从频域求解最小均方误差；很快我们将看到如何从时域求解最小均方误差。

3. 稳态条件下的维纳滤波

我们先做以下几个假定：

1) 滤波器输入为信号和噪声的线性叠加，两者均为协方差平稳过程，且自相关和互相关已知；

2) 滤波器是线性时不变的；

3) 输出是协方差平稳的（即不考虑系统初启动时的状态，而认为系统已经稳定运行了较长时间）；

4) 误差记为：

     $e(t)=x(t+alpha)-y(t)$.                         (3.1)

     若$alpha >0$，这就是一个预测。

这次我们从时域来求解最小均方误差。

$e^2(t)=x^2(t+alpha)-2x(t+alpha)y(t)+y^2(t)$                   (3.2)

$y(t)$可写为卷积积分：

$y(t) = int_{-infty}^{+infty} g(u)[x(t-u)+n(t-u)] du$             (3.3)

将3.3代入3.2并求期望：

$E[e^2]=E left[  x^2(t+alpha) -2x(t+alpha) int_{-infty}^{+infty}g(u)[x(t-u)+n(t-u)]du + int_{-infty}^{+infty}g(u)[x(t-u)+n(t-u)]du int_{-infty}^{+infty}g(v)[x(t-v)+n(t-v)]dv    ight]$

 $=E[x^2(t+a)]  -2 int_{-infty}^{+infty}g(u) E { [x(t-u)+n(t-u)] x(t+alpha) } du + int_{-infty}^{+infty} int_{-infty}^{+infty} g(u)g(v) E {  [ x(t-u) + n(t-u) ] [x(t-v) + n(t-v) ] } dudv $

$=R_x(0) - 2 int_{-infty}^{+infty} g(u) R_{x+n, x}(alpha+u) du + int_{-infty}^{+infty} int_{-infty}^{+infty} g(u) g(v) R_{x+n} (u-v) du dv $    (3.4)

* 如果信号与噪声不相关，上式可以根据下面的式子进一步简化：

$R_{x+n}=R_x+R_n, R_{x+n, x} = R_x$      (3.5)

我们希望找到能够最小化$E[e^2]$的$g(u)$。这是一个典型的变分问题。因此将(3.4)写成(3.6)：

$E[e^2] = R_x(0) - 2 int_{-infty}^{+infty} [g(u) + epsilon eta(u)] R_{x+n, x}(alpha + u) du + int_{-infty}^{+infty}  int_{-infty}^{+infty}  [g(u) + epsilon eta(u)][ g(v) + epsilon eta(v)] R_{x+n}(u - v) du dv $

$frac{d E[e^2]} {depsilon} =  -2 int_{-infty}^ {+infty} eta(u)R_{x+n,x}(alpha + u) du + int_{-infty}^{+infty}  int_{-infty}^{+infty}  [g(u) eta(v) + g(v) eta(u) + 2epsilon eta(u)eta(v)] R_{x+n}(u-v) dudv $

$frac{d E[e^2]} {depsilon} | _{epsilon=0} = -2 int_{-infty}^ {+infty} eta(u)R_{x+n,x}(alpha + u) du + int_{-infty}^{+infty}  int_{-infty}^{+infty}  [g(u)eta(v) + g(v) eta(u)] R_{x+n}(u-v) dudv = 0$            (3.6)

上面第二个积分（二重积分）式子可以写成$int_{-infty}^{+infty}  int_{-infty}^{+infty} g(u)eta(v)R_{x+n}(u-v)dudv +  int_{-infty}^{+infty}  int_{-infty}^{+infty} g(v)eta(u)R_{x+n}(v-u)dvdu = 2 int_{-infty}^{+infty}  int_{-infty}^{+infty} g(v)eta(u)R_{x+n}(v-u)dvdu $

于是就得到（作一下简单的变量替换：$v ightarrow u, u ightarrow au$）

$ int_{-infty}^{+infty} eta  ( au) [-R_{x+n,x}(alpha + au) + int_{-infty}^{+infty} g(u)R_{x+n}(u- au)du]d au=0$          (3.7)

对上面这个方程，我们将首先求非因果解，然后求因果条件下的解。

3.1 非因果解

非因果解主要用于离线处理数据的情况。

若式(3.7)成立，则根据变分法基本引理：

$int_{-infty}^{+infty} g(u) R_{x+n}(u- au) du = R_{x+n,x}(alpha + au) $            (3.8)

因为自相关是偶对称的，所以上面等式的左边可看做卷积。两边对$ au$作Laplace变换：

$G(s)S_{x+n}(s) = S_{x+n,x}(s) e^{alpha s}$            (3.9)

或写为

$G(s) = frac {S_{x + n, x}(s) e^{alpha s} } { S_{x+n}(s) } $        (3.10)

这样我们就得到了非因果条件下的传递函数。

*上面的传递函数可以有一个直观的理解：

假定信号与噪声不相关，且$alpha=0$，那么：

$G(s)=frac{S_x(s)}{S_x(s)+S_n(s)}$

分子是有效信号功率谱，分母是信号+噪声的功率谱，这几乎是一个不能再直观的“功率谱误差最小”的式子了。

均方误差为（对(3.4)作一个小小的变形）：

$E[e^2] = R_x(0) - int_{-infty}^{+infty} g(u) R_{x+n,x} (alpha + u) du + int_{-infty}^{+infty} g(u) [-R_{x+n, x} (alpha + u) + int_{-infty}^{+infty} g(v) R_{x+n} (v-u) dv] du $           (3.11)

根据(3.8)，上面中括号内的式子等于零，所以：

$E[e^2] = R_x(0) - int_{-infty}^{+infty} g(u) R_{x+n, x}(alpha + u) du$

3.2 因果解

我们依然从(3.7)开始。

对因果解，可得(3.18)：

$int_{-infty}^{+infty}  g(u) R_{x+n}(u- au)  du  -  R_{x+n,x}(alpha + au) = 0$.        $ au geq 0 $            (3.18)

上式称为Wiener-Hopf（维纳－霍普夫）方程，它只需考虑$ au geq 0$的情况。

为了解这个方程，将右边的0替换为一个形式未知的函数$a( au)$，其在负轴取值不确定，在正轴取值为0。于是将(3.18)改写为：

$int_{-infty}^{+infty}  g(u) R_{x+n}(u- au)  du  -  R_{x+n,x}(alpha + au) = a( au)$.        $-infty lt au lt +infty $            (3.19)

两边做Laplace变换：

$G(s)S_{x+n}(s) - S_{x+n,x}(s) e ^{alpha s} = A(s)$                 (3.20)

将$S_{x+n}$分解为两个分式的乘积，这两个分式分别包含左半平面和右半平面的零极点：

$ [ G(s) S_{x+n}^{+} (s) ] S_{x+n}^{-} (s)   - S_{x+n,x} (s) e ^{alpha s} = A(s)  $

$G(s) S_{x+n}^+(s) = frac{A(s)}{S_{x+n}^- (s) } + frac{S_{x+n,x} e^{alpha s}} {S_{x+n}^- (s)} $            (3.21)

$S^+$包含左半平面零极点，$S^-$包含右半平面零极点。

因为$g(u)$是因果的，所以$G(s)S_{x+n}^+(s)$的零极点都在左半平面，而$a( au)$又是一个非因果信号。于是(3.21)可用如下框图表示：

[正时间函数] = [负时间函数] + [正负时间函数]

在(3.21)中令两边正时间函数相等：

$G(s)S_{x+n}^+(s)=$ positive-time part of $frac {S_{x+n,x}(s)e^{alpha s}} {S_{x+n}^-(s)}$

即

$G(s)=frac{1}{S_{x+n}^+}$  $ig[$ positive-time part of  $frac{S_{x+n,x}(s) e^{alpha s} } {S_{x+n}^-(s)} $ $ ig ] $        (3.22)

这个式子括号里的部分可以按如下方法求解：

1) 找出其反变换

2) 移位$alpha$

3) 求单边Laplace变换

4. 例子

设有效信号为Gauss-Markov过程，方差$sigma^2$，时间常数$eta$；噪声为白噪声，功率谱为常数$A$。两者不相关。我们分别使用参数滤波器、非因果维纳滤波和因果维纳滤波三种方法处理。

4.1 单参数滤波器

我们使用一阶滤波器$G(s)=frac{1}{1+Ts}$。

$G(s)=frac{1/T}{s+1/T}$

$1-G(s)=frac{s}{s+1/T}$

$S_x(s)=frac {2sigma^2eta} {-s^2+eta ^2} = {frac {sqrt {2sigma^2eta} } {s+eta} } cdot {frac {sqrt {2sigma^2eta} } {-s+eta}  } $

$S_n(s)=A=sqrt{A} cdot sqrt{A}$

上面几个式子代入(2.5)，查积分表可得

$E[e^2] = frac {sigma^2 eta T} { 1 + eta T } + frac {A} {2T} $

对上式求最小值，就得到

$T = frac { sqrt A } {sigma sqrt {2eta} - eta sqrt {A} } $

4.2 非因果滤波器

为简化计算，我们令$sigma^2=eta=A=1$，$alpha=0$。因为有效信号和噪声不相关，所以

$S_{x+n}=S_x+S_n=frac{2}{-s^2+1}+1=frac{-s^2+3}{-s^2+1}$

$S_{x+n,x}=S_x=frac{2}{-s^2+1}$

$e^{alpha s}=1$

根据(3.10)，

$G(s)=frac {frac {2} {-s^2+1}} {frac{-s^2+3}{-s^2+1}}=frac{2}{-s^2+3}$

部分分式展开为

$G(s)=frac{1/sqrt 3}{s+sqrt 3} + frac {1/sqrt 3} {-s + sqrt 3} $

$g(t)=frac{1}{sqrt 3}e^{-sqrt 3 t} u(t) + frac{1}{sqrt 3} e^{sqrt 3 t}u(-t)$

4.3 因果滤波器

$S_{x+n}=S_x+S_n = frac{2}{-s^2+1}+1=frac{-s^2+3}{-s^2+1}$

$S_{x+n, x}=S_x=frac{2}{-s^2+1}$

$e^{alpha s}=1$

首先，将$S_{x+n}$按零极点在正负平面的分布进行分解：

$S_{x+n}=S_{x+n}^+S_{x+n}^-=[frac{s+sqrt 3}{s+1}][frac{-s+sqrt 3}{-s+1}]$

下一步，求得$S_{x+n,x}/S_{x+n}^-$

$frac{S_{x+n,x}}{S_{x+n}^-} = frac{frac{2}{-s^2+1}} { frac{-s+sqrt 3 }{-s+1} } = frac{2} {(-s+sqrt 3)(s+1)} $

接着，对上面的式子进行分式分解：

$frac{S_{x+n,x}} {S_{x+n}^-}=frac{sqrt 3 - 1} {s+1} + frac{sqrt 3 - 1} {-s+sqrt 3}$

上面式子的第一个分式就是正时间部分。按(3.22)，

$G(s)=frac{1}{frac{s+sqrt 3}{s+1} } cdot frac{sqrt 3 - 1} {s+1} = frac {sqrt 3 - 1}{s + sqrt 3}$

$g(t)=(sqrt 3 -1)e^{-sqrt 3 t} u(t)$

5. 总结

单参数滤波器需要事先确定滤波器形式，例子中仅考虑一阶情况，效果在上面三种滤波器中是最差的；非因果滤波器在维纳滤波器下几乎没有限制，效果最佳；因果滤波器不能利用“未来”的值，效果在中间。

维纳滤波器要求至少了解信号和噪声的自相关函数，如果信号和噪声相关，还要了解它们的互相关，这对维纳滤波的应用范围造成了一定的限制。