维纳滤波器实验

实验一 维纳滤波器的计算机实现 一、实验目的 1．利用计算机编程实现加性噪声信号的维纳滤波。 2．将计算机模拟实验结果与理论分析结果相比较，分析影响维纳滤波效果的各 种因素，从 而加深对维纳滤波的理解。 利用维纳一步纯预测方法实现对信号生成模型的参数估计。 二、实验原理 维纳滤波是一种从噪声背景中提取信号的最佳线性滤波方法，假定一个随 机信号 x(n)具有以下形式： ( ) ( ) ( )nvnsnx += （1-1） 其中，s(n) 为有用信号，v(n) 为噪声干扰，将其输入一个单位脉冲响应为 h(n) 的 线性系统，其输出为： ( ) ( ) ( )mnxmhny m −= ∑ ∞ −∞= （1-2） 我们希望 x(n) 通过这个系统后得到的 y(n) 尽可能接近于 s(n) ，因此，称 y(n) 为 信号 s(n)的估值。按照最小均方误差准则，h(n) 应满足下面的正则方程： ( ) ( ) ( )mkmhk xx m xs −= ∑ ∞ −∞= φφ （1-3） 这就是著名的维纳－霍夫方程，其中 ( )mxxφ 是 x(n)的自相关函数，定义为： ( ) ( ) ( )][ mnxnxEm xx +=φ （1-4） ( )m xs φ 是 x(n)与 s(n)的互相关函数，定义为 ( ) ( ) ( )][ mnsnxEm xs +=φ （1-5） 这里，E[·] 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示求数学期望。 在要求 h(n)满足因果性的条件下，求解维纳-霍夫方程是一个典型的难题。虽然 目前有几种求解 h(n)的解析方法，但它们在计算机上实现起来非常困难。因此， 本实验中，利用近似方法，即最佳 FIR 维纳滤波方法，在计算机上实现随机信 号的维纳滤波。 设 h(n)为一因果序列，其长度为 N，则 ( ) ( ) ( )mnxmhny N m −=∑ − = 1 0 （1-6） 同样利用最小均方误差准则，h(n)满足下面方程： xsxx rhR = （1-7） 其中， ( ) ( ) ( ) T Nhhhh ]1,,1,0[ −= ⋯ （1-8） ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ −− +− +−− = 021 201 110 xxxxxx xxxxxx xxxxxx xx NN N N R φφφ φφφ φφφ ⋯ ⋯⋯⋯⋯ ⋯ ⋯ （1-9） ( ) ( ) ( )[ ]T xsxsxsxs Nr 110 −= φφφ ⋯ （1-10） 这里 T 表示转置运算。称 xxR 为信号 x(n)的 N 阶自相关矩阵， xsr 为 x(n)与 s(n)的互相关函数向量。当 xxR 为满秩矩阵时，由公式(1-7)可得 xsxx rRh 1−= （1-11） 由此可见，利用有限长的 h(n)实现维纳滤波器，只要已知 xxR 和 xsr ,就可以 按上式解得满足因果性的 h。只要 N 选择的足够大，它就可以很好地逼近理想的 维纳滤波器，这一点我们可以在下面实验中得到证实。 在本实验中，s(n)由下式来确定： ( ) ( ) ( ) nwnasns +−= 1 (1-12) 我们把式（1-12）称为信号的生成模型，其中 ,95.0=a , w(n)是零均值方差为 22 1 a w −=σ 的高斯白噪声，v(n)是与 s(n)互不相关的高斯白噪声，其均值为零， 方差 1 2 = v σ 。 根据理论推导，此时维纳最佳滤波器为: ( ) 17239.01 2379.0 −− = z zH (1-13) 单位脉冲响应为: ( ) ( ) ( )nunh n7239.02379.0= (1-14) 由此可以实现对信号 x(n)的最佳过滤，即 ( ) ( ) ( ) ( ) nxnsnsny 2379.01ˆ7239.0ˆ +−== （1-15） 其中 ( ) ( )nns sˆ 为 的最佳估值。同时，可以推出，经过理想维纳滤波后，均方误差 应为: ( )[ ] ( ) ( )( )[ ] 2379.0ˆ 22 =−= nsnsEneE （1-16） 在实验中，我们利用下面公式来统计均方误差： ( ) ( )[ ]∑ = −= L i isis L E 1 2ˆ 1ˆ （1-17） 其中 L 为维纳滤波数据长度。 通过理论推导，我们可以得 s(n)的自相关函数 ( ) m ss m 95.0=φ ,进而得到 x(n)的 自 相 关 函 数 ( ) ( )mm m xx δφ += 95.0 以 及 s(n) 与 x(n) 的 互 相 关 函 数 ( ) ( )mm ssxs φφ = . 实际中，一般很难确切地知道 ( )mxxφ 和 ( )mxsφ ，通常是利用有限个 x(n)和 s(n) 的样本来估计它们： ( ) ( ) ( ) mixix mL m mL i xx + − = ∑ − =1 1ˆ φ （1-18） ( ) ( ) ( )misix mL m mL i xs + − = ∑ − =1 1ˆ φ （1-19） 为了在检验实际中某次产生序列的自相关特性与理论值的近似程度，我们可以采 用下式进行度量： ( ) ( )( ) ( )∑ ∑ −= −= − = K Km xx K Km xxxx xx m mm φ φφ ρ 2ˆ （1-20） ( ) ( )( ) ( )∑ ∑ −= −= − = K Km xs K Km xsxs xs m mm φ φφ ρ 2ˆ （1-21） 该式表示了自相关函数的理论值与某次实现的实际值的相对平方误差。实验中为 了得到与自相关特性理论值相符的观测序列，往往需要多次产生序列，直到两者 的 相 对 平 方 误 差 ρ 足 够 小 。 本 实 验 中 ， 我 们 取 K=50 ， 并 认 为 01.003.0 s ≤≤ xxx ρρ 且 的序列才是满足要求的。 在上面这部分实验中，s(n)是已知的。但是在实际中如果已知 s(n)，维纳滤波也 就没有多少意义了。因此，上述实验纯粹是为了理解维纳滤波原理而 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 的。下 面我们考虑维纳滤波器的一种应用，即利用维纳预测方法来估计 AR 模型的参 数。 假定 s(n)是一个 p 阶 AR 模型，即： ( ) ( ) ( ) ( )nwpnsansans p =−+−+ ⋯11 （1-21） 其 中 w(n)是均值为零方差等于 2 wσ 的高斯白噪声。在已知准确自相关函数 ( )nssφ 的 情况下，由下面 Yule-Walker 方程可以得到信号生成模型参数 ( ) 2 w,1 σ和piai ⋯= ， εσ 2 wss AR = （1-22） 其中 ssR 为 ( ) ( )11 +×+ pp 的自相关矩阵，其意义类似于(1-9)式，只是将 N 换成 p+1, ( ) ( )nnxx ssφφ 换成 而已，A 为 ( ) 11 ×+p 的系数列向量，定义为 [ ]21,,1 paaA ⋯= （1-23） ( ) 11 ×+p为ε 的单位列向量，除第一个元素等于 1 外，其余元素均为零，即 [ ]T0,,0,1 ⋯=ε （1-24） 三、 实验步骤及分析 实验流程图如图 1 所示。其中 ( ) ( )nSˆˆ R、nS I 分别是理想维纳滤波和 FIR 滤波的输 出， 2 R 2 I 2 ˆˆˆ eee x 、、 分别为滤波前，理想滤波后和 FIR 滤波后的噪声方差。 开 始 产 生 L个 v(n),w( n) ,s (n )和 x(n ),利 用 L个 s ( n)和 x (n) ,估 计 R x x 和 r x s 。 输入信号样本个数 L ， F I R 滤 波 器 阶 数 N 利 用 式 （ 1-2 0 ） 检验产生序列 x (n)的自相关和互相关函数是否与理 论 值 相 符 在同一坐标内绘出最后 1 0 0 个 s( n) 和 x (n )。 调矩阵求逆子程序计算 ，将 N个 理 想 的 h( n )和估计的绘于同一坐标内 利 用 公 式 （ 1-1 5 ）进行理想的维纳滤波得 L个 ，将 最 后 1 0 0 个 s (n)和 绘于同一坐标内 ̂ ( )lS n ̂ ( )lS n 在同一坐标下绘出 x (n)的 自 相 关函数的理论值和实际值 由 ， 根 据 公 式 （ 1-6）对 x (n)进行过 滤得 L 个， 将 最 后 1 0 0 个 s (n )和 绘于同一坐标内 。 ̂ ( )h n ̂ ( )h n L个 x( n) ,s (n), 和 。 根 据 公 式 （ 1- 1 7 ） ， 统计并打印出 ̂ ( ) l S n ̂ ( ) R S n 2 2 2 x l R e e e ̂ ̂ ̂ 结 束 Y 图 1 维纳滤波器编制流程 1、运行维纳滤波程序，选择 L=5000,N=10，记录实验结果并分析比较下列 三个问题： ① 与 s(n)比较，信号 x(n)在维纳滤波前后有何差别？滤波效果如何？（注 意：本实验中比较噪声方差时，我们统一取 100 次实现的平均） 答：图 2、图 3 为滤波前后的序列 x(n)（红）和 s(n)（蓝）。可知滤波前后 x(n) 围绕 s(n)的波动比较大，这种变化是由滤波前有很大噪声造成的。滤波后 x(n) 变得比较光滑，与 s(n)更为接近。这说明维纳滤波器对广义平稳输入的滤波效 果是相当明显。从波形上可以看出经过维纳滤波获得的信号比原信号有一定的滞 后，这是维纳滤波器的因果性造成的。 图 2 滤波前的序列 x(n)（红）和 s(n)（蓝） 图 3 滤波后的序列 x(n)（红）和 s(n)（蓝） ② 由图 5 估计出的 hˆ(n)和理想的 h(n)比较，很相似。分析如下：A、L=2000， N=10 时，估计出的 hˆ(n)和理想的 h(n)见图 6-1。B、固定 N=10 不变，L=5000 时，估计出的 hˆ(n)和理想的 h(n)见图 5。C、固定 L=5000，当 N=3，h(n)的估计 值与真值的图 6-2所示。通过对图 5 与图 6-1的信号进行对比可以看出，当滤波 器的阶数一定，观测数据的长度增加时，可以减小输出信号与期望信号之间的差 值。因此，观测数据的长度对实验有着重要的影响，即增加输入信号样本个数可 以提高维纳滤波的性能。对比图 5 与图 6-2中的信号可知，数据的长度一定时， 可以通过改变滤波器的阶数来减小最小误差，从而改善滤波器的滤波效果。但是 通过大量的实验可知，当阶数达到某个值时，误差的改善不再明显，因此滤波器 的阶数对实验结果有很大影响，增加滤波器的阶数可以提高滤波器的性能。通过 以上分析可知，增加输入信号样本个数 L 和滤波器的阶数 N，可以显著地提高维 纳滤波的性能。本实验样本个数取 L=5000，滤波器阶数 N=10,检验得到符合要 求的序列 x(n)（ 01.003.0 s ≤≤ xxx ρρ 且 ）。序列 x(n)自相关函数及和 x(n)与 s(n)的 互相关函数理论值相似程度最接近时分别见图 4a，图 4b。这样，我们可以通过 矩阵求逆运算，得到维纳滤波器的最佳解。从而达到我们设计维纳滤波器的任务， 即选择 h(n)，使其输出信号 y(n)与期望信号 d(n)误差的均方值为最小。 图 4a x(n)自相关函数 图 4b x(n)与 s(n)的互相关函数 图 5 估计值 hˆ(n)（红色）和理想 h(n)—N=5000，L=10 图 6-1 估计值 hˆ(n)（红色）和理想 h(n)—L 固定，N=2000 图 6-2 估计值 hˆ(n)（红色）和理想 h(n)—N 固定，L=3 ③ 理想的维纳滤波和 FIR 维纳滤波效果有何差异？ 答：二者滤波效果分别见图 7、图 8。从仿真 100 次后的均方误差值上看， 0.2494 ˆ0.2506, ˆ 22 == ri ee ，可知 FIR 维纳滤波的均方误差比理想维纳滤波的均方误 差非常接近。由于每次产生的信号都不一样，仿真结果会有一些差异，但 FIR 维纳滤波效果总是逼近与维纳滤波。 图 7 维纳滤波的 s(n)（红）理想的 s(n)(蓝 图 8 FIR 滤波后得到的 s(n)（红）的估计值与理想的 s（n）（蓝） ④ 若去掉流程图 1.1 中判断数据自相关和互相关特性的步骤，可能得出理想 维纳滤波效果不如 FIR 滤波的结论，思考原因。 答：因为如果去掉流程图中判断数据自相关和互相关特性的步骤会使自相关和互 相关中的估计值和真值的误差较大，可能得出理想维纳滤波器效果不如 FIR 滤波 器的结论，其原因是：实际上在用 1 xx xs h R r −= 估计有限长因果序列 h(n)时，由于 ̂ ( ) xx mφ 和 ̂ ( ) xs mφ 是利用有限个 x(n)和 s(n)来估计的，每次仿真 x(n)和 s(n)也都 不同，因此 h(n)会有差异。这样我们利用维纳滤波器的得到了信号会于原来的 输入信号 x(n)误差较大，只是滤波效果很差。 2、维纳预测 开 始 输 入 ：A R 模 型 的 阶 数 p, A R 模 型 的 参 数 ，信 号 s(n)样 本 个 数 L2, 1, , i w a i p σ= ⋯ 根 据 公 式 （1-1 8 ） （ 将x换 为 s),由L个s(n)估 计 p+1个 自 相 关 函 数，解 方 程 （1-2 2 ）得 ̂ ̂1 , , pa a⋯ 结 束 当 p=1,a1=-0.6,L=100， 2 w σ =1，我们所的到的实验结果如图所示： （1）蓝色表示信号 sn，红色表示噪声 wn (2) 由图 2 所示我们记录估计值 A=[1.1,-0.7],我们可以发现估计值 1 ˆ a 与理论值之比有 一定的误差. (3) 由图 3 所示我们记录下估计值 2ˆ w σ 约为 1.1，和理论值有一定的误差。 3、固定 p=1, 6.01 −=a , 1 2 = w σ ,当 L=50 和 500 时我们观测到的信号生成模型的参 数相比较如下图所示： 由上图所示我们可以的出 L越大我们所得到的估计值就越接近理论值，误差就越 小。 四、思考题 解：1、(1) ( ) ( ) ( ) ( ) ( )( ) + −− − + − ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − ⋅ −− − ⋅ − − ⋅= ⎥ ⎦ ⎤ ⎢ ⎣ ⎡Φ = z z zzz z zB z zB zH xs w opt 723947.01 95.01 95.0195.01 95.01 723947.01 95.01 312249.1 1 1 1 2 1 1 12 σ ( )( ) ( ) ( ) ( ) ( )138132.1312249.095.0312249.0 38132.1 312249.0 95.0 312249.0 723947.01 95.01 95.0195.01 95.01 1 1 2 1 −−−= ⎥⎦ ⎤ ⎢⎣ ⎡ − − + − =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − − ⋅ −− − − − − nunu zz Z z z zz Z nn ( ) 1 11 1 723947.01 237949.0 95.01 312249.0 . 723947.01 95.01 312249.1 1 − −− − − = −− − ⋅=∴ z zz z zH opt 这里保留了 6位有效数字，所以公式（1-13）得证。 (2) ( )[ ] ( ) ( ) ( ){ } ( ) ( ) ( )[ ] ( )( ) ( )( ) ( )( ) ( ) 237950.0 05263.1723947.0 782105.0 723947.005263.1 0782105.0 2 1 723947.0195.01 0742999.0 2 1 723947.01 237949.0 1. 95.010.95z-1 0.0975 2 1 v 2 1 1 1 1 11min 2 1 11 min 2 = − − = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ −− − = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ −− = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − − − =∴ Φ=Φ∴ Φ−Φ= ∫ ∫ ∫ ∫ − − − −− − −− dz zzj dzz zzj dzz zzj neE zzs dzzzzHz j neE xsss xsoptss π π π π 不相关，与∵ 综上式（1-16）得证。 2、我们由（1-21）式， ( ) ( ) ( ) ( )nwpnsansans p =−++−+ ⋯11 得到 Yule—Walker 方程为 εσ 2 wss AR = .Yule—Walker 方程式为 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ − +− −− 0 0 1 01 101 10 2 1 ⋮⋮ ⋯ ⋮⋱⋮⋮ ⋯ ⋯ p pp p xxxxxx xxxxxx xxxxxx a a pp p p σ φφφ φφφ φφφ 从一阶信号生成模型开始我们求一阶参数 1a 和 2 p σ 。 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )( )⎩⎨ ⎧ −= −= ⇒−=⇒ ⎩ ⎨ ⎧ += −= ⇒⎥ ⎦ ⎤ ⎢ ⎣ ⎡ =⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ 2 0 2 11 2 1 2 0112 11 2 1 11 2 1 11 2 1 11 1 /1 01 10 0/1 0 1 01 10 σσ σφ φσ φφσ φφ σ φφ φφ a a a a a a xx xx xxxx xxxx xxxx xxxx 在看二阶信号生成模型 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ 0 0 1 012 101 210 22 22 21 σ φφφ φφφ φφφ a a xxxxxx xxxxxx xxxxxx 用克莱姆法则解得出 ( ) ( ) ( ) ( ) ( ) ( ) ( )[ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ −= += − +− = +−= − +− = 2 1 2 22 2 2 1122112221 2 11122 2 22 1 10 2110 /12 10 120 σσ φφ φφφφ σφφ φφ φφφ a aaaa aa xxxx xxxxxxxx xxxx xxxx xxxxxx 依次我们可以类推得到递推公式： ( ) ( ) ( ) ( )0,1 / 2 0 2 1 22 ,1 2 1 1 1 ,1 xxkkkk ikkkikki kxx k l lkxxkk a aaaa lkaka φσσσ σφφ =−= += ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ −+−= − −− − − = −∑ 所以我们就求得到了信号生成模型的系数 2w1, σ和paa ⋯ . 五、实验结论 通过以上实验分析，增加输入信号样本个数 L 和滤波器的阶数 N，可以显 著提高维纳滤波器的性能。维纳滤波可以很好的减小误差，但是实现维纳滤波的 要求是：①输入过程是广义平稳的；②输入过程的统计特性是已知的。然而，由 于输入过程取决于外界的信号、干扰环境，这种环境的统计特性常常是未知的、 变化的，因而难以满足上述两个要求。这就促使人们研究自适应滤波器，也即我 们进一步要学习的内容。 附录（源程序） L=5000; a=0.95; K=50; N=10; sigma_a2=1-a^2; a_=[1,-a]; while(1) wn =sqrt(sigma_a2)*(randn(L,1)); sn=filter(1,a_,wn); vn=randn(L,1); xn=sn+vn; r_xx=xcorr(xn,'unbiased'); r_xx_t=a.^abs([-K:K]); r_xx_t(K+1)=r_xx_t(K+1)+1; p=xcorr(sn,xn,'unbiased'); r_xs=p(L:L+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-a.^[0:K]').^2)/sum(a.^[0:K].^2); if rou_xx<0.03&rou_xs<0.01 break; end end figure(1),clf stem(r_xx(L-K:L+K),'r') hold on stem(r_xx_t,'b') title('r_xx的估计值（红色）及真值（蓝色）') figure(2),clf stem(r_xs,'r') hold on stem(a.^[0:K]) title('r_xs的估计值（红色）及真值（蓝色）') figure(3),clf stem(xn(L-100:L),'r') hold on stem(sn(L-100:L),'b') title('最后 100个 s(n)（蓝色）和 x(n)（红色）'); %%%% 计算估计的 h(n) N=10; n=0:N-1; R_xx=zeros(N,N);%定义 N阶方阵 for i=1:N R_xx(i,:)=r_xx(L+i-1:-1:L-N+i); end %生成自相关矩阵 hopt=inv(R_xx)*p(L:L+N-1); %h(n)的估计值 hopt_t=0.2379*(0.7239).^n; %h(n)的理想值 figure(4),clf stem(hopt,'-','b'); hold on stem(hopt_t,'*','r'); title('h(n)估计值(蓝色-),与真值(红色*)的比较'); %%%% 进行理想的维纳滤波得 L 个 sn_i,将最后 100个 s(n)和 sn_i 绘于同一坐标 sn_i=filter(hopt_t,1,xn); figure(5),clf stem(sn_i(L-100:L),'*','b'); hold on stem(sn(L-100:L),'-','r') title('最后 100个 sn（红色）和由维纳滤波器得到的 sn_i（蓝色）的比较'); sn_r=filter(hopt,1,xn); figure(6),clf stem(sn_r(L-100:L),'*','b'); hold on; stem(sn(L-100:L),'-','r'); title('最后 100个 sn（红色）和对 x(n)进行过滤得到的 sn_r（蓝色）的比较'); %%%%计算 e_x,e_i,e_r e_x=sum((xn(1:L)-sn(1:L)).^2)/L e_i=sum((sn_i(1:L)-sn(1:L)).^2)/L e_r=sum((sn_r(1:L)-sn(1:L)).^2)/L 2、%% 信号生成模型参数估计 L=500; p=1; a1=-0.6; sigma_w2=1; wn=sqrt(sigma_w2)*(randn(L,1)); a=[1,a1]; sn=filter(1,a,wn); r_ss=xcorr(sn,'unbiased'); figure(1),clf stem(wn,'r') hold on stem(sn,'b') %%%% 由 Yule-Walker 方程求得模型参数 ai R_ss=zeros(p+1,p+1); for i=1:p+1 R_ss(i,:)=r_ss(L+i-1:-1:L-(p+1)+i); end I=zeros(1,p+1); for i=1:p+1 if i==1 I(i)=1; else I(i)=0; end end A=inv(R_ss)*I'*sigma_w2; figure(2),clf stem(A,'r') hold on stem(a,'b') sigma_w2_t=r_ss(L)+I(p+1)*r_ss(L+1);%估计值 figure(3),clf stem(sigma_w2_t,'r') hold on stem(sigma_w2,'b') 实验一维纳滤波器的计算机实现