DSP 维纳滤波器的计算机实现.doc

维纳滤波器的计算机实现
一、 实验目的
1、 利用计算机编程实现加性噪声信号的维纳滤波。
2、 将计算机模拟实验结果与理论分析结果相比较，分析影响维纳滤波效果的各种因素，从而加深对维纳滤波的理解。
3、 利用维纳滤波一步纯预测方法实现对信号生成模型的参数估计。
二、 实验原理
1、 维纳滤波器是一种从噪声中提取信号的最佳线性估计方法，假定一个随机信号形式为：x(n)=s(n)+v(n),其中s(n)为有用信号，v(n)为噪声信号。而维纳滤波的作用就是让x(n)通过一个系统h(n)尽可能滤掉噪声，提取近似s(n),h(n)的选择以最小均方误差为准则。由维纳-霍夫方程知，只要求出φxx 及φxs就可求出h（h=φ-1xxφxs）。但要求h(n)满足因果性要求，维纳-霍夫方程便是一个难题，这里利用最佳FIR维纳滤波方法求解h(n)的近似，这也便于在计算机上实现，公式为：h =R-1xx rxs。实验中s(n)由信号生成模型：s(n)=as(n-1)+w(n)确定，其中a=0.95,w(n)是均值为0，方差为бw2=1的高斯白噪声，v(n)为均值为0，方差为1的高斯白噪声，且s(n)与v(n)不相关。实验中s(n)是已知的，但实际中如果s(n)已知，维纳滤波也就失去意义了，因此实验纯粹是为了理解维纳滤波原理而设计。
2、 维纳一步纯预测问题
S(n)的生成模型：s(n)+a1(n-1)+…+aps(n-p)=w(n),已知φxx（n）,利用Yule-walker方程即可得到信号生成模型参数ai(i=1,2…p)和
б2w 。
三、 实验步骤及结果分析
1、 根据维纳滤波原理绘制程序流程图
开始

输入样本个数L,FIR滤波器阶数N

产生L个v(n),w(n),s(n)和x(n)，
利用L个s(n)和x(n)，估计RSS和rxs

检验产生序列x(n的自相关和互相关函数
是否与理论值相符

N

在同一坐标内绘出x(n)自相关函数的
理论值和实际值
Y

在同一坐标内绘出最后100个s(n)和
x(n)。

调矩阵求逆子程序计算,将N个理想的h(n)和估计的h(n) 绘于同一坐标内

进行理想的维纳滤波得L个SI (n),和最后100个s(n)绘制于同一坐标

对x(n)进行过滤得L个SR(n)，和最后100个s(n)和绘于同一坐标内

结束
L个x(n),s(n), SI (n), SR(n),统计ex2,eI2,,eR2

2、 根据流程图编写程序（见附录1）并分析运行结果：
选择L=5000,N=10观察并记录、分析实验结果。
1） 与s(n)相比，信号x(n)在维纳滤波前后效果比较：

图1 图2
图1为维纳滤波后的s(n)与最后100个s(n)比较图
图2为未经维纳滤波的x(n)与最后100个s(n)比较图。
分析：显然与s(n)相比，x(n)在维纳滤波前与s(n)相差很大，维纳滤波后较接近s(n),可见滤波效果比较好。
2） 估计(n)与理想h(n)的比较：

图3
图3为估计(n)与理想h(n)的对比图。
分析：由图可见，二者近似程度除最后几个点外，其他近似度还是满高的，总体而言，近似效果不错。
3） 理想的维纳滤波与FIR维纳滤波效果对比：

图 4 图5
图4为理想维纳滤波效果，图5为FIR维纳滤波效果
分析：直接从图形观察，差异太小，无法观察其精度。只能通过最小均方差来比较其差异，结果为：理想维纳滤波ei= 0.2287,FIR维纳滤波ef=0.2254。可见，理想维纳滤波效果要好过FIR维纳滤波。
4） 自相关与互相关数据判断对效果的影响
分析：若去掉流程图中自相关与互相关数据判断步骤，可能会得到理想维纳滤波不如FIR滤波的效果，如其中一个结果：ei= 0.2503，ef= 0.2495。这里的判断步骤就是为了检测实际产生序列的自相关或互相关特性与理论值的近似程度，若误差很小且通过我们设定的某一下限则认为二者近似，所以最终的滤波效果才很近似。如果没有这里的判断，实际自相关或互相关则是任意的，完全有可能出现比理想维纳滤波更好的效果。
3、 固定L=5000,分别取N=3、20，根据实验结果，观察N的大小对(n)的估计和滤波效果的影响并记录实验结果。
实验结果：
图6为N=3时估计(n)与理想h(n)的对比图。
图7为N=20时估计(n)与理想h(n)的对比图。
图8为N=3的FIR滤波后所得(n)与实际S(n)后100位的比较图。
图9为N=20的FIR滤波后所得(n)与实际S(n)后100位的比较图。
其均方误差分别为：ei= 0.3175(N=3), 0.2500（N=20）
ef= 0.2762(N=3), 0.2488 (N=20)

图6 图7

图8 图9
分析：由图6、7可知，N的大小决定(n)与h(n)取值的个数，并通过观察并结合N=10的情况可知，N越大(n)与h(n)越接近。从最终均方误差的比较可知，N越大，滤波效果越好。
4、 固定N=10,改变L=1000、5000，根据实验结果，观察并记录L的大小对(n)的精度和滤波效果的影响。
实验结果：
图10为L=1000时估计(n)与理想h(n)的对比图。
图11为L=5000时估计(n)与理想h(n)的对比图。
图12为L=1000的FIR滤波后所得(n)与实际S(n)后100位的比较图。
图13为L=5000的FIR滤波后所得(n)与实际S(n)后100位的比较图。
其均方误差分别为：ei=0.2400 (L=1000), 0.2381（L=5000）
ef= 0.2390 (L=1000), 0.2375 (L=5000)

图10 图11

图12 图13
分析：由图10、11可知，L越大(n)与h(n)越接近，(n)的精度越高。由均方误差可知，L越大，滤波效果越高。这也容易理解，样本越大，精度自然越高。
5、 维纳一步纯预测
1） 画出信号生成模型参数估计的流程图

开始

输入信号生产模型的阶数p, AR模型的参数ai(i=1,2…p),бw2,信号s(n)的样本数L

利用randn函数产生L个w(n),并产生L个s(n)

利用Yule-Walker 方程，求出1….p

结束

2）根据流程图编写程序（见附录2）
3）运行信号生成模型程序，选择p=1,a1=-0.6,L=100.
理论值： бw2=1-a12=0.6400 a1=-0.6
估计值： w2= 0.9860 1= -0.5876
相对误差：error-бw2= -0.0140
error-a1= -0.0206
4)固定p=1,a1=-0.6, бw2=1,改变L=50、500，观察L的大小对信号生成模型参数估计精度的影响。
实验结果：
理论值： бw2=1-a12=0.6400 a1=-0.6
估计值：L=50 w2= 0.9730 1= -0.5760
估计值：L=500 w2= 0.9967 1= -0.5965
相对误差：L=50 error-бw2= -0.0270
error-a1= -0.0401
相对误差：L=500 error-бw2= -0.0033
error-a1= -0.0058
分析：显然样本个数L的增大，使得信号模型参数精度明显提高。
四、 实验总结
通过实验结果及分析可得出以下结论：
1、 样本个数越大，参数精度越高。
2、 影响维纳滤波效果的因素包括样本个数L、FIR滤波阶数，且均成正比关系。
3、 维纳一步纯预测，只要调整ai(1,2…p)即可实现最小均方误差。
五、 思考题答案
1、推导公式,验证式
推导：已知a= 0.95,w(n)为零均值方差为的高斯白噪声，v(n)是与s(n)互不相关的高斯白噪声，其均值为零，方差。
A(z)=,所有S(z)=W(z)A(z)=

2、由公式s(n)+(n-1)+…s(n-p)=w(n),怎样得到…和？
分析： 理论w(n)已知，即均值及已知，那么根据Yule-Walker方程有
,其中为(p+1)*(p+1)的s(n)自相关矩阵,A为（p+1）*1的系数列向量
及A=,而
由给出的理论，解方程即可得到估计值…；用估计值…代入方程即可得到估计值。
．
六、 源程序见附录1、2

附录1
clear all
L=input('请输入信号样本个数L=');
N=input('请输入FIR滤波器的阶数N=');
a = 0.95;
K = 50;
sigma_a2 = 1-a^2;
a_ = [ 1, -a];
while(1)% 利用randn()函数产生白噪声
wn = sqrt(sigma_a2)*( randn(L,1));
sn = filter(1, a_, wn);%H(Z)=1/(1-az_-1)
vn = randn(L,1);
xn = sn + vn;
r_xx = xcorr(xn,'unbiased'); %x(n)z自相关估计值
r_xx_t = a.^abs([-K:K]); % x(n)自相关理论值
r_xx_t(K+1)=r_xx_t(K+1)+1;
p = xcorr(sn,xn,'unbiased');%x(n)与s(n)互相关估计值
r_xs = p(L : L+K);
r_xs_t=a.^[0:K];%x(n)与s(n)互相关理论值
%检测实际值与理论值的近似程度
rou_xx = sum((r_xx(L-K:L+K)-r_xx_t').^2)/sum(r_xx_t.^2);
rou_xs = (sum(r_xs-r_xs_t').^2)/sum(r_xs_t.^2);
if rou_xx < 0.03 & rou_xs < 0.01
break;
end
end
%同一坐标绘制x(n)自相关函数理论值与实际值
figure(1),clf
stem(r_xx(L-K:L+K),'r')
hold on
stem(r_xx_t,'k')
title('r_xx的实际值（红色）及理论值（黑色）')
%同一坐标绘制x(n)与s(n)互相关函数理论值与实际值
figure(2),clf
stem(r_xs,'r')
hold on
stem(r_xs_t,'k')
title('r_xs的实际值（红色）及理论值（黑色）')
%同一坐标绘制最后100个s(n)与x(n)
figure(3),clf
stem(xn(L-99:L),'r')
hold on
stem(sn(L-99:L),'k')
title('最后100个s(n)（黑色）和x(n)（红色）');

%计算h(n)的估计值并与理想值比较
%构造x(n)的N阶自相关矩阵R_xx

n=0:N-1;
for i=1:N
for j=1:N
R_xx(i,j)=r_xx(i-j+L);
end
end
hopt=inv(R_xx)*r_xs(1:N);%利用维纳—霍夫方程求h
hopt_t=0.238*(0.724).^n; %理想h
%同一坐标绘制h(n)的实际值与理想值
figure(4),clf
stem(hopt,'-','r');
hold on
stem(hopt_t,'*','k');
title('h(n)估计值（红色），与真值（黑色）的比较');
%同一座标绘制理想维纳滤波后的s(n)和最后100个s(n)
sn_w=filter(hopt_t,1,xn);%理想维纳滤波
figure(5),clf
stem(sn_w(L-99:L),'*','r');
hold on
stem(sn(L-99:L),'.','k')
title('最后100个sn（黑色）和由维纳滤波器得到的sn_w（红色）的比较');
%同一坐标绘制实际维纳滤波后s(n)与最后100个s(n)
sn_f=filter(hopt,1,xn);
figure(6),clf
stem(sn_f(L-99:L),'*','r');
hold on
stem(sn(L-99:L),'.','k')
title('最后100个sn（黑色）和由FIR滤波器得到的sn_w（红色）的比较');
%求并比较各个均方差
e_x=sum((xn(1:L)-sn(1:L)).^2)/L
e_i=sum((sn_w(1:L)-sn(1:L)).^2)/L
e_f=sum((sn_f(1:L)-sn(1:L)).^2)/L

附录二
clear all;
%输入AR模型的参数
p=1;
a1=-0.6;
sigma_w2=1;
L=500;%s(m)的样本数
s=1;%统计5000次，最终估计值取平均
while s<=5000
wn=sqrt(sigma_w2)*randn(1,L);%w(n)为均值为0，方差为1的白噪声
sn=wn(1);%s(n)+a1s(n-1)=w(n),n=1时，s(n)=w(n)
for smp=2:L
sn=[sn,wn(smp)-a1*sn(smp-1)];%n>=2时，s(n)=w(n)-a1s(n-1)
end
r_ss=xcorr(sn,'unbiased');%求s(n)的自相关
est_a1(s)=-r_ss(L+1)/r_ss(L);%由r_ss(l)+a1r_ss(l-1)=0可得a1的估计值
%由r_ss(l)+a1r_ss(l-1)=sigma_w2可得sigma_w2的估计值
est_sigma(s)=r_ss(L)+est_a1(s)*r_ss(L-1);
s=s+1;
end
%估计值取统计平均值
estimate_a1=mean(est_a1)
estimate_sigma=mean(est_sigma)
%求估计值与理论的相对误差
error_a1=(estimate_a1-a1)/a1
error_sigma=(estimate_sigma-sigma_w2)/sigma_w2