维纳滤波器的应用

H1
噪声源

v ( n)
v2 (n)
维纳滤波器
R v2 h = rv1v2
假设v 假设v2(n)与d(n)不相关
R v2 h = rxv2
5
clear all; M=200; nn=randn(); v1(1)=nn; v2(1)=nn; d(1)=sin(0.05*pi); for ii=2:M nn=randn(); v1(ii)=0.8*v1(ii-1)+nn; v2(ii)=-0.6*v2(ii-1)+nn; d(ii)=sin(0.05*pi*ii); end x=d+v1; figure; plot(x);
14
X
D
15
应用：FIR维纳反卷积：MMSE均衡 应用：FIR维纳反卷积：MMSE均衡 维纳反卷积
x ( n) = d ( n) * g ( n ) + v ( n )
ˆ d ( n)
？
卷积失真：失焦的摄像机、运动模糊，频率选择性的通信信道 卷积失真：失焦的摄像机、运动模糊，
ˆ d ( n) = x ( n) * h( n)
信道 G(z) 不一定是最小相位的，且经常表示为FIR滤波器，而 不一定是最小相位的，且经常表示为FIR滤波器 滤波器， 均衡器 H(z) 也希望设计为FIR滤波器，即使逆滤波器G-1(z)存 也希望设计为FIR滤波器 即使逆滤波器G 滤波器， 在且是常态的，可能使得噪声被严重放大，导致显著的误差。 在且是常态的，可能使得噪声被严重放大，导致显著的误差。
确定滤波器阶数 ˆ 估计协方差矩阵 R
v2
P=32; N Re=zeros(1,P); 1 ˆ rxy (k ) = ∑+1 x(n)y* (n − k ) , k ≥ 0 N − k n=k for i=1:P k = i-1 ifor j=1:M n = j+i-1 j+iif j+i-1<M Re(1,i)=v2(1,j)*v2(1,j+i-1)'/(M-i+1)+Re(1,i); end end end R=toeplitz(Re.');
R x h = rdx
3
h=R r
−1 x dx
实际应用中的问题： 实际应用中的问题： R x , rdx 标准无偏估计 标准无偏估计： 估计：
？
x ( n) y * ( n − k ) , k ≥ 0
1 ˆ rxy (k ) = N −k
n = k +1
∑
N
N：用于估计的序列长度 最小均方误差： 最小均方误差：
10
维纳滤波器的基础：平稳条件 维纳滤波器的基础：
很多应用中都无法满足！ 很多应用中都无法满足！
11
令FIR维纳滤波器的阶数为 p=1, FIR维纳滤波器的阶数为
ˆ d (n) = hn d (n)
R x h = rdx
2 rdx (0) σd hn = = 2 rx (0) σ d + σ v2
估计被噪声污染的随机参数a P65页 P65页 例2： 估计被噪声污染的随机参数a a是高斯型： N (0, σ a ) 是高斯型： n
N (0, σ n )
进行N次观测： 进行N次观测： ri
= a + ni , i = 1, 2,L , N
2 σa 1 N ˆ ams (R ) = 2 ∑ Ri 2 σ a + σ n / N N i =1
i =1
可以利用训练序列的方法 盲估计（均衡）方法：采用训练序列并利用某种误差准则，或不需 盲估计（均衡）方法：采用训练序列并利用某种误差准则， 训练序列而基于信号的统计特性的方法。 训练序列而基于信号的统计特性的方法。
盲信号处理…… 盲信号处理……
17
I= imread('peppers.png'); figure;imshow(I); LEN = 31; THETA = 11; PSF = fspecial('motion',LEN,THETA); Blurred = imfilter(I,PSF,'circular','conv'); figure; imshow(Blurred); wnr1 = deconvwnr(Blurred,PSF); figure;imshow(wnr1);
8
rx (0) rx* (1) L rx* ( p − 1) rx (1) rx (0) L rx* ( p − 2) Rx = M M M rx (0) rx ( p − 1) rx ( p − 2) L
估计互相关矩阵 rxv ˆ 2
N 1 ˆ rxy (k ) = ∑+1 x(n)y* (n − k ) , k ≥ 0 N − k n=k
使得ξ 使得ξ最小的充分必要是相对 h*(k) 的偏导数为零
E{e(n) x* (n − k )} = 0 ,
2
k = 0,1,L , p − 1
由联合平稳条件： E{x (n − k ) x* (n − l )} = rx (l - k ) 由联合平稳条件：
且 rx (l ) = rx* ( −l )
s=conv(v2, Wopt.'); S=s(1:M); figure; plot(v1); axis([1 M -5 5]); hold on; plot(S,'r'); hold off figure; plot(x-S, 'LineWidth',2); hold on; plot(x,'r'); axis([1 M -5 5]);
检测、估计与调制理论
Detection, Estimation and Modulation Theory
维纳滤波器的应用
北京航空航天大学信号与信息处理教研室
1
FIR维纳滤波器－ FIR维纳滤波器－最佳抽头滤波器 维纳滤波器
x ( n) = d ( n ) + v ( n )
p −1
x(n)和d(n)联合平稳
16
ˆ d ( n ) = ∑ h( k ) x ( n − k )
k =0
p −1
使得如下均方误差最小 使得如下均方误差最小
R x h = rdx
G
ˆ ξ = E{| e(n) |2 } = E{| d (n) − d (n) |2 }
rdx (k ) = E{d (n) x* (n − k )} = ∑ g (i )rd (k + i )
设计一个FIR滤波器 设计一个FIR滤波器 h(k) , k=0,…p-1 ，对 d(n) 进行估计 =0,…p
ˆ d ( n ) = ∑ h( k ) x ( n − k )
k =0
使得如下均方误差最小 使得如下均方误差最小
ˆ ξ = E{| e(n) |2 } = E{| d (n) − d (n) |2 }
6
d (n) = sin(0.05π n)
v1 (n) = 0.8v1 (n − 1) + w(n)
v2 (n) = −0.6v2 (n − 1) + w(n)
w(n)是零均值、单位方差的白噪声 是零均值、
x(n) = d (n) + v1 (n)
d ( n)
v1 (n)
v2 (n)
7
x ( n)
r=zeros(1,P); for i=1:P for j=1:M if j-i+1>0 r(1,i)=x(1,j)*v2(1,j-i+1)'/(M-i+1)+r(1,i); end end end Wopt=inv(R)*r';
9
− 求滤波器系数 h = R v 1rxv 2 2
ˆ 滤波输出 d ( n) = x(n) − v2 ( n) * h( n)
* H ξ min = rd (0) − ∑ h(k )rdx (k ) = rd (0) − rdx R −1rdx x k =0 p −1
4
应用：基于维纳滤波器的噪声抑制 应用：
信号源
d ( n)
x( n) = d ( n) + v1 (n)
+
ˆ d ( n)
− v ( n) ˆ1
v1 (n)
18
12
实际又需要估计系数
σ x2 − σ v2 ˆ hn = 2 σx
ˆ µ= 1 NM
图像的维纳滤波器 wiener2.m
n1 , n2 ∈
∑η x(n , n )
1 2
ˆ σ x2 =
1 NM
n1 , n2 ∈
∑η
ˆ x 2 (n1 , n2 ) − µ 2
σ x2 − σ v2 ˆ (n , n ) = µ + ˆ ˆ ˆ d 1 2 ( x(n1 , n2 ) − µ ) 2 ˆ σx
ˆ ˆ σ v2 = mean(σ x2 )
13
RGB = imread('saturn.png'); I = rgb2gray(RGB); X = imnoise(I,'gaussian',0,0.005); D = wiener2(J,[5 5]); imshow(X, 'InitialMagnification','fit') figure, imshow(D, 'InitialMagnification','fit');
E{d (n) x* (n − l )} = rdx (l )
∑ h(k )r (l − k ) = r