现代数字信号处理实验报告.doc

1、现代数字信号处理实验报告1、估计随机信号的样本自相关序列。先以白噪声 为例。()xn(a) 产生零均值单位方差高斯白噪声的 1000 个样点。(b) 用公式： 901()()xnrkxk估计 的前 100 个自相关序列值。与真实的自相关序列 相比()xn ()xrk较，讨论你的估计的精确性。(c) 将样本数据分成 10 段，每段 100 个样点，将所有子段的样本自相关的平均值作为 自相关的估值，即：()x901(1)(10) ,1.9mnrkxxnkmk与(b)的结果相比，该估计值有什么变化？它更接近真实自相关序列吗？()xr(d) 再将 1000 点的白噪声 通过滤波器 产生 1000 点的

2、 y(n)，()x1()0.9Hzz试重复(b) 的工作，估计 y(n)的前 100 个自相关序列值，并与真实的自相关序列 相比较，讨论你的估计的精确性。()yrk仿真结果：(a)k0 100 200 300 400 500 600 700 800 900 1000xk-4-3-2-10123 零均值单位方差高斯宝噪声， 1000个样本点图 1.1 零均值单位方差高斯白噪声的 1000 个样本点分析图 1.1：这 1000 个样本点是均值近似为 0，方差为 1 的高斯白噪声。(b) k0 10 20 30 40 50 60 70 80 90 100ressk-0.200.20.40.60.81

3、1.2 根据样本点估计出的前 100自相关序列值根据样本点估计出的前 100自相关序列值真实的自相关序列图 1.2 的前 100 个自相关序列值()xn分析上图可知：当 k=0 时取得峰值，且峰值大小比较接近于 1，而当 k0时估计的自相关值在 0 附近有小幅度的波动，这与真实自相关序列 rx(k)=(k)比较接近，k0 时估计值非常接近 0，说明了估计的结果是比较精确的。(c) k0 10 20 30 40 50 60 70 80 90 100ress2k-0.200.20.40.60.811.2 Bartlett法估计功率谱方法得出的前 100个自相关序列值Bartlett法估计功率谱方法

4、得出的前 100个自相关序列值真实的自相关序列图 1.3 基于 Bartlett 法的前 100 个自相关序列值与(b)的结果相比，同样在 k=0 时达到峰值，k0 时 0 值附近上下波动；估计值的方差比较小，随着 k 的增大波动幅度逐渐变小，在 k 较大时它更接近真实自相关序列 。即采用分段方法得到的自相关序列的估计值更加接()xr近 rx(k)=(k)。分析仿真图也可以看出：将样本数据分段，将所有子段的样本自相关的平均值作为 自相关的估值时，可以有效的降低自相关估计的方差，n而分段样本估计的优点在于，估计自相关序列与实际自相关序列的方差减小，且当分段数越大，估计值越趋向于无偏估计。(d)

5、k0 10 20 30 40 50 60 70 80 90 100ress3(k)0123456 yn前 100个自相关序列估计值yn前 100个自相关序列估计值yn的真实自相关序列图 1.4 y(n)的前 100 个自相关序列值与真实值的对比从图中可以看出在 k=0 时估计与真实的自相关序列之间有较小的误差，随着 k 的增大，估计得到的值有较大的波动，存在一定误差。源程序clcclear%产生 1000 个高斯白噪声的样本点x=randn(1,1000);K=1000;figure(1);k=0:K-1; stem(k,x,.); %绘制 1000 个高斯白噪声title(零均值单位方差高斯

6、宝噪声，1000 个样本点 );xlabel(k);ylabel(xk);mean_x=mean(x) var_x=var(x) %for k=0:99 for n=k+1:1000y_ess(n)=x(n)*x(n-k); endr_ess(k+1)=sum(y_ess)/1000; end figure(2);k=0:99;stem(k,r_ess,.); title(根据样本点估计出的前 100 自相关序列值);xlabel(k);ylabel(r_essk);hold on; realvalue=1,zeros(1,99); stem(k,realvalue,r,.); legend(

7、根据样本点估计出的前 100 自相关序列值,真实的自相关序列); error1=r_ess-realvalue;mean_error_b=mean(error1) var_error_b=var(error1) %for k=0:99 for m=0:9 for n=k+1:100 y_ess2(m+1,n)=x(n+100*m)*x(n-k+100*m); endendr_ess2(k+1)=sum(sum(y_ess2)/1000; endfigure(3);k=0:99;stem(k,r_ess2,b.); hold on; realvalue2=1,zeros(1,99);stem(k

8、,realvalue2,r.,.);title(Bartlett 法估计功率谱方法得出的前 100 个自相关序列值 );xlabel(k);ylabel(r_ess2k);legend(Bartlett 法估计功率谱方法得出的前 100 个自相关序列值,真实的自相关序列); error2=r_ess2-realvalue2;mean_error_c=mean(error2) var_error_c=var(error2) %y=zeros(1,1000);B=1;A=1,-0.9;y=filter(B,A,x); r_ess3=zeros(1,100); for k=0:99for n=(k+

9、1):1000r_ess3(k+1)=r_ess3(k+1)+y(n)*y(n-k); endr_ess3(k+1)=r_ess3(k+1)/1000;endfigure(4);stem(r_ess3,.);title(yn前 100 个自相关序列估计值);xlabel(k),ylabel(r_ess3(k);hold on;p=1,zeros(1,99);h=filter(B,A,p);for i=1:100h1(i)=h(101-i);end rh=conv(h,h1);rh=rh(100:199);realvalue3=conv(p,rh);realvalue3=realvalue3(1

10、:100); stem(realvalue3,r.,.);legend(yn前 100 个自相关序列估计值,yn的真实自相关序列);2、计算机练习 2：AR 过程的线性建模与功率谱估计。考虑 AR 过程：()1()(2)(3)(4)(0)xnaaxnaxnaxnbvn是单位方差白噪声。v(a) 取 b(0)=1, a(1)=0.-7348, a(2)=-1.8820, a(3)=-0.7057, a(4)=-0.8851，产生 x(n)的N=256 个样点。(b) 计算其自相关序列的估计 ，并与真实的自相关序列值相比较。()xrk(c) 将 的 DTFT 作为 x(n)的功率谱估计，即：()x

11、rk。1 21()|()|Nj jkjxkPereXe(d) 利用所估计的自相关值和 Yule-Walker 法( 自相关法)，估计和 的值，并讨论估计的精度。(1), 2(3), 4aa(0)b(e) 用(d)中所估计的 和 来估计功率谱为： 。()ak0b 241|(0)|()jx jkkbPeae(f) 将(c) 和 (e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。(g) 重复上面的(d)(f)，只是估计 AR 参数分别采用如下方法：(1) 协方差法；(2) Burg 方法； (3) 修正协方差法。试比较它们的功率谱估计精度。仿真结果：（a）n0 50 100 150

12、200 250 300xn-15-10-5051015 xn的 256个样本点图 2.1 x(n)的 N=256 个样点（b）k0 50 100 150 200 250Rx(k)-40-30-20-10010203040 xn的 256个样本点估计出的前 256个自相关序列和真实值估计值真实值图 2.2 自相关序列的估计值与真实的对比图 2.2 中估计的自相关序列与真实的自相关序列差异较大；在 k100 后，真实的自相关序列就波动得很小，而估计的自相关序列则仍有较大的波动，但总体上来言两者都随着 k 的增大而逐渐衰减，当 k150 时，真实自相关序列逐渐趋于 0，而估计出的自相关序列却仍有较大

13、的波动，这是因为估计的点数较少，使得估计精度不够。（c）w/pi0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2P(dB)-15-10-505101520253035 功率谱估计比较Rx(k)的 DTFT真实频谱图 2.3 估计的功率谱与真实功率谱对比（d）Yule-Walker 法(自相关法)估计的参数值为：b(0)= 1.1537a(1) a(2) a(3) a(4)=- 0.7174 -1.7952 -0.6387 -0.8167w/pi0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2P(dB)-15-10-50510152025303

14、5 Yule-Walker法功率谱估计比较Yule-Walke法真实频谱图 2.4 Yule-Walker 法估计的功率谱与真实功率谱对比Yule-Walker 法(自相关法)估计的参数，其系数的符号正确，但数值大小相差较大，并且多次实验的相差值也很大，参数估计的精度远远不够。因此从图2.4 中也能看出，该方法得到功率谱与真实的谱相差很大（e）协方差法w/pi0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2P(dB)-15-10-505101520253035 协方差法功率谱估计比较协方差法真实频谱图 2.5 协方差法估计的功率谱与真实功率谱对比采用协方差法估计的参数

15、，其系数与真实的系数非常接近，该方法能够较精确的估计出功率谱。（f）修正协方差w/pi0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2P(dB)-15-10-505101520253035 修正协方差法功率谱估计比较修正协方差法真实频谱图 2.6 修正的协方差法估计的功率谱与真实功率谱对比采用修正的协方差法估计的参数，其系数虽然没有协方差法和 burg 法那么精确，但是估计误差也很小。从图 2.6 中也能看出，该方法能够较精确的估计出功率谱。（g）Burg 算法w/pi0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2P(dB)-15-10-5

16、05101520253035 Burg法功率谱估计比较Burg法真实谱图 2.7 burg 法估计的功率谱与真实功率谱对比采用 burg 估计的参数，其系数几乎等于真实的系数，分析图 2.7 中也能看出，该方法非常精确的估计出功率谱，几乎与真实的功率谱曲线重合。源程序：clc;clear;N=256;NN=20000;v1=normrnd(0,1,50,NN);v=v1(:,1:N);r=zeros(1,N);r1=zeros(1,N);w=0:2*pi/100:2*pi;P=zeros(1,length(w);PP1=zeros(1,length(w);PP2=zeros(1,length(

17、w);PP3=zeros(1,length(w);PP4=zeros(1,length(w);for s=1:50x1=filter(1,1,0.7348,1.8820,0.7057,0.8851,v1(s,:);x=x1(1:N);for k=1:Nrx(k)=0;for n=k:Nrx(k)=rx(k)+x(n)*x(n-(k-1);endrx(k)=rx(k)/(N);endr=r+rx;for k=1:Nrx1(k)=0;for n=k:NNrx1(k)=rx1(k)+x1(n)*x1(n-(k-1);endrx1(k)=rx1(k)/(NN);endr1=r1+rx1;for i=1

18、:length(w)P(i)=P(i)+rx(1);for n=2:NP(i)=P(i)+rx(n)*exp(-j*(n-1)*w(i)+rx(n)*exp(j*(n-1)*w(i);endendA=toeplitz(rx(1:4),rx(1:4);B=-rx(2:5);Ap=AB;b0=rx(1);for i=1:4b0=b0+Ap(i)*rx(i+1);endb0=sqrt(b0);for i=1:length(w)P1(i)=1;for k=1:4P1(i)=P1(i)+Ap(k)*exp(-j*k*w(i);endP1(i)=b02/abs(P1(i)2;endPP1=PP1+P1;A

19、=;for k=1:4c=;for l=1:4rr=0;for n=5:Nrr=rr+x(n-l)*x(n-k);endc=c;rr;endA=A c;endB=;for l=1:4rr=0;for n=5:Nrr=rr+x(n-l)*x(n);endB=B;rr;endB=B*(-1);Ap=AB;for i=1:length(w)P2(i)=1;for k=1:4P2(i)=P2(i)+Ap(k)*exp(-j*k*w(i);endP2(i)=x(1)2/abs(P2(i)2;endPP2=PP2+P2;A=;for k=1:4c=;for l=1:4rr=0;for n=5:Nrr=rr

20、+x(n-l)*x(n-k)+x(n-4+l)*x(n-4+k);endc=c;rr;endA=A c;endB=;for l=1:4rr=0;for n=5:Nrr=rr+x(n-l)*x(n)+x(n-4+l)*x(n-4);endB=B;rr;endB=B*(-1);Ap=AB;for i=1:length(w)P3(i)=1;for k=1:4P3(i)=P3(i)+Ap(k)*exp(-j*k*w(i);endP3(i)=x(1)2/abs(P3(i)2;endPP3=PP3+P3;p=4;ef=zeros(1+p,N);eb=zeros(1+p,N);ef(1,:)=x;eb(1,

21、:)=x;km=;for m=2:p+1mol=0;den=0;for n=m:Nmol=mol+(-2)*ef(m-1,n)*eb(m-1,n-1);den=den+(ef(m-1,n)2+(eb(m-1,n-1)2;endkm(m-1)=mol/den;for n=m:Nef(m,n)=ef(m-1,n)+km(m-1)*eb(m-1,n-1);eb(m,n)=eb(m-1,n-1)+km(m-1)*ef(m-1,n);endendaa=1;for i=1:4aa=aa;0;bb=aa(length(aa):-1:1);aa=aa+bb*km(i);endfor i=1:length(w)

22、P4(i)=1;for k=2:5P4(i)=P4(i)+aa(k)*exp(-j*(k-1)*w(i);endP4(i)=1/abs(P4(i)2;endPP4=PP4+P4;endfigure(1)stem(1:N,x,.); title(xn的 256 个样本点);xlabel(n);ylabel(xn);figure(2)plot(0:N-1,r/50); hold on;plot(0:N-1,r1/50,r);title(xn的 256 个样本点估计出的前 256 个自相关序列和真实值);ylabel(Rx(k);xlabel(k);legend(估计值,真实值);grid on;a

23、xis(0 256 -40 40);hold off;figure(3)plot(w/pi,10*log10(P/50); hold on;title(功率谱估计);ylabel(P(dB);xlabel(w/pi);plot(w/pi,10*log10(PP1/50),r);plot(w/pi,10*log10(PP2/50),g);plot(w/pi,10*log10(PP3/50),y);plot(w/pi,10*log10(PP4/50),k);aap=0.7348,1.8820,0.7057,0.8851;for i=1:length(w)P5(i)=1;for k=1:4P5(i)

24、=P5(i)+aap(k)*exp(-j*k*w(i);endP5(i)=1/abs(P5(i)2;endplot(w/pi,10*log10(P5),:);legend(Rx(k)的 DTFT,Yule-Walker);grid on;hold off;figure(4)plot(w/pi,10*log10(P/50); hold on;title(功率谱估计比较);ylabel(P(dB);xlabel(w/pi);plot(w/pi,10*log10(P5),r);legend(Rx(k)的 DTFT,真实频谱);grid on;hold off;figure(5)plot(w/pi,1

25、0*log10(PP1/50); hold on;title(Yule-Walker 法功率谱估计比较);ylabel(P(dB);xlabel(w/pi);plot(w/pi,10*log10(P5),r);legend(Yule-Walke 法,真实频谱);grid on;hold off;figure(6)plot(w/pi,10*log10(PP2/50); hold on;title(协方差法功率谱估计比较);ylabel(P(dB);xlabel(w/pi);plot(w/pi,10*log10(P5),r);legend(协方差法,真实频谱);grid on;hold off;f

26、igure(7)plot(w/pi,10*log10(PP3/50); hold on;title(修正协方差法功率谱估计比较);ylabel(P(dB);xlabel(w/pi);plot(w/pi,10*log10(P5),r);legend(修正协方差法,真实频谱);grid on;hold off;figure(8)plot(w/pi,10*log10(PP4/50); hold on;title(Burg 法功率谱估计比较);ylabel(P(dB);xlabel(w/pi);plot(w/pi,10*log10(P5),r);legend(Burg 法,真实谱);grid on;h

27、old off;3、计算机练习 3：维纳噪声抑制(例 6.6 的扩展实验) 。假定图 6.8 中所需的信号 是一()dn个正弦序列 , , 噪声序列 和 都是 AR(1) 过程，0)sind0.21()vn2分别由如下的一阶差分方程产生： 11()8()vng22.6(其中 是零均值、单位方差的白噪声，与 不相关。()gn()dn(a) 试用 Matlab 程序产生 x(n)和 的 500 个样点，画出波形图。2v(b) 基于 x(n)和 的 500 个样点，设计 p 阶的最优 FIR 维纳滤波器，由 估计 ，2v 2()vn1()进而估计出 ，其中阶数 p 分别取为 p=3，6，9，12，试

28、计算各种情况下估计d时的平均平方误差(均方误差的样本估计，要叙述估计方案)，并画出对 d(n)估计1的结果。(c) 有时辅助观测数据中也会漏入一些 d(n)信号，即辅助观测信号不仅是 ，而是2v02()v试针对 p=12 的情况，分别取几个不同的 值( 如 0.1, 0.5, 1.0)，研究这时的噪声抑制性能。(d) 若只有一路观测 的 1000 个样点，你能想办法近似完成对噪声1()xndn的有效抑制吗？试解释你的方法的基本原理，叙述你的实现方案。1()vn图 6.8 有辅观测数据的维纳噪声抑制器的原理图仿真结果：（a）n0 50 100 150 200 250 300 350 400 45

29、0 500v2(n)-4-3-2-1012345 辅助噪声观测 v2(n)图 3.1 的波形2()vn0 50 100 150 200 250 300 350 400 450 500x(n)-6-5-4-3-2-101234 观测信号 x(n)和正弦信号 d(n)观测信号 x(n)正弦信号 d(n)图 3.2 x(n)的 500 个样点的波形（b）基于 和 的 500 个样点，可以得到()Xn2()V、120()Nvnrkk1220()()NxvnrxVk求解 Wiener-Hopf 方程可以得到最优 FIR 维纳滤波器。22vxvRw均方误差的样本估计可以用 计算得当 p=3、6、9、12

30、时，估1210()(Nnv计 时的平均平方误差分别为 0.7849、0.2173、0.0747、0.0453。 1()vn0 50 100 150 200 250 300 350 400 450 500-4-3-2-10123 FIR维纳滤波器阶数 p=3输出结果图 3.3 滤波器阶数 p=3 时的估计值与真实值对比0 50 100 150 200 250 300 350 400 450 500-2-1.5-1-0.500.511.52 FIR维纳滤波器阶数 p=6输出结果图 3 .4 滤波器阶数 p=6 时的估计值与真实值对比0 50 100 150 200 250 300 350 400

31、450 500-1.5-1-0.500.511.5 FIR维纳滤波器阶数 p=9输出结果图 3.5 滤波器阶数 p=9 时的估计值与真实值对比0 50 100 150 200 250 300 350 400 450 500-1.5-1-0.500.511.5 FIR维纳滤波器阶数 p=12输出结果图 3.6 滤波器阶数 p=12 时的估计值与真实值对比分析以上 4 个图：红色曲线代表真实值。蓝色曲线代表估计值。由结果来看，随着滤波器阶数的提高，误差越来越小，对 的估计也更加精确了，这是因)(nd为阶数越高，使用的自相关序列的值的个数就越多，估计的值也就越准确了。（c）0 50 100 150

32、200 250 300 350 400 450 500-4-3-2-101234 辅助观测数据中漏入 dn的 0.1估计结果期望值图 3.7 漏入的参数 a=0.1 时的估计值与真实值对比0 50 100 150 200 250 300 350 400 450 500-4-3-2-101234 辅助观测数据中漏入 dn的 0.5估计结果期望值图 3.8 漏入的参数 a=0.5 时的估计值与真实值对比0 50 100 150 200 250 300 350 400 450 500-4-3-2-101234 辅助观测数据中漏入 dn的 1.0估计结果期望值图 3.9 漏入的参数 a=1.0 时的估

33、计值与真实值对比漏泄的参数为 0.1、0.5、1.0 时，估计误差依次为 0.2312、 1.8892、 2.4204，当辅助观测信号 受到 干扰时，由仿真结果可以看出，干扰的)(2nv)(d程度越深，即 的值越大，估计的性能越差。当漏泄系数 时，还可分辨 5.0是正弦波形；当漏泄系数 时，波形已经基本失真，不能起到分辨效果。5.0（d）原理框图如下：如图，采用“自适应滤波”的方法近似完成对噪声 的抑制；假设)(1nv和 都是实值的、零均值过程，且互不相关，另外假设 是窄带过)(nd1v d程， 是宽带过程，对 ，可以近似认为 与 自相关序列为0k)(1v)1k0，此时：)()( )()(2|

34、21 112nydEnv nydvEve将观测信号进行延迟产生“参考信号” ，将这个参考信号通过自适应滤波器来估计出 。因为 与 不相关，则有1)(2)()(211 nyvEnydvE另外，由于到自适应滤波器的输入是 ，因此其输出0kx pk pknn knvkdwxwy00 0100 )()()()(从而 pkn kvEvEyvE0 01011 )()()()(因为 和 是互不相关，与 也是近似不相关的，所以)(nd1 01k，从而最后的均方误差变为1yv )()(|)(| 2212 nydEnveE可见，最小化 等效于最小化 ，所以自适应滤波器|n的输出 就是 的最小均方估计，即实现了噪声

35、抑制。)(ny)(d源程序：clc;clear;N=500;g=normrnd(0,1,1,N); v1=filter(1,1,-0.8,g); v2=filter(1,1,0.6,g); figure(1)plot(0:N-1,v2); title(辅助噪声观测 v2(n);ylabel(v2(n);xlabel(n);grid on;d=sin(0:N-1*0.02*pi-pi/2); x=d+v1; figure(2)plot(0:N-1,x,r); hold on;plot(0:N-1,d,-.);title(观测信号 x(n)和正弦信号 d(n);ylabel(x(n);xlabel

36、(n);legend(观测信号 x(n),正弦信号 d(n);grid on;hold off;%d(n),p=3,6,9,12for k=1:13rv2(k)=0;for n=k:Nrv2(k)=rv2(k)+v2(n)*v2(n-(k-1);endrv2(k)=rv2(k)/N;end%Rxv2(k)for k=1:13rxv2(k)=0;for n=k:Nrxv2(k)=rxv2(k)+x(n)*v2(n-(k-1);endrxv2(k)=rxv2(k)/N;end%p=zeros(1,4);p(1)=3;p(2)=6;p(3)=9;p(4)=12;A=toeplitz(rv2(1:p(

37、1)+1),rv2(1:p(1)+1);B=rxv2(1:p(1)+1);w1=AB;A2=toeplitz(rv2(1:p(2)+1),rv2(1:p(2)+1);B2=rxv2(1:p(2)+1);w2=A2B2;A3=toeplitz(rv2(1:p(3)+1),rv2(1:p(3)+1);B3=rxv2(1:p(3)+1);w3=A3B3;A4=toeplitz(rv2(1:p(4)+1),rv2(1:p(4)+1);B4=rxv2(1:p(4)+1);w4=A4B4;%d(n)v11=filter(w1,1,v2);v12=filter(w2,1,v2);v13=filter(w3,

38、1,v2);v14=filter(w4,1,v2);d1=x-v11; d2=x-v12; d3=x-v13; d4=x-v14; (v1-v11)*(v1-v11)/N(v1-v12)*(v1-v12)/N(v1-v13)*(v1-v13)/N(v1-v14)*(v1-v14)/Nfigure(3)plot(0:N-1,d1); hold on;plot(0:N-1,d,r);title(FIR 维纳滤波器阶数 p=3 输出结果);grid on;hold off;figure(4)plot(0:N-1,d2); hold on;plot(0:N-1,d,r);title(FIR 维纳滤波器

39、阶数 p=6 输出结果);grid on;hold off;figure(5)plot(0:N-1,d3); hold on;plot(0:N-1,d,r);title(FIR 维纳滤波器阶数 p=9 输出结果);grid on;hold off;figure(6)plot(0:N-1,d4); hold on;plot(0:N-1,d,r);title(FIR 维纳滤波器阶数 p=12 输出结果);grid on;hold off;w1w2w3w4N=500;g=normrnd(0,1,1,N); v1=filter(1,1,-0.8,g); v2=filter(1,1,0.6,g); d=

40、sin(0:N-1*0.05*pi-pi/2); x=d+v1; d1=;a=0.1 0.5 1.0;for l=1:length(a)v0=v2+a(l)*d;for k=1:13rv2(k)=0;for n=k:Nrv2(k)=rv2(k)+v0(n)*v0(n-(k-1);endrv2(k)=rv2(k)/N;end%for k=1:13rxv2(k)=0;for n=k:Nrxv2(k)=rxv2(k)+x(n)*v0(n-(k-1);endrxv2(k)=rxv2(k)/N;end%pp=12;A=toeplitz(rv2(1:pp+1),rv2(1:pp+1);B=rxv2(1:p

41、p+1);w1=AB;%v11=filter(w1,1,v0);d1=d1 x-v11; (v1-v11)*(v1-v11)/Nendfigure(7)plot(0:N-1,d1(1:N); hold on;plot(0:N-1,d,r);title(辅助观测数据中漏入 dn的 0.1);legend(估计结果,期望值);axis(0 N -4 4);grid on;hold off;figure(8)plot(0:N-1,d1(N+1:2*N); hold on;plot(0:N-1,d,r);title(辅助观测数据中漏入 dn的 0.5);legend(估计结果,期望值);axis(0 N -4 4);grid on;hold off;figure(9)plot(0:N-1,d1(2*N+1:3*N); hold on;plot(0:N-1,d,r);title(辅助观测数据中漏入 dn的 1.0);legend(估计结果,期望值);axis(0 N -4 4);grid on;hold off;