1、维纳噪声抑制(例 6.6 的扩展实验)一、实验内容：假定图 6.8 中所需的信号 是一个正弦序列 , , ()dn0()sin)d0.5噪声序列 和 都是 AR(1) 过程，分别由如下的一阶差分方程产生：1()vn2110.8vng22()6()(其中 是零均值、单位方差的白噪声，与 不相关。()g d(a) 试用 Matlab 程序产生 x(n)和 的 500 个样点，画出波形图。2(b) 基于 x(n)和 的 500 个样点，设计 p 阶的最优 FIR 维纳滤波器，由2v估计 ，进而估计出 ，其中阶数 p 分别取为 p=3，6，9，12，2()v1试计算各种情况下估计 时的平均平方误差(均
2、方误差的样本估计，要叙1)述估计方案)，并画出对 d(n)估计的结果。(c) 有时辅助观测数据中也会漏入一些 d(n)信号，即辅助观测信号不仅是 ，2()vn而是 02v试针对 p=12 的情况，分别取几个不同的 值( 如 0.1, 0.5, 1.0)，研究这时的噪声抑制性能。(d) 若只有一路观测 的 1000 个样点，你能想办法近似完成对1()()xndn噪声 的有效抑制吗？试解释你的方法的基本原理，叙述你的实现方案。1()vn图 1 有辅观测数据的维纳噪声抑制器的原理图二、实验结果及分析问题一：试用 Matlab 程序产生 x(n)和 的 500 个样点，画出波形图。2v思路：噪声产生使
3、用 MATLAB 中的随机函数 Randn();结果分析：图 2、图 3 分别为观测信号 x(n)和辅助观测信号 的 500 个采样值，2()vn将用于估计 。其中 x(n)的相位取 0.25pi。2()vn1()vn0 50 100 150 200 250 300 350 400 450 500-6-4-202468估估估估x(n)0 50 100 150 200 250 300 350 400 450 500-4-3-2-1012345估估估估v2(n)图 2 观测信号 图 3 辅助观测信号问题二：基于 x(n)和 的 500 个样点，设计 p 阶的最优 FIR 维纳滤波器，由2v估计 ，
4、进而估计出 ，其中阶数 p 分别取为 p=3，6，9，12，试计2()v1()dn算各种情况下估计 时的平均平方误差(均方误差的样本估计，要叙述估计方1)案)，并画出对 d(n)估计的结果。思路辅助观测信号的噪声抑制可利用 Wiener-Hopf 方程 求出维纳滤波2xvvrwR器的系数 ，式中 是 的自相关阵， 是输入信号 和滤波器输入w2vR()n2xvr()n之间的互相关矢量；自相关序列的估计值采用公式2()vn得到；互相关序列的估计值采用公式221Nvnrkvk得到；估计 时的平均平方误差的样本估计为2 21xvnx1()vn，式中 N 为样本点数， 为噪声 的真实值，1Nnerv11
5、vn为噪声 的估计值；1v1结果分析：图 4图 7 是辅助观测信号的噪声抑制在最优 FIR 维纳滤波器的阶数 p分别为 3，6，9，12 的输出结果，其中蓝色为估计值，红色为真实值。由图可知，阶数 p 越高，平均平方误差就减小，FIR 的滤波效果就越好，d(n) 的估计值越接近于真实值。0 50 100 150 200 250 300 350 400 450 500-4-3-2-10123d(n)估估估d(n)估估估0 50 100 150 200 250 300 350 400 450 500-2.5-2-1.5-1-0.500.511.52d(n)估估估d(n)估估估图 4 维纳滤波器阶数
6、 p=3 图 5 维纳滤波器阶数 p=60 50 100 150 200 250 300 350 400 450 500-1.5-1-0.500.511.5d(n)估估估d(n)估估估0 50 100 150 200 250 300 350 400 450 500-1.5-1-0.500.511.5d(n)估估估d(n)估估估图 6 维纳滤波器阶数 p=9 图 7 维纳滤波器阶数 p=12表 1 为各种情况下估计 时的平均平方误差，从表 1 中可定量的知道，阶1()vn数 p 越高，FIR 的滤波效果就越好。这是因为随着 p 的增加，估计 时利用1()vn的数据越多，因此滤波效果越强。表 1
7、不同阶数下估计 时的平均平方误差1()vn阶数 p 平均平方误差3 0.60506 0.15809 0.042312 0.01233 4 5 6 7 8 9 10 11 1200.10.20.30.40.50.60.7图 8 不同阶数下估计 时的平均平方误差1()vn问题三：有时辅助观测数据中也会漏入一些 d(n)信号，即辅助观测信号不仅是 ，2()vn而是 试针对 p=12 的情况，分别取几个不同的 值(如 0.1, 0.5, 02()()vndn1.0)，研究这时的噪声抑制性能。思路：辅助观测信号漏入一些所需信号 d(n)时，维纳噪声抑制器的 Wiener-Hopf方程变为 ，式中 是 的
8、自相关阵， 是输入信号 和滤波器0xvvrwR0vR()n0xvr()xn输入 之间的互相关矢量；0()n结果分析：图 9图 11 是当辅助观测信号漏入一些所需信号 d(n)时候维纳噪声抑制器的输出结果，由图可知，泄漏系数 a 越大，平均平方误差越大，滤波效果越差，d(n)的估计值离真实值越远。0 50 100 150 200 250 300 350 400 450 500-1-0.8-0.6-0.4-0.200.20.40.60.81d(n)估估估d(n)估估估0 50 100 150 200 250 300 350 400 450 500-4-3-2-101234d(n)估估估d(n)估估
9、估图 9 泄漏系数 a=0.1 图 10 泄漏系数 a=0.50 50 100 150 200 250 300 350 400 450 500-2-1.5-1-0.500.511.522.53d(n)估估估d(n)估估估0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.511.522.5图 11 泄漏系数 a=1.0 图 12 不同的泄漏系数对应的平均平方误差表 2 为不同的泄漏系数对应的平均平方误差，从表 2 中可定量的知道，泄漏系数 a越大，滤波效果越差。因为泄漏系数 a 的增加意味着辅助观测信号中有用信号 d(n)的成份越多，这样在维纳滤波器的输出端含有的有
10、用信号 d(n)的成份也就越多，而 d(n)的估计值， 为维纳滤波器的输出端信号，相应地减去了更多的有用信号 d(n)，1()dnxvn1从而导致 d(n)的估计值离真实值越远。表 2 不同的泄漏系数对应的平均平方误差泄露系数 a 平均平方误差0.1 0.24320.5 1.77881.0 2.4148问题四：若只有一路观测 的 1000 个样点，你能想办法近似完成1()()xndvn对噪声 的有效抑制吗？试解释你的方法的基本原理，叙述你的实现方案。1()vn图 13 有辅观测数据的维纳噪声抑制器的原理图思路：没有辅助观测信号时，可利用信号的延时去掉噪声序列的相关性，从而进行维纳滤波，估计出所
11、需信号 d(n)。结果分析：当只有一路观测 的 1000 个样点，能近似完成对噪声1()()xndvn的有效抑制，其原理框图见图 14。其基本原理是 为窄带过程，)(1nv )(nd为宽带过程，对延迟 ， 的自相关序列近似为 0，另外假设Lk)(1和 互不相关，求解框图中维纳滤波器的系数 的 Wiener-Hopf 方程为)(d1v w，式中 为 的自相关阵， 是 和 之间的互相关矢量。dyrwRyR)(ndyr)(n由于 11 111 *ydvxrkEykdnLvnkLdnvnLkdEkEvvrk11 *dydrkEdnykLvnkLkEdvnkLrkL当延迟 时， 的自相关序列近似为 0，
12、所以 ，Lk)(1nv )()()rrxddy Wiener-Hopf 方程变为 ，式中 为 的自相关阵， 为以 为起LxrwRxR)(nL始值的自相关矢量。由于 可由一路观测 的 1000 个样点进行估计，因)(k此维纳滤波器的系数 可以求出，即能近似完成对噪声 的有效抑制。)(1nv维纳滤波器 W ( z )1xndvnLzynxLdn图 14 有辅观测数据的维纳噪声抑制器的原理图附录：MATLAB 代码clear all;clc;N=500;n=1:N;n=n;d=sin(0.05*pi*n+0.25*pi);g=randn(N,1);v1=filter(1,1,-0.8,g);v2=f
13、ilter(1,1,0.6,g);x=d+v1;figure(1);plot(n,x);ylabel(观察信号 x(n);grid on;figure(2);plot(n,v2);ylabel(噪声信号 v2(n);grid on;px=3 6 9 12;err=;for l=1:length(px)p=px(l);rv2=rxk2(v2,p);rxv2=rxyk(x,v2,p);tp=zeros(2*p+1,1);tp(p+1:2*p+1)=rv2(1:p+1);for i=1:p;tp(i)=rv2(p+2-i);endRv2=zeros(p+1,p+1);for i=1:p+1;Rv2(
14、i,:)=tp(p+2-i:2*p+2-i);endw=inv(Rv2)*rxv2(1:p+1);v1g=filter(w,1,v2);dg=x-v1g;figure(2+l);plot(n,dg,b-,n,d,r-);legend(d(n)估计值,d(n)真实值);grid on;errg=sum(v1-v1g).2)/N;err=err,errg;end;figure(7);stem(px,err);ax=0.1 0.5 1.0;errx=;for l=1:length(ax)p=12;v0=v2+ax(l)*d;rv2x=rxk2(v0,p);rxv2x=rxyk(x,v0,p);tpx
15、=zeros(2*p+1,1);tpx(p+1:2*p+1)=rv2x(1:p+1);for i=1:p;tpx(i)=rv2x(p+2-i);endRv2x=zeros(p+1,p+1);for i=1:p+1;Rv2x(i,:)=tpx(p+2-i:2*p+2-i);endw=inv(Rv2x)*rxv2x(1:p+1);v1gx=filter(w,1,v0);dgx=x-v1gx;figure(7+l);plot(n,dgx,b-,n,d,r-);legend(d(n)估计值,d(n)真实值);grid on;errgx=sum(v1-v1gx).2)/N;errx=errx,errgx;end;figure(11)stem(ax,errx);
