第二章 维纳滤波 随机信号笔记
第二章 维纳滤波
2.1 引言
信号检测目标：从噪声中提取信号。
对于一个线性系统来说，如果其冲激响应为
h
(
n
)
h(n)
h(n) ，则当输入某一随机信号
x
(
n
)
x(n)
x(n) 时，它的输入为
y
(
n
)
=
∑
m
h
(
m
)
x
(
n
−
m
)
y(n)=\sum_{m}h(m)x(n-m)
y(n)=m∑h(m)x(n−m)
这里的输入往往是掺杂着噪声的信号，噪声是我们无法避免的，我们只能尽量减小噪声对信号的影响。
x
(
n
)
=
s
(
n
)
+
v
(
n
)
x(n)=s(n)+v(n)
x(n)=s(n)+v(n)
式中
s
(
n
)
s(n)
s(n) 表示信号，
v
(
n
)
v(n)
v(n) 表示噪声。
这时我们希望这种线性系统的输出是尽量逼急
s
(
n
)
s(n)
s(n) 的，即
y
(
n
)
=
s
^
(
n
)
y(n) = \widehat{s}(n)
y(n)=s
(n)
而我们的目标就是求出能够使均方误差（
E
[
e
2
(
n
)
]
=
E
[
(
s
(
n
)
−
s
^
(
n
)
)
2
]
,
e
(
n
)
=
s
(
n
)
−
s
^
(
n
)
E[e^2(n)]=E[(s(n)-\widehat{s}(n))^2],\quad e(n)=s(n)-\widehat{s}(n)
E[e2(n)]=E[(s(n)−s
(n))2],e(n)=s(n)−s
(n) ）最小的冲激响应函数
h
(
n
)
h(n)
h(n) 。
这种估计器的主要功能是利用当前的观测值以及一系列过去的观测值来完成对当前信号值的某种估计。
已知：
x
(
n
)
,
x
(
n
−
1
)
,
x
(
n
−
2
)
,
⋯
x(n),x(n-1),x(n-2),\cdots
x(n),x(n−1),x(n−2),⋯ 期望输出：
d
(
n
)
=
{
s
^
(
n
)
滤
波
s
^
(
n
+
N
)
预
测
s
^
(
n
−
N
)
内
插
/
平
滑
d(n)= \left\{ \begin{aligned} &\widehat{s}(n)\quad &滤波\\ &\widehat{s}(n+N)\quad &预测\\ &\widehat{s}(n-N)\quad &内插/平滑 \end{aligned} \right.
d(n)=⎩⎪⎨⎪⎧s
(n)s
(n+N)s
(n−N)滤波预测内插/平滑
2.2 维纳滤波器的时域解
2.2.1 教材内容
y ( n ) = s ^ ( n ) = ∑ m h ( m ) x ( n − m ) y(n)=\widehat{s}(n)=\sum_{m}h(m)x(n-m) y(n)=s (n)=m∑h(m)x(n−m)
如果系统是物理可实现（因果系统）的，即
h
(
n
)
=
0
,
n
<
0
h(n)=0,\qquad n<0
h(n)=0,n<0
则此时的
y
(
n
)
=
s
^
(
n
)
=
∑
m
=
0
∞
h
(
m
)
x
(
n
−
m
)
y(n)=\widehat{s}(n)=\sum_{m=0}^{\infty}h(m)x(n-m)
y(n)=s
(n)=m=0∑∞h(m)x(n−m)
维纳滤波器的设计则是要确定均方误差
E
[
e
2
(
n
)
]
=
E
[
(
s
(
n
)
−
s
^
(
n
)
)
2
]
E[e^2(n)]=E[(s(n)-\widehat{s}(n))^2]
E[e2(n)]=E[(s(n)−s
(n))2]
最小意义下的冲激响应
h
o
p
t
(
n
)
h_{opt}(n)
hopt(n) 。
为了方便得出矩阵表示，我们现将
s
^
(
n
)
=
∑
m
=
0
∞
h
(
m
)
x
(
n
−
m
)
\widehat{s}(n)=\sum_{m=0}^{\infty}h(m)x(n-m)
s
(n)=m=0∑∞h(m)x(n−m)
改写成
s
^
(
n
)
=
∑
i
=
1
∞
h
i
x
i
i
=
m
+
1
,
或
m
=
i
−
1
h
i
=
h
(
m
)
=
h
(
i
−
1
)
x
i
=
x
(
n
−
m
)
=
x
(
n
−
i
−
1
)
}
\widehat{s}(n)=\sum_{i=1}^{\infty}h_ix_i \\ \left. \begin{matrix} i=m+1,或m=i-1\\ h_i=h(m)=h(i-1)\\ x_i=x(n-m)=x(n-i-1) \end{matrix} \right\}
s
(n)=i=1∑∞hixii=m+1,或m=i−1hi=h(m)=h(i−1)xi=x(n−m)=x(n−i−1)⎭⎬⎫
这时我们将均方误差对
h
j
h_j
hj 求偏导
∂
E
[
e
2
(
n
)
]
∂
h
j
=
E
[
2
(
s
−
(
h
1
x
1
+
h
2
x
2
+
⋯
)
)
(
−
x
j
)
]
,
j
≥
1
\frac{\partial E[e^2(n)]}{\partial h_j}= E[2(s-(h_1x_1+h_2x_2+\cdots))(-x_j)], \quad j \geq 1
∂hj∂E[e2(n)]=E[2(s−(h1x1+h2x2+⋯))(−xj)],j≥1
令偏导数等于0得
E
[
(
s
−
∑
i
=
1
∞
h
i
x
i
)
x
j
]
=
0
,
j
≥
1
E
[
e
x
j
]
=
0
,
j
≥
1
E[(s-\sum_{i=1}^{\infty}h_ix_i)x_j]=0, \quad j \geq 1\\ E[ex_j]=0, \quad j \geq 1
E[(s−i=1∑∞hixi)xj]=0,j≥1E[exj]=0,j≥1
由于上式符合正交性原理，即满足两个矢量的内积为0，故满足正交性原理与满足均方误差最小的条件是一致的。由于
E
[
x
i
x
j
]
=
ϕ
x
i
,
x
j
E[x_ix_j]=\phi_{x_i,x_j}
E[xixj]=ϕxi,xj 以及
E
[
x
i
s
]
=
ϕ
x
i
,
s
E[x_is]=\phi_{x_i,s}
E[xis]=ϕxi,s 可得
ϕ
x
j
s
=
∑
i
=
1
∞
h
i
ϕ
x
i
x
j
\phi_{x_js}=\sum_{i=1}^{\infty}h_i\phi_{x_ix_j}
ϕxjs=i=1∑∞hiϕxixj
所以有
E
{
[
s
(
n
)
−
∑
m
=
0
∞
h
o
p
t
(
m
)
x
(
n
−
m
)
]
x
(
n
−
k
)
}
=
0
,
k
≥
0
E\{[s(n)-\sum_{m=0}^{\infty}h_{opt}(m)x(n-m)]x(n-k)\}=0,\quad k \geq 0
E{[s(n)−m=0∑∞hopt(m)x(n−m)]x(n−k)}=0,k≥0
即（维纳-霍夫方程）
ϕ
x
s
(
k
)
=
∑
m
=
0
∞
h
o
p
t
(
m
)
ϕ
x
x
(
k
−
m
)
,
k
≥
0
\phi_{xs}(k)=\sum_{m=0}^{\infty}h_{opt}(m)\phi_{xx}(k-m), \quad k \geq 0
ϕxs(k)=m=0∑∞hopt(m)ϕxx(k−m),k≥0
如果没有物理条件限制（因果系统
k
≥
0
k\geq 0
k≥0）那么非因果的维纳-霍夫方程可以变换到
z
z
z 域，
Φ
x
s
(
z
)
=
H
o
p
t
(
z
)
Φ
x
x
(
z
)
H
o
p
t
(
z
)
=
Φ
x
s
(
z
)
Φ
x
x
(
z
)
\Phi_{xs}(z)=H_{opt}(z)\Phi_{xx}(z)\\ H_{opt}(z)=\frac{\Phi_{xs}(z)}{\Phi_{xx}(z)}
Φxs(z)=Hopt(z)Φxx(z)Hopt(z)=Φxx(z)Φxs(z)
从而很方便地解出
h
o
p
t
(
n
)
=
Z
−
1
[
Φ
x
s
(
z
)
Φ
x
x
(
z
)
]
h_{opt}(n)=Z^{-1}[\frac{\Phi_{xs}(z)}{\Phi_{xx}(z)}]
hopt(n)=Z−1[Φxx(z)Φxs(z)]
对于有物理条件限制的情况，对于
x
(
n
)
x(n)
x(n) 而言相当于只能取当前的与过去的观测数据。如果从
{
x
(
n
)
}
\{x(n)\}
{x(n)} 求
s
^
(
n
)
\widehat{s}(n)
s
(n) 时，并不严格要求实时得到有关结果，二是允许在将来（例如可以经过一定等待或延迟）得到相应的结果，这在数字信号处理中实现起来比较方便，那么
y
(
n
)
=
s
^
(
n
)
=
∑
m
=
0
∞
h
(
m
)
x
(
n
−
m
)
y(n)=\widehat{s}(n)=\sum_{m=0}^{\infty}h(m)x(n-m)
y(n)=s
(n)=m=0∑∞h(m)x(n−m)
中的求和下限不必限于
m
=
0
m=0
m=0 ，而是可以让
m
<
0
m< 0
m<0 了，也没有了
h
(
m
<
0
)
=
0
h(m<0)=0
h(m<0)=0 的限制，因而前述的维纳-霍夫方程就也不再有
k
≥
0
k \geq 0
k≥0 的限制了。
在要求严格实时处理的场合，不允许有过多的等待和延迟，则必须考虑因果性限制，在工程实际中常在时域逼近的方法解上述方程，这时如果具有因果性约束的
h
(
n
)
h(n)
h(n) 可用长度为
N
N
N 的有限长序列来逼近，则上述方程可变成如下形式：
y
(
n
)
=
s
^
(
n
)
=
∑
m
=
0
∞
h
(
m
)
x
(
n
−
m
)
s
^
(
n
)
=
∑
m
=
0
N
−
1
h
(
m
)
ϕ
x
x
(
k
−
m
)
,
k
=
0
,
1
,
2
,
⋯
,
N
y(n)=\widehat{s}(n)=\sum_{m=0}^{\infty}h(m)x(n-m)\\ \widehat{s}(n)=\sum_{m=0}^{N-1}h(m)\phi_{xx}(k-m), \quad k=0, 1,2,\cdots,N
y(n)=s
(n)=m=0∑∞h(m)x(n−m)s
(n)=m=0∑N−1h(m)ϕxx(k−m),k=0,1,2,⋯,N
和
E
[
(
s
−
∑
i
=
1
N
h
i
x
i
)
x
j
]
=
0
,
j
=
0
,
1
,
2
,
⋯
,
N
E
[
e
x
j
]
=
0
,
j
=
0
,
1
,
2
,
⋯
,
N
E[(s-\sum_{i=1}^{N}h_ix_i)x_j]=0, \quad j=0, 1,2,\cdots,N\\ E[ex_j]=0, \quad j =0, 1,2,\cdots,N
E[(s−i=1∑Nhixi)xj]=0,j=0,1,2,⋯,NE[exj]=0,j=0,1,2,⋯,N
以及
ϕ
x
j
s
=
∑
i
=
1
∞
h
i
ϕ
x
i
x
j
,
j
=
0
,
1
,
2
,
⋯
,
N
ϕ
x
s
(
k
)
=
∑
m
=
0
∞
h
o
p
t
(
m
)
ϕ
x
x
(
k
−
m
)
,
k
=
0
,
1
,
2
,
⋯
,
N
\phi_{x_js}=\sum_{i=1}^{\infty}h_i\phi_{x_ix_j}, \quad j=0, 1,2,\cdots,N\\ \phi_{xs}(k)=\sum_{m=0}^{\infty}h_{opt}(m)\phi_{xx}(k-m), \quad k=0, 1,2,\cdots,N
ϕxjs=i=1∑∞hiϕxixj,j=0,1,2,⋯,Nϕxs(k)=m=0∑∞hopt(m)ϕxx(k−m),k=0,1,2,⋯,N
矩阵表达形式
[
ϕ
x
s
]
=
[
ϕ
x
x
]
[
h
]
[
ϕ
x
s
]
=
[
ϕ
x
1
s
ϕ
x
2
s
⋮
ϕ
x
N
s
]
[
ϕ
x
x
]
=
[
ϕ
x
1
x
1
ϕ
x
1
x
2
⋯
ϕ
x
1
x
N
ϕ
x
2
x
1
ϕ
x
2
x
2
⋯
ϕ
x
2
x
N
⋮
⋮
⋮
ϕ
x
N
x
1
ϕ
x
N
x
2
⋯
ϕ
x
N
x
N
]
[
h
]
=
[
h
1
h
2
⋮
h
N
]
[\phi_{xs}]=[\phi_{xx}][h]\\ \quad \\ [\phi_{xs}]= \left[ \begin{matrix} \phi_{x_1s} \\ \phi_{x_2s} \\ \vdots \\ \phi_{x_Ns} \\ \end{matrix} \right] \quad [\phi_{xx}]= \left[ \begin{matrix} \phi_{x_1x_1} & \phi_{x_1x_2} & \cdots & \phi_{x_1x_N} \\ \phi_{x_2x_1} & \phi_{x_2x_2} & \cdots & \phi_{x_2x_N} \\ \vdots & \vdots & & \vdots \\ \phi_{x_Nx_1} & \phi_{x_Nx_2} & \cdots & \phi_{x_Nx_N} \\ \end{matrix} \right] \quad [h]= \left[ \begin{matrix} h_{1} \\ h_{2} \\ \vdots \\ h_{N} \\ \end{matrix} \right]
[ϕxs]=[ϕxx][h][ϕxs]=⎣⎢⎢⎢⎡ϕx1sϕx2s⋮ϕxNs⎦⎥⎥⎥⎤[ϕxx]=⎣⎢⎢⎢⎡ϕx1x1ϕx2x1⋮ϕxNx1ϕx1x2ϕx2x2⋮ϕxNx2⋯⋯⋯ϕx1xNϕx2xN⋮ϕxNxN⎦⎥⎥⎥⎤[h]=⎣⎢⎢⎢⎡h1h2⋮hN⎦⎥⎥⎥⎤
故可以解得
[
h
]
=
[
h
o
p
t
]
=
[
ϕ
x
x
]
−
1
[
ϕ
x
s
]
[h]=[h_{opt}]=[\phi_{xx}]^{-1}[\phi_{xs}]
[h]=[hopt]=[ϕxx]−1[ϕxs]
2.2.2 课堂笔记
从另一个角度上来看 E [ e ( n ) x ( n − k ) ] = 0 , k ≥ 0 E[e(n)x(n-k)]=0,\quad k \geq 0 E[e(n)x(n−k)]=0,k≥0 满足正交性原理与满足最小均方误差是等价的，我们可以假设 h ( n ) h(n) h(n) 的长度为 N 。令：
输入矢量： X ( n ) = [ x ( n ) , x ( n − 1 ) , ⋯ , x ( n − k + 1 ) ] T X(n)=[x(n),x(n-1),\cdots,x(n-k+1)]^T X(n)=[x(n),x(n−1),⋯,x(n−k+1)]T
权向量： W = [ w 0 , w 1 , ⋯ , w N − 1 ] T = [ h ( 0 ) , h ( 1 ) , ⋯ , h ( N − 1 ) ] T W=[w_0,w_1,\cdots,w_{N-1}]^T=[h(0),h(1),\cdots,h(N-1)]^T W=[w0,w1,⋯,wN−1]T=[h(0),h(1),⋯,h(N−1)]T
滤波器输出： y ( n ) = W T X ( n ) = X T ( n ) W y(n)=W^TX(n)=X^T(n)W y(n)=WTX(n)=XT(n)W
均方误差：
E
[
e
2
(
n
)
]
=
E
[
(
d
(
n
)
−
W
T
X
(
n
)
)
(
d
(
n
)
−
X
(
n
)
T
W
)
]
=
E
[
d
2
(
n
)
]
−
E
[
d
(
n
)
X
(
n
)
T
W
]
−
E
[
d
(
n
)
W
T
X
(
n
)
]
−
W
T
E
[
X
(
n
)
X
(
n
)
T
]
W
\begin{aligned} E[e^2(n)]&=E[(d(n)-W^TX(n))(d(n)-X(n)^TW)]\\ &=E[d^2(n)]-E[d(n)X(n)^TW]-E[d(n)W^TX(n)]-W^TE[X(n)X(n)^T]W \end{aligned}
E[e2(n)]=E[(d(n)−WTX(n))(d(n)−X(n)TW)]=E[d2(n)]−E[d(n)X(n)TW]−E[d(n)WTX(n)]−WTE[X(n)X(n)T]W
令互相关矢量
P
=
E
[
X
(
n
)
d
(
n
)
]
=
{
ϕ
x
d
(
0
)
ϕ
x
d
(
1
)
⋮
ϕ
x
d
(
N
−
1
)
P=E[X(n)d(n)]= \left\{ \begin{aligned} &\phi_{xd}(0) \\ &\phi_{xd}(1) \\ &\vdots \\ &\phi_{xd}(N-1) \\ \end{aligned} \right.
P=E[X(n)d(n)]=⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧ϕxd(0)ϕxd(1)⋮ϕxd(N−1)
d
(
n
)
d(n)
d(n) 的平均功率：
σ
d
2
=
E
[
d
2
(
n
)
]
\sigma_{d}^2=E[d^2(n)]
σd2=E[d2(n)]
代价函数
J
(
w
)
=
E
[
e
2
(
n
)
]
=
σ
d
2
−
P
T
W
−
W
T
P
+
W
T
R
x
W
=
σ
d
2
−
2
W
T
P
+
W
T
R
x
W
\begin{aligned} J(w)&=E[e^2(n)] \\ &=\sigma_d^2-P^TW-W^TP+W^TR_xW \\ &=\sigma_d^2-2W^TP+W^TR_xW \end{aligned}
J(w)=E[e2(n)]=σd2−PTW−WTP+WTRxW=σd2−2WTP+WTRxW
J
(
w
)
J(w)
J(w) 是
w
w
w 的二次函数，有唯一的极小值点，下面用两种方法求解使得
J
(
w
)
J(w)
J(w) 取最小值时
w
w
w 的值：
- 配方法：
令 U = R − 1 P ⟶ ( W − U ) T R ( W − U ) = W T R W − 2 W T P + P T R − 1 P U=R^{-1}P\quad\longrightarrow (W-U)^TR(W-U)=W^TRW-2W^TP+P^TR^{-1}P U=R−1P⟶(W−U)TR(W−U)=WTRW−2WTP+PTR−1P
故 J ( w ) = σ d 2 − P T R − 1 P + ( W − R − 1 P ) T R ( W − R − 1 P ) J(w)=\sigma_d^2-P^TR^{-1}P+(W-R^{-1}P)^TR(W-R^{-1}P) J(w)=σd2−PTR−1P+(W−R−1P)TR(W−R−1P)
R
R
R 是正定的，
R
−
1
R^{-1}
R−1 也是正定的，所以
h
o
p
t
=
W
o
p
t
=
R
−
1
P
h_{opt}=W_{opt}=R^{-1}P
hopt=Wopt=R−1P ，最小值
J
m
i
n
=
J
(
W
o
p
t
)
=
σ
d
2
−
P
T
R
−
1
P
=
σ
d
2
−
P
T
W
o
p
t
\begin{aligned} J_{min}&=J(W_{opt}) \\ &=\sigma_d^2-P^TR^{-1}P \\ &=\sigma_d^2-P^TW_{opt} \end{aligned}
Jmin=J(Wopt)=σd2−PTR−1P=σd2−PTWopt
- 求导法
设 X = [ x 1 , x 2 , ⋯ , x n ] f ( X ) = f ( x 1 , x 2 , ⋯ , x n ) X=[x_1,x_2,\cdots,x_n] \quad f(X)=f(x_1,x_2,\cdots,x_n) X=[x1,x2,⋯,xn]f(X)=f(x1,x2,⋯,xn)
则
f
f
f 的微分形式如下
∇
f
(
X
)
=
[
∂
f
∂
x
1
,
∂
f
∂
x
2
,
⋯
,
∂
f
∂
x
n
]
\nabla f(X)=[\frac{\partial f}{\partial x_1},\frac{\partial f}{\partial x_2,} \cdots,\frac{\partial f}{\partial x_n}]
∇f(X)=[∂x1∂f,∂x2,∂f⋯,∂xn∂f]
∇ 2 f ( X ) = [ ∂ 2 f ∂ x 1 ∂ x 1 ∂ 2 f ∂ x 1 ∂ x 2 ⋯ ∂ 2 f ∂ x 1 ∂ x n ⋮ ⋮ ⋮ ∂ 2 f ∂ x n ∂ x 1 ∂ 2 f ∂ x n ∂ x 2 ⋯ ∂ 2 f ∂ x n ∂ x n ] (12) \nabla^2f(X)= \left[ \begin{array}{cc} \frac{\partial^2 f}{\partial x_1\partial x_1} & \frac{\partial^2 f}{\partial x_1\partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1\partial x_n} \\ \vdots & \vdots & & \vdots \\ \frac{\partial^2 f}{\partial x_n\partial x_1} & \frac{\partial^2 f}{\partial x_n\partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n\partial x_n} \\ \end{array}\tag{12} \right] ∇2f(X)=⎣⎢⎢⎡∂x1∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f⋮∂xn∂x2∂2f⋯⋯∂x1∂xn∂2f⋮∂xn∂xn∂2f⎦⎥⎥⎤(12)
若
f
(
X
)
=
1
2
X
T
A
X
+
b
T
X
+
C
f(X)=\frac{1}{2}X^TAX+b^TX+C
f(X)=21XTAX+bTX+C 其中
A
A
A 是对称阵，
b
b
b 是列向量，
C
C
C 是常数。则
∇
f
(
X
)
=
A
X
+
b
∇
2
f
(
X
)
=
A
\nabla f(X)=AX+b\\ \nabla^2f(X)=A
∇f(X)=AX+b∇2f(X)=A
所以
J
(
W
)
=
W
T
R
x
W
−
2
P
T
W
+
σ
d
2
J(W)=W^TR_xW-2P^TW+\sigma_d^2
J(W)=WTRxW−2PTW+σd2 对
W
W
W 求导，得
∇
J
(
W
)
=
2
R
W
−
2
P
=
0
\nabla J(W)=2RW-2P=0
∇J(W)=2RW−2P=0
解得
h
o
p
t
=
W
o
p
t
=
R
−
1
P
h_{opt}=W_{opt}=R^{-1}P
hopt=Wopt=R−1P
- 坐标系的平移旋转
考虑到在计算均方误差时二次项 W T R x W W^TR_xW WTRxW 的结果中含有交叉乘积项，不便于计算分析，这里我们对坐标系进行平移旋转，使得表达式尽量的简单便于计算。
首先我们根据上述的
W
W
W 的最优解，对
W
W
W 进行平移：
W
~
=
W
−
W
o
p
t
\widetilde{W}=W-W_{opt}
W
=W−Wopt
这样将最优解的点平移到了 零向量处，考虑讲矩阵
R
R
R 对角化
R
=
Q
Λ
Q
−
1
=
Q
Λ
Q
T
(
Q
是
正
交
阵
)
R=Q\Lambda Q^{-1}=Q\Lambda Q^T(Q是正交阵)
R=QΛQ−1=QΛQT(Q是正交阵)
原式变形为
J
=
J
m
i
n
+
W
~
T
Q
Λ
Q
T
W
~
=
J
m
i
n
+
(
Q
T
W
~
)
T
Λ
(
Q
T
W
~
)
\begin{aligned} J&=J_{min}+\widetilde{W}^TQ\Lambda Q^T\widetilde{W} \\ &=J_{min}+(Q^T\widetilde{W})^T\Lambda (Q^T\widetilde{W}) \\ \end{aligned}
J=Jmin+W
TQΛQTW
=Jmin+(QTW
)TΛ(QTW
)
令
V
=
Q
T
W
~
V=Q^T\widetilde{W}
V=QTW
得：
J
=
J
m
i
n
+
V
T
Λ
V
=
J
m
i
n
+
∑
i
=
0
N
−
1
λ
i
v
i
2
\begin{aligned} J&=J_{min}+V^T \Lambda V \\ &=J_{min}+\sum_{i=0}^{N-1}\lambda_{i}v_{i}^2 \\ \end{aligned}
J=Jmin+VTΛV=Jmin+i=0∑N−1λivi2
从数学上考虑当
Q
Q
Q 为单位阵时，即
R
R
R 本身就是对角阵 ，求解最简单。这时我们考虑物理情景，自相关函数只有对角线上不为零，这时显然该信号为冲激函数。当
R
R
R 本身不为对角阵时 考虑讲信号通过白噪声的逆滤波器（信号经过该滤波器后输出高斯白噪声），后再进行上述过程，这便是下一节处理维纳滤波的 z 域解的主要处理方法。
2.3 维纳滤波器的 z 域解
根据上一节所提到的方法，任何一个具有功率谱密度的随机信号都可以看作是由一个白噪声
w
(
n
)
w(n)
w(n) 激励某个物理网络所得到的，若信号用
s
(
n
)
s(n)
s(n) 表示，则有
S
(
z
)
=
A
(
z
)
W
(
z
)
S(z)=A(z)W(z)
S(z)=A(z)W(z) ，白噪声的自相关函数以及功率谱密度分别为
ϕ
w
w
(
n
)
=
σ
w
2
δ
(
n
)
Φ
w
w
(
z
)
=
σ
w
2
\phi_{ww}(n)=\sigma_w^2\delta(n) \\ \Phi_{ww}(z)=\sigma_w^2
ϕww(n)=σw2δ(n)Φww(z)=σw2
当信号模型的冲激响应为实序列时，
s
(
n
)
s(n)
s(n) 的功率谱密度为
Φ
s
s
(
z
)
=
σ
w
2
A
(
z
)
A
(
z
−
1
)
\Phi_{ss}(z)=\sigma_w^2A(z)A(z^{-1})
Φss(z)=σw2A(z)A(z−1)
对于实际信号
x
(
n
)
=
s
(
n
)
+
v
(
n
)
x(n)=s(n)+v(n)
x(n)=s(n)+v(n) 也可以使用这种模型来表示。
Φ
x
x
(
z
)
=
σ
w
2
B
(
z
)
B
(
z
−
1
)
X
(
z
)
=
B
(
z
)
W
(
z
)
\Phi_{xx}(z)=\sigma_w^2B(z)B(z^{-1})\\ X(z)=B(z)W(z)
Φxx(z)=σw2B(z)B(z−1)X(z)=B(z)W(z)
对于一个系统
B
(
z
)
B(z)
B(z) 稳定的条件为极点在单位圆内，零点无要求。但是我们为了实现逆滤波，需要构建物理系统
1
B
(
z
)
\frac{1}{B(z)}
B(z)1 ，因为 系统
B
(
z
)
B(z)
B(z) 的零点就是系统
1
B
(
z
)
\frac{1}{B(z)}
B(z)1 的极点， 所以要求系统
B
(
z
)
B(z)
B(z) 的零点也要位于单位圆内，这样的的系统，我们叫做最小相位系统。
W
(
z
)
=
X
(
z
)
B
(
z
)
W(z)=\frac{X(z)}{B(z)}
W(z)=B(z)X(z)
为了方便得到
H
o
p
t
(
z
)
H_{opt}(z)
Hopt(z) 我们通过将信号白化来实现
x
(
n
)
=
s
(
n
)
+
v
(
n
)
⟶
H
(
z
)
⟶
y
(
n
)
=
s
^
(
n
)
x(n)=s(n)+v(n)\longrightarrow H(z) \longrightarrow y(n)=\widehat{s}(n)
x(n)=s(n)+v(n)⟶H(z)⟶y(n)=s
(n)
x ( n ) ⟶ 1 B ( z ) ⟶ w ( n ) G ( z ) ⟶ y ( n ) = s ^ ( n ) x(n)\longrightarrow \frac{1}{B(z)} \stackrel{w(n)}{\longrightarrow} G(z)\longrightarrow y(n)=\widehat{s}(n) x(n)⟶B(z)1⟶w(n)G(z)⟶y(n)=s (n)
这里我们就能相对简单地得到
H
(
z
)
=
G
(
z
)
B
(
z
)
H(z)=\frac{G(z)}{B(z)}
H(z)=B(z)G(z)
2.3.1 非因果维纳滤波器
根据前面的讨论，可得
S
^
(
z
)
=
G
(
z
)
W
(
z
)
⟶
s
^
(
n
)
=
∑
k
=
−
∞
∞
g
(
k
)
w
(
n
−
k
)
\widehat{S}(z)=G(z)W(z) \longrightarrow \widehat{s}(n)=\sum_{k=-\infty}^{\infty}g(k)w(n-k)
S
(z)=G(z)W(z)⟶s
(n)=k=−∞∑∞g(k)w(n−k)
所以均方误差
E
[
e
2
(
n
)
]
=
E
{
[
s
(
n
)
−
∑
k
=
−
∞
∞
g
(
k
)
w
(
n
−
k
)
]
2
}
=
E
{
s
2
(
n
)
−
2
∑
k
=
−
∞
∞
g
(
k
)
w
(
n
−
k
)
s
(
n
)
+
∑
k
=
−
∞
∞
∑
r
=
−
∞
∞
g
(
k
)
g
(
r
)
w
(
n
−
k
)
w
(
n
−
r
)
}
=
E
[
s
2
(
n
)
]
−
2
∑
k
=
−
∞
∞
g
(
k
)
E
[
w
(
n
−
k
)
s
(
n
)
]
+
∑
k
=
−
∞
∞
∑
r
=
−
∞
∞
g
(
k
)
g
(
r
)
E
[
w
(
n
−
k
)
w
(
n
−
r
)
]
=
ϕ
s
s
(
0
)
−
2
∑
k
=
−
∞
∞
g
(
k
)
ϕ
w
s
(
k
)
+
∑
k
=
−
∞
∞
∑
r
=
−
∞
∞
g
(
k
)
g
(
r
)
ϕ
w
w
(
k
−
r
)
\begin{aligned} E[e^2(n)]&=E\{[s(n)-\sum_{k=-\infty}^{\infty}g(k)w(n-k)]^2\} \\ &=E\{ s^2(n)-2\sum_{k=-\infty}^{\infty}g(k)w(n-k)s(n) \\ &\quad +\sum_{k=-\infty}^{\infty} \sum_{r=-\infty}^{\infty}g(k)g(r)w(n-k)w(n-r)\}\\ &=E[s^2(n)]-2\sum_{k=-\infty}^{\infty}g(k)E[w(n-k)s(n)]\\ &\quad +\sum_{k=-\infty}^{\infty} \sum_{r=-\infty}^{\infty}g(k)g(r)E[w(n-k)w(n-r)]\\ &=\phi_{ss}(0)-2\sum_{k=-\infty}^{\infty}g(k)\phi_{ws}(k) +\sum_{k=-\infty}^{\infty}\sum_{r=-\infty}^{\infty} g(k)g(r)\phi_{ww}(k-r) \end{aligned}
E[e2(n)]=E{[s(n)−k=−∞∑∞g(k)w(n−k)]2}=E{s2(n)−2k=−∞∑∞g(k)w(n−k)s(n)+k=−∞∑∞r=−∞∑∞g(k)g(r)w(n−k)w(n−r)}=E[s2(n)]−2k=−∞∑∞g(k)E[w(n−k)s(n)]+k=−∞∑∞r=−∞∑∞g(k)g(r)E[w(n−k)w(n−r)]=ϕss(0)−2k=−∞∑∞g(k)ϕws(k)+k=−∞∑∞r=−∞∑∞g(k)g(r)ϕww(k−r)
由于
ϕ
w
w
(
n
)
=
σ
w
2
δ
(
n
)
\phi_{ww}(n)=\sigma_{w}^2\delta(n)
ϕww(n)=σw2δ(n) 所以上式可表示成
E
[
e
2
(
n
)
]
=
ϕ
s
s
(
0
)
−
2
∑
k
=
−
∞
∞
g
(
k
)
ϕ
w
s
(
k
)
+
∑
k
=
−
∞
∞
g
2
(
k
)
σ
w
2
=
ϕ
s
s
(
0
)
+
∑
k
=
−
∞
∞
[
σ
w
g
(
k
)
−
ϕ
w
s
(
k
)
σ
w
]
2
−
∑
k
=
−
∞
∞
ϕ
w
s
2
(
k
)
σ
w
2
\begin{aligned} E[e^2(n)]&=\phi_{ss}(0)-2\sum_{k=-\infty}^{\infty}g(k)\phi_{ws}(k) +\sum_{k=-\infty}^{\infty}g^2(k)\sigma_w^2\\ &=\phi_{ss}(0)+\sum_{k=-\infty}^{\infty} [\sigma_wg(k)-\frac{\phi_{ws}(k)}{\sigma_w}]^2 -\sum_{k=-\infty}^{\infty}\frac{\phi_{ws}^2(k)}{\sigma_w^2} \end{aligned}
E[e2(n)]=ϕss(0)−2k=−∞∑∞g(k)ϕws(k)+k=−∞∑∞g2(k)σw2=ϕss(0)+k=−∞∑∞[σwg(k)−σwϕws(k)]2−k=−∞∑∞σw2ϕws2(k)
对
g
(
k
)
g(k)
g(k) 求导，并令导数等于0
∂
E
[
e
2
(
n
)
]
∂
g
(
k
)
=
0
\frac{\partial E[e^2(n)]}{\partial g(k)}=0
∂g(k)∂E[e2(n)]=0
解得
g
o
p
t
(
k
)
=
ϕ
w
s
(
k
)
σ
w
2
,
−
∞
<
k
<
∞
G
o
p
t
(
z
)
=
g
o
p
t
(
k
)
=
1
σ
w
2
Φ
w
s
(
z
)
g_{opt}(k)=\frac{\phi_{ws}(k)}{\sigma_w^2},\qquad -\infty<k<\infty\\ G_{opt}(z)=g_{opt}(k)=\frac{1}{\sigma_w^2}\Phi_{ws}(z)
gopt(k)=σw2ϕws(k),−∞<k<∞Gopt(z)=gopt(k)=σw21Φws(z)
由于上面提到
H
(
z
)
=
G
(
z
)
B
(
z
)
H(z)=\frac{G(z)}{B(z)}
H(z)=B(z)G(z) ，则
H
o
p
t
(
z
)
=
G
o
p
t
(
z
)
B
(
z
)
=
1
σ
w
2
Φ
w
s
(
z
)
B
(
z
)
H_{opt}(z)=\frac{G_{opt}(z)}{B(z)}=\frac{1}{\sigma_w^2}\frac{\Phi_{ws}(z)}{B(z)}
Hopt(z)=B(z)Gopt(z)=σw21B(z)Φws(z)
如果
x
(
n
)
x(n)
x(n) 是由白噪声激励一个系统函数为
B
(
z
)
B(z)
B(z) 的线性系统得到的话，
x
(
n
)
x(n)
x(n) 与
s
(
n
)
s(n)
s(n) 的互相关函数
ϕ
s
x
(
m
)
=
E
[
s
(
n
)
x
(
n
+
m
)
]
=
E
[
s
(
n
)
∑
k
=
−
∞
∞
b
(
k
)
w
(
n
+
m
−
k
)
]
=
∑
k
=
−
∞
∞
b
(
k
)
E
[
s
(
n
)
w
(
n
+
m
−
k
)
]
=
∑
k
=
−
∞
∞
b
(
k
)
ϕ
s
w
(
m
−
k
)
=
b
(
m
)
∗
ϕ
s
w
(
m
)
\begin{aligned} \phi_{sx(m)}&=E[s(n)x(n+m)]=E[s(n)\sum_{k=-\infty}^\infty b(k)w(n+m-k)] \\ &=\sum_{k=-\infty}^\infty b(k)E[s(n)w(n+m-k)]\\ &=\sum_{k=-\infty}^\infty b(k)\phi_{sw}(m-k)=b(m)*\phi_{sw}(m) \end{aligned}
ϕsx(m)=E[s(n)x(n+m)]=E[s(n)k=−∞∑∞b(k)w(n+m−k)]=k=−∞∑∞b(k)E[s(n)w(n+m−k)]=k=−∞∑∞b(k)ϕsw(m−k)=b(m)∗ϕsw(m)
由于
ϕ
s
x
(
m
)
=
ϕ
x
s
(
−
m
)
\phi_{sx}(m)=\phi_{xs}(-m)
ϕsx(m)=ϕxs(−m) 故有
ϕ
x
s
(
−
m
)
=
b
(
m
)
∗
ϕ
w
s
(
−
m
)
\phi_{xs}(-m)=b(m)*\phi_{ws}(-m)
ϕxs(−m)=b(m)∗ϕws(−m) 或
ϕ
x
s
(
m
)
=
b
(
−
m
)
∗
ϕ
w
s
(
m
)
\phi_{xs}(m)=b(-m)*\phi_{ws}(m)
ϕxs(m)=b(−m)∗ϕws(m) 故
Φ
x
s
(
z
)
=
B
(
z
−
1
)
Φ
w
s
(
z
)
Φ
w
s
(
z
)
=
Φ
x
s
(
z
)
B
(
z
−
1
)
\Phi_{xs}(z)=B(z^{-1})\Phi_{ws}(z)\\ \Phi_{ws}(z)=\frac{\Phi_{xs}(z)}{B(z^{-1})}
Φxs(z)=B(z−1)Φws(z)Φws(z)=B(z−1)Φxs(z)
代入
H
o
p
t
(
z
)
=
G
o
p
t
(
z
)
B
(
z
)
=
1
σ
w
2
Φ
w
s
(
z
)
B
(
z
)
H_{opt}(z)=\frac{G_{opt}(z)}{B(z)}=\frac{1}{\sigma_w^2}\frac{\Phi_{ws}(z)}{B(z)}
Hopt(z)=B(z)Gopt(z)=σw21B(z)Φws(z) 得
H
o
p
t
(
z
)
=
1
σ
w
2
Φ
x
s
(
z
)
B
(
z
−
1
)
B
(
z
)
=
Φ
x
s
(
z
)
Φ
x
x
(
z
)
H_{opt}(z)=\frac{1}{\sigma_w^2}\frac{\Phi_{xs}(z)}{B(z^{-1})B(z)}=\frac{\Phi_{xs}(z)}{\Phi_{xx}(z)}
Hopt(z)=σw21B(z−1)B(z)Φxs(z)=Φxx(z)Φxs(z)
这与时域没有因果性约束所得的表达式完全一样。若噪声与信号无关的话则有
Φ
x
x
(
z
)
=
Φ
s
s
(
z
)
+
Φ
v
v
(
z
)
H
o
p
t
(
z
)
=
Φ
x
s
(
z
)
Φ
s
s
(
z
)
+
Φ
v
v
(
z
)
\Phi_{xx}(z)=\Phi_{ss}(z)+\Phi_{vv}(z)\\ H_{opt}(z)=\frac{\Phi_{xs}(z)}{\Phi_{ss}(z)+\Phi_{vv}(z)}
Φxx(z)=Φss(z)+Φvv(z)Hopt(z)=Φss(z)+Φvv(z)Φxs(z)
或
H
o
p
t
(
e
j
ω
)
=
P
x
s
(
ω
)
P
s
s
(
ω
)
+
P
v
v
(
ω
)
H_{opt}(e^{j\omega})=\frac{P_{xs}(\omega)}{P_{ss}(\omega)+P_{vv}(\omega)}
Hopt(ejω)=Pss(ω)+Pvv(ω)Pxs(ω)
2.3.2 因果滤波器
因为有物理可实现约束，这是有
g
(
k
)
=
0
,
k
<
0
g(k)=0,\quad k < 0
g(k)=0,k<0 则
s
^
(
n
)
=
∑
k
=
0
∞
g
(
k
)
w
(
n
−
k
)
E
[
e
2
(
n
)
]
=
ϕ
s
s
(
0
)
+
∑
k
=
0
∞
[
σ
w
g
(
k
)
−
ϕ
w
s
(
k
)
σ
w
]
2
−
1
σ
w
2
∑
k
=
0
∞
ϕ
w
s
2
(
k
)
\widehat{s}(n)=\sum_{k=0}^\infty g(k)w(n-k)\\ E[e^2(n)]=\phi_{ss}(0)+\sum_{k=0}^\infty[\sigma_wg(k)-\frac{\phi_{ws}(k)}{\sigma_w}]^2 -\frac{1}{\sigma_w^2}\sum_{k=0}^\infty \phi_{ws}^2(k)
s
(n)=k=0∑∞g(k)w(n−k)E[e2(n)]=ϕss(0)+k=0∑∞[σwg(k)−σwϕws(k)]2−σw21k=0∑∞ϕws2(k)
而且有
g
o
p
t
(
n
)
=
{
1
σ
w
2
ϕ
w
s
(
n
)
,
n
≥
0
0
,
n
<
0
g_{opt}(n)= \left\{ \begin{aligned} \frac{1}{\sigma_w^2}\phi_{ws(n)},\qquad&n \geq 0\\ 0,\qquad &n < 0 \end{aligned} \right.
gopt(n)=⎩⎪⎨⎪⎧σw21ϕws(n),0,n≥0n<0
这里
g
o
p
t
(
n
)
g_{opt}(n)
gopt(n) 显然是一个因果序列，所以系统的极点均在单位圆内。
G
o
p
t
(
z
)
=
1
σ
w
2
[
Φ
w
s
(
z
)
]
+
G_{opt}(z)=\frac{1}{\sigma_w^2}[\Phi_{ws}(z)]_+
Gopt(z)=σw21[Φws(z)]+
H o p t ( z ) = 1 B ( z ) G o p t = 1 σ w 2 B ( z ) [ Φ w s ( z ) ] + = 1 σ w 2 B ( z ) [ Φ x s ( z ) B ( z − 1 ) ] + \begin{aligned} H_{opt}(z)&=\frac{1}{B(z)}G_{opt}=\frac{1}{\sigma_w^2B(z)}[\Phi_{ws}(z)]_+\\ &=\frac{1}{\sigma_w^2B(z)}[\frac{\Phi_{xs}(z)}{B(z^{-1})}]_+ \end{aligned} Hopt(z)=B(z)1Gopt=σw2B(z)1[Φws(z)]+=σw2B(z)1[B(z−1)Φxs(z)]+
最小均方差
E
[
e
2
(
n
)
]
m
i
n
=
ϕ
s
s
(
0
)
−
1
σ
w
2
∑
k
=
0
∞
ϕ
w
s
2
(
k
)
=
ϕ
s
s
(
0
)
−
1
σ
w
2
∑
k
=
0
∞
[
ϕ
w
s
(
k
)
u
(
k
)
]
ϕ
w
s
(
k
)
\begin{aligned} E[e^2(n)]_{min}&=\phi_{ss}(0)-\frac{1}{\sigma_w^2}\sum_{k=0}^\infty\phi_{ws}^2(k)\\ &=\phi_{ss}(0)-\frac{1}{\sigma_w^2}\sum_{k=0}^\infty [\phi_{ws}(k)u(k)]\phi_{ws}(k) \end{aligned}
E[e2(n)]min=ϕss(0)−σw21k=0∑∞ϕws2(k)=ϕss(0)−σw21k=0∑∞[ϕws(k)u(k)]ϕws(k)
在 z 域中可表示为：
E
[
e
2
(
n
)
]
m
i
n
=
1
2
π
j
∮
C
{
Φ
s
s
(
z
)
−
1
σ
w
2
[
Φ
w
s
(
z
)
]
+
Φ
w
s
(
z
−
1
)
}
z
−
1
d
z
=
1
2
π
j
∮
C
{
Φ
s
s
(
z
)
−
1
σ
w
2
[
Φ
x
s
(
z
)
B
(
z
−
1
)
]
+
Φ
x
s
(
z
−
1
)
B
(
z
)
}
z
−
1
d
z
=
1
2
π
j
∮
C
{
Φ
s
s
(
z
)
−
H
o
p
t
(
z
)
Φ
x
s
(
z
−
1
)
}
z
−
1
d
z
\begin{aligned} E[e^2(n)]_{min}&=\frac{1}{2\pi j}\oint_C\Big\{ \Phi_{ss}(z)-\frac{1}{\sigma_w^2}[\Phi_{ws}(z)]_+\Phi_{ws}(z^{-1}) \Big\}z^{-1}dz \\ &=\frac{1}{2\pi j}\oint_C\Big\{ \Phi_{ss}(z)-\frac{1}{\sigma_w^2} \Big[\frac{\Phi_{xs}(z)}{B(z^{-1})}\Big]_+ \frac{\Phi_{xs}(z^{-1})}{B(z)} \Big\}z^{-1}dz\\ &=\frac{1}{2\pi j}\oint_C\Big\{ \Phi_{ss}(z)-H_{opt}(z)\Phi_{xs}(z^{-1}) \Big\}z^{-1}dz \end{aligned}
E[e2(n)]min=2πj1∮C{Φss(z)−σw21[Φws(z)]+Φws(z−1)}z−1dz=2πj1∮C{Φss(z)−σw21[B(z−1)Φxs(z)]+B(z)Φxs(z−1)}z−1dz=2πj1∮C{Φss(z)−Hopt(z)Φxs(z−1)}z−1dz
例题：
注意：做环路积分时要注意极点，只积单位圆内的。
2.4 维纳预测器
2.4.1 预测的可能性
随机信号的特点，任何时刻 x ( n ) x(n) x(n) 和 s ( n ) s(n) s(n) 均具有偶然性，即使该序列的过去值都知道，仍然很难确定当前时刻的真实取值。但是利用 x ( n ) x(n) x(n) 和 s ( n ) s(n) s(n) 的统计特性，借助估计方法，依然可以正确预测当前以及今后某一时刻最具有可能出现的取值。
对于任何一种信号 x ( n ) x(n) x(n) ，均可根据其自相关函数，了解任何两点的相关程度。数据间的关联越紧密，预测越可靠；完全不关联，当然就无法无误地准确预测；而白噪声由于前后毫无关联，因此也就无法预测。
随机信号预测的主要特点
1）以信号的统计特性作为预测的主要依据；
2）不可能预测误差为零的绝对精准的预测；
3）实际获得的随机信号通常是带有噪声干扰的，它使预测常与滤波联系在一起，成为带滤波的预测。不考虑噪声干扰或不带滤波的预测就是纯预测。
我们对序列的预测是严重依靠系统的惯性，如果当预测范围超过系统惯性作用范围（自相关函数为0）那么预测便不再可靠。
2.4.2 预测计算公式
对于预测功能我们期望输入输出有如下关系（N > 0）：
y
(
n
)
=
s
^
(
n
+
N
)
=
∑
m
h
(
m
)
x
(
n
−
m
)
=
∑
i
h
i
x
i
y(n)=\widehat{s}(n+N)=\sum_{m}h(m)x(n-m)=\sum_{i}h_ix_i
y(n)=s
(n+N)=m∑h(m)x(n−m)=i∑hixi
与维纳滤波器的设计相似，维纳预测器的设计实现也是所作预测的均方误差
E
[
e
2
(
n
+
N
)
]
=
E
[
(
s
(
n
+
N
)
−
s
^
(
n
+
N
)
)
2
]
E[e^2(n+N)]=E[(s(n+N)-\widehat{s}(n+N))^2]
E[e2(n+N)]=E[(s(n+N)−s
(n+N))2]
最小条件下的
h
(
n
)
h(n)
h(n) 或
H
(
z
)
H(z)
H(z) 的问题。
∂
E
[
e
2
(
n
+
N
)
]
∂
h
j
=
0
,
j
≥
1
E
[
(
s
(
n
+
N
)
−
∑
i
=
1
h
i
x
i
)
x
j
]
=
0
,
j
≥
1
\frac{\partial E[e^2(n+N)]}{\partial h_j}=0, \quad j \geq 1\\ E[(s(n+N)-\sum_{i=1}h_ix_i)x_j]=0, \quad j \geq 1\\
∂hj∂E[e2(n+N)]=0,j≥1E[(s(n+N)−i=1∑hixi)xj]=0,j≥1
即
ϕ
x
i
y
d
=
∑
i
=
1
∞
h
i
ϕ
x
i
x
j
,
j
≥
0
ϕ
x
y
d
=
∑
m
=
0
∞
h
o
p
t
(
m
)
ϕ
x
x
(
k
−
m
)
,
k
≥
0
\phi_{x_iy_d}=\sum_{i=1}^\infty h_i\phi_{x_ix_j},\quad j\geq 0\\ \phi_{xy_d}=\sum_{m=0}^\infty h_{opt}(m)\phi_{xx}(k-m),\quad k \geq 0
ϕxiyd=i=1∑∞hiϕxixj,j≥0ϕxyd=m=0∑∞hopt(m)ϕxx(k−m),k≥0
如果我们在维纳滤波器中也以
y
d
y_d
yd 表示希望得到的输出，即
y
d
(
n
)
=
s
(
n
)
y_d(n)=s(n)
yd(n)=s(n) 则有
ϕ
x
y
d
(
k
)
=
ϕ
x
s
(
k
)
=
E
[
x
(
n
)
s
(
n
+
k
)
]
\phi_{xy_d}(k)=\phi_{xs}(k)=E[x(n)s(n+k)]
ϕxyd(k)=ϕxs(k)=E[x(n)s(n+k)]
二维纳预测器的
y
d
(
n
)
=
s
(
n
+
N
)
y_d(n)=s(n+N)
yd(n)=s(n+N) ，所以有
ϕ
x
y
d
(
k
)
=
ϕ
x
s
(
k
+
N
)
=
E
[
x
(
n
)
s
(
n
+
N
+
k
)
]
\phi_{xy_d}(k)=\phi_{xs}(k+N)=E[x(n)s(n+N+k)]
ϕxyd(k)=ϕxs(k+N)=E[x(n)s(n+N+k)]
z 变换得
Φ
x
y
d
(
z
)
=
z
N
Φ
x
s
(
z
)
Φ
x
y
d
(
z
−
1
)
=
z
−
N
Φ
x
s
(
z
−
1
)
\Phi_{xy_d}(z)=z^{N}\Phi_{xs}(z) \\ \Phi_{xy_d}(z^{-1})=z^{-N}\Phi_{xs}(z^{-1})
Φxyd(z)=zNΦxs(z)Φxyd(z−1)=z−NΦxs(z−1)
无物理约束
H
o
p
t
(
z
)
=
Φ
x
y
d
(
z
)
Φ
x
x
(
z
)
=
z
N
Φ
x
s
(
z
)
Φ
x
x
(
z
)
E
{
[
s
(
n
+
N
)
−
s
^
(
n
+
N
)
]
2
}
m
i
n
=
1
2
π
j
∮
C
{
Φ
s
s
(
z
)
−
H
o
p
t
(
z
)
z
−
N
Φ
x
s
(
z
−
1
)
}
z
−
1
d
z
H_{opt}(z)=\frac{\Phi_{xy_d}(z)}{\Phi_{xx}(z)}=\frac{z^N\Phi_{xs}(z)}{\Phi_{xx}(z)}\\ E\{[s(n+N)-\widehat{s}(n+N)]^2\}_{min}=\frac{1}{2\pi j}\oint_C \Big\{ \Phi_{ss}(z)-H_{opt}(z)z^{-N}\Phi_{xs}(z^{-1}) \Big\}z^{-1}dz
Hopt(z)=Φxx(z)Φxyd(z)=Φxx(z)zNΦxs(z)E{[s(n+N)−s
(n+N)]2}min=2πj1∮C{Φss(z)−Hopt(z)z−NΦxs(z−1)}z−1dz
有物理约束
H
o
p
t
(
z
)
=
1
σ
w
2
B
(
z
)
[
Φ
x
y
d
(
z
)
B
(
z
−
1
)
]
+
=
1
σ
w
2
B
(
z
)
[
z
N
Φ
x
s
(
z
)
B
(
z
−
1
)
]
+
E
{
[
s
(
n
+
N
)
−
s
^
(
n
+
N
)
]
2
}
m
i
n
=
1
2
π
j
∮
C
{
Φ
s
s
(
z
)
−
H
o
p
t
(
z
)
z
−
N
Φ
x
s
(
z
−
1
)
}
z
−
1
d
z
H_{opt}(z) =\frac{1}{\sigma_w^2B(z)}\Big[\frac{\Phi_{xy_d}(z)}{B(z^{-1})}\Big]_+ =\frac{1}{\sigma_w^2B(z)}\Big[\frac{z^N\Phi_{xs}(z)}{B(z^{-1})}\Big]_+ \\ E\{[s(n+N)-\widehat{s}(n+N)]^2\}_{min}=\frac{1}{2\pi j}\oint_C \Big\{ \Phi_{ss}(z)-H_{opt}(z)z^{-N}\Phi_{xs}(z^{-1}) \Big\}z^{-1}dz
Hopt(z)=σw2B(z)1[B(z−1)Φxyd(z)]+=σw2B(z)1[B(z−1)zNΦxs(z)]+E{[s(n+N)−s
(n+N)]2}min=2πj1∮C{Φss(z)−Hopt(z)z−NΦxs(z−1)}z−1dz
2.4.3 N步纯预测器
我们假设纯预测是在
v
(
n
)
=
0
v(n)=0
v(n)=0 的情况下进行的预测，这时
x
(
n
)
=
s
(
n
)
x(n)=s(n)
x(n)=s(n) ，因而
Φ
x
x
(
z
)
=
Φ
s
s
(
z
)
=
Φ
x
s
=
σ
w
2
B
(
z
)
B
(
z
−
1
)
\Phi_{xx}(z)=\Phi_{ss}(z)=\Phi_{xs}=\sigma_w^2B(z)B(z^{-1})
Φxx(z)=Φss(z)=Φxs=σw2B(z)B(z−1)
在工程上，当噪声小到一定程度以后可以将噪声忽略或不再作更多考虑，为此，经常可以先用维纳滤波器等滤波器先将噪声抑制到一定水平之后再以纯预测方法作相应的预测处理。
对于因果系统
H
o
p
t
(
z
)
=
1
σ
w
2
B
(
z
)
[
Φ
x
y
d
(
z
)
B
(
z
−
1
)
]
+
=
1
σ
w
2
B
(
z
)
[
z
N
σ
w
2
B
(
z
)
B
(
z
−
1
)
B
(
z
−
1
)
]
+
=
1
B
(
z
)
[
z
N
B
(
z
)
]
+
\begin{aligned} H_{opt}(z)&=\frac{1}{\sigma_w^2B(z)} \Big[\frac{\Phi_{xy_d}(z)}{B(z^{-1})}\Big]_+ \\ &=\frac{1}{\sigma_w^2B(z)} \Big[\frac{z^N\sigma_w^2B(z)B(z^{-1})}{B(z^{-1})}\Big]_+\\ &=\frac{1}{B(z)} \Big[z^NB(z)\Big]_+\\ \end{aligned}
Hopt(z)=σw2B(z)1[B(z−1)Φxyd(z)]+=σw2B(z)1[B(z−1)zNσw2B(z)B(z−1)]+=B(z)1[zNB(z)]+
且
E
[
e
2
(
n
+
N
)
]
m
i
n
=
1
2
π
j
∮
C
{
Φ
s
s
(
z
)
−
H
o
p
t
(
z
)
Φ
x
y
d
(
z
−
1
)
}
z
−
1
d
z
E[e^2(n+N)]_{min}=\frac{1}{2\pi j}\oint_C \Big\{ \Phi_{ss}(z)-H_{opt}(z)\Phi_{xy_d}(z^{-1}) \Big\}z^{-1}dz
E[e2(n+N)]min=2πj1∮C{Φss(z)−Hopt(z)Φxyd(z−1)}z−1dz
考虑到
Φ
x
y
d
(
z
−
1
)
=
z
−
N
Φ
x
s
(
z
−
1
)
,
Φ
x
s
(
z
)
=
σ
w
2
B
(
z
)
B
(
z
−
1
)
=
Φ
x
s
(
z
−
1
)
\Phi_{xy_d}(z^{-1})=z^{-N}\Phi_{xs}(z^{-1}),\Phi_{xs}(z)=\sigma_w^2B(z)B(z^{-1})=\Phi_{xs}(z^{-1})
Φxyd(z−1)=z−NΦxs(z−1),Φxs(z)=σw2B(z)B(z−1)=Φxs(z−1) 故
E
[
e
2
(
n
+
N
)
]
m
i
n
=
1
2
π
j
∮
C
{
σ
w
2
B
(
z
)
B
(
z
−
1
)
−
σ
w
2
z
−
N
B
(
z
−
1
)
[
z
N
B
(
z
)
]
+
}
z
−
1
d
z
E[e^2(n+N)]_{min}=\frac{1}{2\pi j}\oint_C \Big\{ \sigma_w^2B(z)B(z^{-1})-\sigma_w^2z^{-N}B(z^{-1}) [z^{N}B(z)]_+ \Big\}z^{-1}dz
E[e2(n+N)]min=2πj1∮C{σw2B(z)B(z−1)−σw2z−NB(z−1)[zNB(z)]+}z−1dz
对于帕塞伐尔公式，当
x
(
n
)
,
y
(
n
)
x(n),y(n)
x(n),y(n) 为实序列时，有
∑
n
=
−
∞
∞
x
(
n
)
y
(
n
)
=
1
2
π
j
∮
c
X
(
z
)
Y
(
z
−
1
)
z
−
1
d
z
\sum_{n=-\infty}^\infty x(n)y(n)=\frac{1}{2\pi j}\oint_cX(z)Y(z^{-1})z^{-1}dz
n=−∞∑∞x(n)y(n)=2πj1∮cX(z)Y(z−1)z−1dz
由于
B
(
z
)
,
b
(
n
)
B(z),b(n)
B(z),b(n) 代表的是因果系统，所以
E
[
e
2
(
n
+
N
)
]
m
i
n
=
σ
w
2
∑
n
=
0
∞
b
2
(
n
)
−
σ
w
2
∑
n
=
0
∞
b
2
(
n
+
N
)
E[e^2(n+N)]_{min}=\sigma_w^2\sum_{n=0}^\infty b^2(n)- \sigma_w^2\sum_{n=0}^\infty b^2(n+N)
E[e2(n+N)]min=σw2n=0∑∞b2(n)−σw2n=0∑∞b2(n+N)
即
E
[
e
2
(
n
+
N
)
]
m
i
n
=
σ
w
2
∑
n
=
0
N
−
1
b
2
(
n
)
E[e^2(n+N)]_{min}=\sigma_w^2\sum_{n=0}^{N-1} b^2(n)
E[e2(n+N)]min=σw2n=0∑N−1b2(n)
这说明预测距离越长预测误差也就越大。
2.4.4 一步线性预测器的时域计算公式
我们对于一步线性预测器的讨论条件也是在无噪声的条件下进行的。在线性预测中，常常利用过去的
p
p
p 个样本
x
(
n
−
1
)
,
x
(
n
−
2
)
,
⋯
,
x
(
n
−
p
)
x(n-1),x(n-2),\cdots,x(n-p)
x(n−1),x(n−2),⋯,x(n−p) 来预测当前时刻
x
(
n
)
x(n)
x(n) 的问题。
y
(
n
)
=
x
^
(
n
)
=
∑
k
=
0
p
−
1
h
(
k
)
x
(
n
−
1
−
k
)
y(n)=\widehat{x}(n)=\sum_{k=0}^{p-1}h(k)x(n-1-k)
y(n)=x
(n)=k=0∑p−1h(k)x(n−1−k)
这是一个典型的
p
p
p 阶一步预测器，这里我们令
a
p
k
=
−
h
(
k
−
1
)
a_{pk}=-h(k-1)
apk=−h(k−1) 这可以表示为
x
^
(
n
)
=
−
∑
k
=
1
p
a
p
k
x
(
n
−
k
)
\widehat{x}(n)=-\sum_{k=1}^{p}a_{pk}x(n-k)
x
(n)=−k=1∑papkx(n−k)
所以
E
[
e
2
(
n
)
]
=
E
{
[
x
(
n
)
+
∑
k
=
1
p
a
p
k
x
(
n
−
k
)
]
2
}
∂
E
[
e
2
(
n
)
]
∂
a
p
l
=
0
E[e^2(n)]=E\Big\{ \big[ x(n)+\sum_{k=1}^{p}a_{pk}x(n-k) \big]^2 \Big\} \\ \frac{\partial E[e^2(n)]}{\partial a_{pl}}=0
E[e2(n)]=E{[x(n)+k=1∑papkx(n−k)]2}∂apl∂E[e2(n)]=0
故
2
E
{
[
x
(
n
)
+
∑
k
=
1
p
a
p
k
x
(
n
−
k
)
]
x
(
n
−
l
)
}
=
0
,
l
=
1
,
2
,
3
,
⋯
2E\Big\{ \big[ x(n)+\sum_{k=1}^{p}a_{pk}x(n-k) \big]x(n-l) \Big\}=0,\quad l=1,2,3,\cdots
2E{[x(n)+k=1∑papkx(n−k)]x(n−l)}=0,l=1,2,3,⋯
考虑到
x
^
(
n
)
=
−
∑
k
=
1
p
a
p
k
x
(
n
−
k
)
\widehat{x}(n)=-\sum_{k=1}^{p}a_{pk}x(n-k)
x
(n)=−∑k=1papkx(n−k) 所以上式可改写为
E
[
e
(
n
)
x
^
(
n
)
)
]
=
0
E
[
e
2
(
n
)
]
=
E
[
e
(
n
)
(
x
(
n
)
−
x
^
(
n
)
)
]
=
E
[
e
(
n
)
x
(
n
)
]
−
E
[
e
(
n
)
x
^
(
n
)
]
E[e(n)\widehat{x}(n))]=0\\ E[e^2(n)]=E[e(n)(x(n)-\widehat{x}(n))]=E[e(n)x(n)]-E[e(n)\widehat{x}(n)]
E[e(n)x
(n))]=0E[e2(n)]=E[e(n)(x(n)−x
(n))]=E[e(n)x(n)]−E[e(n)x
(n)]
故最小均方误差
E
[
e
2
(
n
)
]
m
i
n
=
E
[
e
(
n
)
x
(
n
)
]
=
E
{
[
x
(
n
)
+
∑
k
=
1
p
a
p
k
x
(
n
−
k
)
]
x
(
n
)
}
\begin{aligned} E[e^2(n)]_{min}&=E[e(n)x(n)] \\ &=E\Big\{ \big[ x(n)+\sum_{k=1}^{p}a_{pk}x(n-k) \big]x(n) \Big\} \end{aligned}
E[e2(n)]min=E[e(n)x(n)]=E{[x(n)+k=1∑papkx(n−k)]x(n)}
如果用自相关函数，并考虑到
ϕ
x
x
(
m
)
=
ϕ
x
x
(
−
m
)
\phi_{xx}(m)=\phi_{xx}(-m)
ϕxx(m)=ϕxx(−m) 等情况，则上式可表示为
ϕ
x
x
(
l
)
+
∑
k
=
1
p
a
p
k
ϕ
x
x
(
k
)
=
0
,
l
=
1
,
2
,
3
,
⋯
E
[
e
2
(
n
)
]
m
i
n
=
ϕ
x
x
(
0
)
+
∑
k
=
1
p
a
p
k
ϕ
x
x
(
k
)
\phi_{xx}(l)+\sum_{k=1}^p a_{pk}\phi_{xx}(k)=0,\quad l=1,2,3,\cdots \\ E[e^2(n)]_{min}=\phi_{xx}(0)+\sum_{k=1}^p a_{pk}\phi_{xx}(k)
ϕxx(l)+k=1∑papkϕxx(k)=0,l=1,2,3,⋯E[e2(n)]min=ϕxx(0)+k=1∑papkϕxx(k)
上述两个方程可以用矩阵形式表示
[
ϕ
x
x
(
0
)
ϕ
x
x
(
1
)
⋯
ϕ
x
x
(
p
)
ϕ
x
x
(
1
)
ϕ
x
x
(
0
)
⋯
ϕ
x
x
(
p
−
1
)
⋮
⋮
⋮
ϕ
x
x
(
p
)
ϕ
x
x
(
p
−
1
)
⋯
ϕ
x
x
(
0
)
]
[
1
a
p
1
⋮
a
p
p
]
=
[
E
[
e
2
(
n
)
]
m
i
n
0
⋮
0
]
\left[ \begin{matrix} \phi_{xx}(0) & \phi_{xx}(1) & \cdots & \phi_{xx}(p) \\ \phi_{xx}(1) & \phi_{xx}(0) & \cdots & \phi_{xx}(p-1) \\ \vdots & \vdots & & \vdots \\ \phi_{xx}(p) & \phi_{xx}(p-1) & \cdots & \phi_{xx}(0) \\ \end{matrix} \right] \left[ \begin{matrix} 1 \\ a_{p1} \\ \vdots \\ a_{pp} \\ \end{matrix} \right] = \left[ \begin{matrix} E[e^2(n)]_{min} \\ 0 \\ \vdots \\ 0 \\ \end{matrix} \right]
⎣⎢⎢⎢⎡ϕxx(0)ϕxx(1)⋮ϕxx(p)ϕxx(1)ϕxx(0)⋮ϕxx(p−1)⋯⋯⋯ϕxx(p)ϕxx(p−1)⋮ϕxx(0)⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡1ap1⋮app⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡E[e2(n)]min0⋮0⎦⎥⎥⎥⎤
常称为伊尔-沃克方程（Yule-Walker），这个方程只需要关注
x
(
n
)
x(n)
x(n) 的自相关函数，而它的自相关又是正定的，从而可以简便地解得
a
p
k
a_{pk}
apk 。