数字信号处理实验——维纳滤波器设计..

实验一 维纳滤波
1. 实验内容
设计一个维纳滤波器：
（1） 产生三组观测数据，首先根据()(1)()s n as n w n =-+产生信号()s n ，将其加噪，（信噪比分别为20,10,6dB dB dB ），得到观测数据123(),(),()x n x n x n 。

（2） 估计()i x n ，1,2,3i =的AR 模型参数。

假设信号长度为L ，AR 模型阶数为N ，分析实验结果，并讨论改变L ，N 对实验结果的影响。

2. 实验原理
滤波目的是从被噪声污染的信号中分离出有用的信号来，最大限度地抑制噪声。

对信号进行滤波的实质就是对信号进行估计。

滤波问题就是设计一个线性滤波器，使得滤波器的输出信号()y n 是期望响应()s n 的一个估计值。

下图就是观测信号的组成和信号滤波的一般模型。

观测信号()()()x n s n v n =+ 信号滤波的一般模型
维纳滤波解决从噪声中提取信号的滤波问题，并以估计的结果与真值之间的误差均方值最小作为最佳准则。

它根据()()(),1,
,x n x n x n m --估计信号的当前
值，它的解以系统的系统函数()H z 或单位脉冲()h n 形式给出，这种系统常称为最佳线性滤波器。

维纳滤波器设计的任务就是选择()h n ，使其输出信号()y n 与期望信号()d n 误差的均方值最小。

假设滤波系统()h n 是一个线性时不变系统，它的()h n 和输入信号都是复函数，设
()()()h n a n jb n =+ 0,1,
n
=
考虑系统的因果性，可得到滤波器的输出
()()()()()0
*m y n h n x n h m x n m +∞
===-∑ 0,1,
n
=
设期望信号()d n ，误差信号()e n 及其均方误差()2
E e n ⎡⎤⎣⎦
分别为
()()()()()e n d n y n s n y n =-=-
()()()()()()22
2
0m E e n E d n y n E d n h m x n m ∞=⎡⎤
⎡⎤⎡⎤=-=--⎢⎥⎣⎦⎣
⎦⎢⎥⎣⎦
∑ 要使均方误差为最小，需满足：
()()
2
0E e n h j ⎡⎤∂⎣⎦=∂ 整理得()()0E x n j e n *⎡⎤-=⎣⎦，等价于()()0E x n j e n *
⎡⎤-=⎣⎦
上式说明，均方误差达到最小值的充要条件使误差信号与任一进入估计的输入信号正交，这就是正交性原理。

将()()0E x n j e n *
⎡⎤-=⎣⎦展开，得
()()()()00m E x n k d n h m x m +∞
***
=⎡⎤⎛⎫--=⎢⎥ ⎪⎝⎭⎣
⎦∑
整理得 ()()()0
dx xx m r k h m r m k +∞
*=-=-∑ 0,1,2,
k
= 等价于()()()()()0
dx xx xx m r k h m r k m h k r k +∞
==-=*∑ 0,1,2,
k
=
此式称为维纳-霍夫（Wiener-Holf ）方程。

解此方程可得到最优权系数
012,,,
h h h ，此式是Wiener 滤波器的一般方程。

定义
121M M h h h h ⨯⎡⎤
⎢⎥
⎢⎥=⎢⎥
⎢⎥
⎣⎦()()()1011xd xd xd xd M r r R r M ⨯⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥-⎢⎥⎣⎦
()()()()()()()()
()011102120xx xx xx xx
xx xx xx xx xx xx M M
r r r M r r r M R r M r M r *
*
*
⨯⎡⎤
-⎢
⎥-⎢⎥=⎢⎥⎢⎥
--⎢⎥⎣⎦
则维纳-霍夫方程可写成矩阵形式
xd xx R R h =求逆，得1
xx xd h R R -=
此式表明，已知期望信号与观测数据的互相关函数及观测数据的自相关函数时，可以通过矩阵求逆运算，得到维纳滤波器的最佳解。

3. 实验结果及分析
（1）当L=200,N=6
信噪比为20dB 的滤波效果
信噪比为10dB的滤波效果
信噪比为6dB的滤波效果（2）当L=200,N=60
信噪比为20dB的滤波效果
信噪比为10dB的滤波效果
信噪比为6dB的滤波效果（3）当L=600,N=6
信噪比为20dB的滤波效果
信噪比为10dB的滤波效果
信噪比为6dB的滤波效果
实验分析：别取信号长度为200、600，滤波器长度为6、60，加噪信噪比为20dB、10dB、6dB，组合进行实验。

每组实验得到的最小均方误差统计如下表。

由此表可以看出，信号长度越长，最小均方误差（绝对值）越大，精度越差；在信噪比较大（误差影响较小）的滤波过程中，滤波器长度约长，最小均方误差（绝对值）越小，精度越好。

表1 最小均方误差统计表
对于相同信号和滤波器（这里取L=200,N=6），信噪比越大，最小距离误差约小；而当信噪比较小时，信号与噪声值接近，导致滤波效果受到影响，最小距离误差变大。

4. 源代码
clear;
clc;
%初始化变量
L=200;
%信号长度
N=6;
%滤波器的阶次
a=0.96;
wn=randn(L,1);
%wn为用于生成信号的噪声信号，随机生成一个L*1矩阵，生成矩阵的元素值在%区间(0.0,1.0)之间
sn=zeros(L,1);
%sn为信号，生成一个L*1的零矩阵
hn=zeros(N,1);
%hn为系统单位脉冲响应生成一个N*1的零矩阵
rxx=zeros(N,1);
%rxx为自相关函数，生成一个N*1的零矩阵
rxd=zeros(N,1);
%rxd为互相关函数，生成一个N*1的零矩阵
yn=zeros(L,1);
%yn为输出信号，生成一个L*1的零矩阵
xt=zeros(L+N,1);
%生成一个(L+N)*1的零矩阵
gn=zeros(L,1);
%gn为yn与sn最小距离误差信号，生成一个L*1的零矩阵
%根据给定公式s(n)=as(n-1)+w(n),生成理想信号
for i=2:L
sn(i,1)=a*sn(i-1,1)+wn(i,1);
end
sn(1,1)=wn(1,1);
subplot(2,2,1);
plot(sn,'r'),axis([0 200 -10 10]),xlabel('时间'),ylabel('幅度'),title('sn')
%生成期望信号方差cd
cd=(var(wn))/(1-a^2);
%对信号加噪
x1=awgn(sn,20);
x2=awgn(sn,10);
x3=awgn(sn,6);
subplot(2,2,2)
plot(x3,'g'),axis([0 200 -10 10]),xlabel('时间'),ylabel('幅度'),title('x3');
%生成输入信号与理想信号的互相关函数,此处x1为输入信号，sn为期望信号for i=1:N,
for m=i:1:L,
rxd(i,1)=rxd(i,1)+x3(m,1)*sn(m-i+1,1);
end
end
%生成输入信号的自相关函数
for i=1:N,
for m=i:1:L,
rxx(i,1)=rxx(i,1)+x3(m,1)*x3(m-i+1,1);
end
end
%将自相关函数生成托普勒斯矩阵
rxx1=toeplitz(rxx);
%生成逆矩阵
irxx=inv(rxx1);
%生成滤波器系数h(n)
hn=irxx*rxd;
for i=1:L
xt(i+N,1)=x3(i,1);
end
%实现滤波
for i=1:L,
for m=1:N,
yn(i,1)=yn(i,1)+xt(i+N+1-m,1)*hn(m,1);
end
end
%计算最小均方误差信号en
en=0;
en=cd-(rxd')*hn;
%生成最小距离误差信号gn
gn=yn-sn;
%画出滤波后的信号时域图
subplot(2,2,3);
plot(yn),axis([0 200 -10 10]),xlabel('时间'),ylabel('幅度'),title('yn');
%画出理想信号与输出信号对比图
subplot(2,2,4);
plot(sn,'r'),axis([0 200 -10 10]),xlabel('时间'),ylabel('幅度'),title('sn与yn对比');
hold on;
plot(yn,'b'),axis([0 200 -10 10]);
hold off;
%画出最小距离误差信号图
figure;
plot(gn),axis([0 200 -2 2]),xlabel('时间'),ylabel('幅度'),title('gn');。