维纳滤波器
1.wiener filter
维纳滤波器是一个均方误差最小准则下的最优滤波器,定义一个离散线性时不变系统,输入 x ( n ) x(n) x(n)、输出 y ^ ( n ) \hat y(n) y^(n),参考信号 d ( n ) d(n) d(n)
误差信号
e
(
n
)
=
y
^
(
n
)
−
d
(
n
)
(1)
e(n)=\hat y(n)-d(n)\tag1
e(n)=y^(n)−d(n)(1)
定义长度
M
M
M的FIR滤波器
h
\bold h
h如下
h
=
[
h
0
,
h
1
.
.
.
.
.
h
M
−
1
]
T
y
^
(
n
)
=
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
(2)
\begin{aligned} \bold h&=[h_0,h_1.....h_{M-1}]^T\\ \hat y(n)&=\sum_{m=0}^{M-1}{x(n-m)*h(m)} \end{aligned}\tag2
hy^(n)=[h0,h1.....hM−1]T=m=0∑M−1x(n−m)∗h(m)(2)
维纳滤波器的目标寻找一组滤波器系数使误差信号的均方误差最小
h
=
a
r
g
m
i
n
E
[
e
2
(
n
)
]
(3)
\bold h=argminE[e^2(n)]\tag3
h=argminE[e2(n)](3)
其中E为数学期望,将(3)式展开
E
[
e
2
(
n
)
]
=
E
[
(
y
^
(
n
)
−
d
(
n
)
2
]
=
E
[
y
^
2
(
n
)
]
−
2
∗
E
[
y
^
(
n
)
∗
d
(
n
)
]
+
E
[
d
2
(
n
)
]
(4)
\begin{aligned} E[e^2(n)]&=E[(\hat y(n)-d(n)^2]\\ &=E[\hat y^2(n)]-2*E[\hat y(n)*d(n)]+E[d^2(n)] \end{aligned}\tag4
E[e2(n)]=E[(y^(n)−d(n)2]=E[y^2(n)]−2∗E[y^(n)∗d(n)]+E[d2(n)](4)
代入
y
^
(
n
)
\hat y(n)
y^(n),继续
E
[
e
2
(
n
)
]
=
E
[
(
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
)
2
]
−
2
∗
E
[
(
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
)
∗
d
(
n
)
]
+
E
[
d
2
(
n
)
]
,
(5)
\begin{aligned} E[e^2(n)]=&E[(\sum_{m=0}^{M-1}{x(n-m)*h(m)})^2]-\\ &2*E[(\sum_{m=0}^{M-1}{x(n-m)*h(m)})*d(n)]+\\ &E[d^2(n)],\\ \end{aligned}\tag5
E[e2(n)]=E[(m=0∑M−1x(n−m)∗h(m))2]−2∗E[(m=0∑M−1x(n−m)∗h(m))∗d(n)]+E[d2(n)],(5)
误差函数对各权值求偏微分,这个复合函数求导可看这里
Δ
E
[
e
2
(
n
)
]
h
i
=
2
E
[
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
∗
x
(
n
−
i
)
]
−
2
E
[
x
(
n
−
i
)
d
(
n
)
]
,
=
2
∑
m
=
0
M
−
1
E
[
x
(
n
−
m
)
∗
x
(
n
−
i
)
]
∗
h
(
m
)
−
2
E
[
x
(
n
−
i
)
d
(
n
)
]
,
i
=
0
,
1....
M
−
1
(6)
\begin{aligned} \frac{\Delta E[e^2(n)]}{h_i}=2&E[\sum_{m=0}^{M-1}{x(n-m)*h(m)}*x(n-i)]-\\2&E[x(n-i)d(n)],\\ =2&\sum_{m=0}^{M-1}E[x(n-m)*x(n-i)]*h(m)-\\2&E[x(n-i)d(n)], i=0,1....M-1\\ \end{aligned}\tag6
hiΔE[e2(n)]=22=22E[m=0∑M−1x(n−m)∗h(m)∗x(n−i)]−E[x(n−i)d(n)],m=0∑M−1E[x(n−m)∗x(n−i)]∗h(m)−E[x(n−i)d(n)],i=0,1....M−1(6)
这里假设
x
(
n
)
和
d
(
n
)
x(n)和d(n)
x(n)和d(n)为平稳随机过程且联合平稳,则有自相关函数
r
x
x
r_{xx}
rxx和互相关函数
r
x
d
r_{xd}
rxd如下
r
x
x
(
m
)
=
E
[
x
(
n
)
x
(
n
−
m
)
]
r
x
d
(
m
)
=
E
[
x
(
n
)
d
(
n
−
m
)
]
(7)
\begin{aligned} r_{xx}(m)&=E[x(n)x(n-m)]\\ r_{xd}(m)&=E[x(n)d(n-m)] \end{aligned}\tag7
rxx(m)rxd(m)=E[x(n)x(n−m)]=E[x(n)d(n−m)](7)
上式中m为第二个信号相对第一个信号的延时,平稳随机信号的相关函数只与信号的延时有关,与信号的起始位置无关,则(6)使可以去掉期望符号重写成
Δ
E
[
e
2
(
n
)
]
h
i
=
2
∑
m
=
0
M
−
1
r
x
x
(
i
−
m
)
∗
h
(
m
)
−
2
r
x
d
(
−
i
)
,
i
=
0
,
1....
M
−
1
(8)
\begin{aligned} \frac{\Delta E[e^2(n)]}{h_i}=2&\sum_{m=0}^{M-1}r_{xx}(i-m)*h(m)-\\ 2&r_{xd}(-i), i=0,1....M-1\\ \end{aligned}\tag8
hiΔE[e2(n)]=22m=0∑M−1rxx(i−m)∗h(m)−rxd(−i),i=0,1....M−1(8)
当均方误差最小时,(8)式中的导数等于0,
∑
m
=
0
M
−
1
r
x
x
(
i
−
m
)
∗
h
(
m
)
=
r
x
d
(
−
i
)
,
i
=
0
,
1....
M
−
1
(9)
\sum_{m=0}^{M-1}r_{xx}(i-m)*h(m)=r_{xd}(-i), i=0,1....M-1\tag9
m=0∑M−1rxx(i−m)∗h(m)=rxd(−i),i=0,1....M−1(9)
(9)式就是Wiener–Hopf equations,注意,写这个式子的时候要注意与相关函数的定义对应起来,有的地方定义
r
x
d
(
m
)
=
E
[
x
(
n
)
x
(
n
−
m
)
]
r_{xd}(m)=E[x(n)x(n-m)]
rxd(m)=E[x(n)x(n−m)],或者
r
d
x
(
m
)
=
E
[
x
(
n
)
x
(
n
+
m
)
]
r_{dx}(m)=E[x(n)x(n+m)]
rdx(m)=E[x(n)x(n+m)],自相关函数是偶对称的(
r
x
x
(
i
−
j
)
=
r
x
x
(
j
−
i
)
r_{xx}(i-j)=r_{xx}(j-i)
rxx(i−j)=rxx(j−i)),而互相关并不对称,
r
x
d
(
k
)
=
r
d
x
(
−
k
)
r_{xd}(k)=r_{dx}(-k)
rxd(k)=rdx(−k),因此要注意(9)式的符号。
继续,将(9)式写成矩阵的形式
R
∗
h
=
r
d
x
(10)
\mathbf{R}*\mathbf{h}=\mathbf{r_{dx}}\tag{10}
R∗h=rdx(10)
展开表示为
[
r
(
0
)
r
(
1
)
⋯
r
(
M
−
1
)
r
∗
(
1
)
r
(
0
)
⋯
r
(
M
−
2
)
⋮
⋮
⋱
⋮
r
∗
(
M
−
1
)
r
∗
(
M
−
2
)
⋯
r
(
0
)
]
∗
[
h
0
h
1
⋮
h
M
−
1
]
=
[
r
d
x
(
0
)
r
d
x
(
1
)
⋮
r
d
x
(
M
−
1
)
]
(11)
\left[\begin{array}{cccc}{r(0)} & {r(1)} & {\cdots} & {r(M-1)} \\ {r^{*}(1)} & {r(0)} & {\cdots} & {r(M-2)} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {r^{*}(M-1)} & {r^{*}(M-2)} & {\cdots} & {r(0)}\end{array}\right]* \begin{bmatrix} h_0\\ h_1\\ {\vdots}\\ h_{M-1} \end{bmatrix}= \begin{bmatrix} r_{dx}(0)\\ r_{dx}(1)\\ {\vdots}\\ r_{dx}(M-1) \end{bmatrix}\tag{11}
⎣⎢⎢⎢⎡r(0)r∗(1)⋮r∗(M−1)r(1)r(0)⋮r∗(M−2)⋯⋯⋱⋯r(M−1)r(M−2)⋮r(0)⎦⎥⎥⎥⎤∗⎣⎢⎢⎢⎡h0h1⋮hM−1⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡rdx(0)rdx(1)⋮rdx(M−1)⎦⎥⎥⎥⎤(11)
这样,维纳滤波在时域的闭式解就为
h
o
p
t
=
R
−
1
∗
r
d
x
(12)
\mathbf{h_{opt}}=\mathbf{R}^{-1}*\mathbf{r_{dx}}\tag{12}
hopt=R−1∗rdx(12)
计算这个式子需要观测信号的
M
∗
M
M*M
M∗M自相关矩阵跟观测信号和期望信号
M
∗
1
M*1
M∗1的互相关函数
同样,将最优解代入到误差函数中,可以得到最小均方误差如下
E
[
e
2
(
n
)
]
m
i
n
=
σ
d
2
−
r
d
x
H
h
o
p
t
=
σ
d
2
−
r
d
x
H
R
−
1
r
d
x
\begin{aligned} E[e^2(n)]_{min} &= \sigma_d^2-\mathbf{r_{dx}}^H\mathbf{h_{opt}}\\ &=\sigma_d^2-\mathbf{r_{dx}}^H\mathbf{R}^{-1}\mathbf{r_{dx}} \end{aligned}
E[e2(n)]min=σd2−rdxHhopt=σd2−rdxHR−1rdx
2.code
在MATLAB里验证一下维纳滤波器闭式解
% clear all
close all
T = 100000;
x = randn(T,1);
%
h = [1.2,4.7,3.6,5.5,0.8]';
M = length(h);
d = filter(h,1,x);
[r,lag1] = xcorr(x,'unbiased');
R = zeros(5,5);
% 利用相关函数的值计算相关矩阵
for i = 1:5
for j = 1:5
R(i,j) = r(find(lag1==(i-j)));
end
end
[p,lag2] = xcorr(x,d,'unbiased');
P = zeros(M,1);
for i = 0:M-1
P(i+1) = p(find(lag2==(i*-1)));
end
h
h_opt = R\P
e_min = var(d)-P'*h_opt
figure,stem(h)
hold on,plot(h_opt)
legend('filter weighs','estimated weights');
从结果上看,闭式解求得的权值和实际权值一致。代码里调用了matlab的xcorr函数,在取互相关函数值的位置(lag)时要,注意互相关函数的定义与闭式解的形式一致
3.求导计算
这里记录一下误差函数的求导
E
[
e
2
(
n
)
]
=
E
[
(
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
)
2
]
−
2
∗
E
[
(
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
)
∗
d
(
n
)
]
+
E
[
d
2
(
n
)
]
,
(13)
\begin{aligned} E[e^2(n)]=&E[(\sum_{m=0}^{M-1}{x(n-m)*h(m)})^2]-\\ &2*E[(\sum_{m=0}^{M-1}{x(n-m)*h(m)})*d(n)]+\\ &E[d^2(n)],\\ \end{aligned}\tag{13}
E[e2(n)]=E[(m=0∑M−1x(n−m)∗h(m))2]−2∗E[(m=0∑M−1x(n−m)∗h(m))∗d(n)]+E[d2(n)],(13)
这里就只看第一个求和式子
f
(
n
)
=
E
[
(
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
)
2
]
=
E
[
g
2
(
n
)
]
\begin{aligned} f(n)& = E[(\sum_{m=0}^{M-1}{x(n-m)*h(m)})^2]\\ &=E[g^2(n)] \end{aligned}
f(n)=E[(m=0∑M−1x(n−m)∗h(m))2]=E[g2(n)]
要计算
∂
f
(
n
)
∂
h
(
m
)
\frac{\partial f(n)}{\partial h(m)}
∂h(m)∂f(n),
∂
f
(
n
)
∂
h
(
m
)
=
E
[
2
∗
g
(
n
)
∗
∂
g
(
n
)
∂
h
(
m
)
]
(14)
\begin{aligned} \frac{\partial f(n)}{\partial h(m)}&=E[2*g(n)*\frac{\partial g(n)}{\partial h(m)}]\\ \end{aligned}\tag{14}
∂h(m)∂f(n)=E[2∗g(n)∗∂h(m)∂g(n)](14)
再计算
∂
g
(
n
)
∂
h
(
m
)
\frac{\partial g(n)}{\partial h(m)}
∂h(m)∂g(n)
g
(
n
)
=
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
h
(
m
)
=
x
(
n
)
h
(
0
)
+
x
(
n
−
1
)
h
(
1
)
⋯
x
(
0
)
h
(
M
−
1
)
∂
g
(
n
)
∂
h
(
0
)
=
x
(
n
)
∂
g
(
n
)
∂
h
(
1
)
=
x
(
n
−
1
)
⋮
∂
g
(
n
)
∂
h
(
i
)
=
x
(
n
−
i
)
⋮
∂
g
(
n
)
∂
h
(
M
−
1
)
=
x
(
n
−
M
+
1
)
\begin{aligned} g(n)&=\sum_{m=0}^{M-1}{x(n-m)*h(m)} \\ &=x(n)h(0)+x(n-1)h(1) \cdots x(0)h(M-1)\\ \frac{\partial g(n)}{\partial h(0)}&=x(n)\\ \frac{\partial g(n)}{\partial h(1)}&=x(n-1)\\ \vdots\\ \frac{\partial g(n)}{\partial h(i)}&=x(n-i)\\ \vdots\\ \frac{\partial g(n)}{\partial h(M-1)}&=x(n-M+1)\\ \end{aligned}
g(n)∂h(0)∂g(n)∂h(1)∂g(n)⋮∂h(i)∂g(n)⋮∂h(M−1)∂g(n)=m=0∑M−1x(n−m)∗h(m)=x(n)h(0)+x(n−1)h(1)⋯x(0)h(M−1)=x(n)=x(n−1)=x(n−i)=x(n−M+1)
所以
∂
g
(
n
)
∂
h
(
m
)
=
x
(
n
−
m
)
(15)
\frac{\partial g(n)}{\partial h(m)}=x(n-m)\tag{15}
∂h(m)∂g(n)=x(n−m)(15)
这样将(15)带入到(14)中,得到
∂
f
(
n
)
∂
h
(
i
)
=
E
[
2
∗
g
(
n
)
∗
∂
g
(
n
)
∂
h
(
i
)
]
=
2
E
[
g
(
n
)
x
(
n
−
i
)
]
=
2
E
[
∑
m
=
0
M
−
1
x
(
n
−
m
)
∗
x
(
m
)
x
(
n
−
i
)
]
(16)
\begin{aligned} \frac{\partial f(n)}{\partial h(i)}&=E[2*g(n)*\frac{\partial g(n)}{\partial h(i)}]\\ &=2E[g(n)x(n-i)]\\ &=2E[\sum_{m=0}^{M-1}{x(n-m)*x(m)x(n-i)]} \end{aligned}\tag{16}
∂h(i)∂f(n)=E[2∗g(n)∗∂h(i)∂g(n)]=2E[g(n)x(n−i)]=2E[m=0∑M−1x(n−m)∗x(m)x(n−i)](16)
这样就完成了这个复合函数的求导,都展开写看其实也挺简单,主要是期望符号和求和符号有些容易混淆
Reference:
1.《数字信号处理:时域离散随机信号处理》
2.《自适应滤波器原理》
3. https://en.wikipedia.org/wiki/Wiener_filter