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

现代数字信号处理实验报告1、估计随机信号的样本自相关序列。先以白噪声为例。(a) 产生零均值单位方差高斯白噪声的1000个样点。(b) 用公式： 估计的前100个自相关序列值。与真实的自相关序列相比较，讨论你的估计的精确性。(c) 将样本数据分成10段，每段100个样点，将所有子段的样本自相关的平均值作为自相关的估值，即： 与(b)的结果相比，该估计值有什么变化？它更接近真实自相关序列吗？(d) 再将1000点的白噪声通过滤波器产生1000点的y(n)，试重复(b)的工作，估计y(n)的前100个自相关序列值，并与真实的自相关序列相比较，讨论你的估计的精确性。仿真结果：(a)图1.1 零均值单位方差高斯白噪声的1000个样本点分析图1.1：这1000个样本点是均值近似为0，方差为1的高斯白噪声。(b) 图1.2 的前100个自相关序列值分析上图可知：当k=0时取得峰值，且峰值大小比较接近于1，而当k0时估计的自相关值在0附近有小幅度的波动，这与真实自相关序列rx(k)=(k)比较接近，k0时估计值非常接近0，说明了估计的结果是比较精确的。(c) 图1.3基于Bartlett法的前100个自相关序列值与(b)的结果相比，同样在k=0时达到峰值，k0时0值附近上下波动；估计值的方差比较小，随着k的增大波动幅度逐渐变小，在k较大时它更接近真实自相关序列。即采用分段方法得到的自相关序列的估计值更加接近rx(k)=(k)。分析仿真图也可以看出：将样本数据分段，将所有子段的样本自相关的平均值作为自相关的估值时，可以有效的降低自相关估计的方差，而分段样本估计的优点在于，估计自相关序列与实际自相关序列的方差减小，且当分段数越大，估计值越趋向于无偏估计。(d) 图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(零均值单位方差高斯宝噪声，1000个样本点);xlabel(k);ylabel(xk);mean_x=mean(x) var_x=var(x) %for k=0:99 for n=k+1:1000 y_ess(n)=x(n)*x(n-k); end r_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(根据样本点估计出的前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); end end r_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,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:99 for n=(k+1):1000 r_ess3(k+1)=r_ess3(k+1)+y(n)*y(n-k); end r_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:100 h1(i)=h(101-i);end rh=conv(h,h1);rh=rh(100:199);realvalue3=conv(p,rh);realvalue3=realvalue3(1:100); stem(realvalue3,r.,.);legend(yn前100个自相关序列估计值,yn的真实自相关序列);2、计算机练习2：AR过程的线性建模与功率谱估计。考虑AR过程：是单位方差白噪声。(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) 计算其自相关序列的估计，并与真实的自相关序列值相比较。(c) 将的DTFT作为x(n)的功率谱估计，即：。(d) 利用所估计的自相关值和Yule-Walker法(自相关法)，估计和的值，并讨论估计的精度。(e) 用(d)中所估计的和来估计功率谱为：。(f) 将(c)和(e)的两种功率谱估计与实际的功率谱进行比较，画出它们的重叠波形。(g) 重复上面的(d)(f)，只是估计AR参数分别采用如下方法：(1) 协方差法；(2) Burg方法；(3) 修正协方差法。试比较它们的功率谱估计精度。仿真结果：（a）图2.1 x(n)的N=256个样点（b）图2.2自相关序列的估计值与真实的对比图2.2中估计的自相关序列与真实的自相关序列差异较大；在k100后，真实的自相关序列就波动得很小，而估计的自相关序列则仍有较大的波动，但总体上来言两者都随着k的增大而逐渐衰减，当k150时，真实自相关序列逐渐趋于0，而估计出的自相关序列却仍有较大的波动，这是因为估计的点数较少，使得估计精度不够。（c）图2.3 估计的功率谱与真实功率谱对比（d）Yule-Walker法(自相关法)估计的参数值为：b(0)= 1.1537a(1) a(2) a(3) a(4)=- 0.7174 -1.7952 -0.6387 -0.8167图2.4 Yule-Walker法估计的功率谱与真实功率谱对比Yule-Walker法(自相关法)估计的参数，其系数的符号正确，但数值大小相差较大，并且多次实验的相差值也很大，参数估计的精度远远不够。因此从图2.4中也能看出，该方法得到功率谱与真实的谱相差很大（e）协方差法图2.5 协方差法估计的功率谱与真实功率谱对比采用协方差法估计的参数，其系数与真实的系数非常接近，该方法能够较精确的估计出功率谱。（f）修正协方差图2.6 修正的协方差法估计的功率谱与真实功率谱对比采用修正的协方差法估计的参数，其系数虽然没有协方差法和burg法那么精确，但是估计误差也很小。从图2.6中也能看出，该方法能够较精确的估计出功率谱。（g）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(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:N rx(k)=0; for n=k:N rx(k)=rx(k)+x(n)*x(n-(k-1); end rx(k)=rx(k)/(N);endr=r+rx;for k=1:N rx1(k)=0; for n=k:NN rx1(k)=rx1(k)+x1(n)*x1(n-(k-1); end rx1(k)=rx1(k)/(NN);endr1=r1+rx1; for i=1:length(w) P(i)=P(i)+rx(1); for n=2:N P(i)=P(i)+rx(n)*exp(-j*(n-1)*w(i)+rx(n)*exp(j*(n-1)*w(i); endend A=toeplitz(rx(1:4),rx(1:4);B=-rx(2:5);Ap=AB;b0=rx(1);for i=1:4 b0=b0+Ap(i)*rx(i+1);endb0=sqrt(b0);for i=1:length(w) P1(i)=1; for k=1:4 P1(i)=P1(i)+Ap(k)*exp(-j*k*w(i); end P1(i)=b02/abs(P1(i)2;endPP1=PP1+P1; A=;for k=1:4 c=; for l=1:4 rr=0; for n=5:N rr=rr+x(n-l)*x(n-k); end c=c;rr; end A=A c;endB=;for l=1:4 rr=0; for n=5:N rr=rr+x(n-l)*x(n); end B=B;rr;endB=B*(-1);Ap=AB;for i=1:length(w) P2(i)=1; for k=1:4 P2(i)=P2(i)+Ap(k)*exp(-j*k*w(i); end P2(i)=x(1)2/abs(P2(i)2;endPP2=PP2+P2; A=;for k=1:4 c=; for l=1:4 rr=0; for n=5:N rr=rr+x(n-l)*x(n-k)+x(n-4+l)*x(n-4+k); end c=c;rr; end A=A c;endB=;for l=1:4 rr=0; for n=5:N rr=rr+x(n-l)*x(n)+x(n-4+l)*x(n-4); end B=B;rr;endB=B*(-1);Ap=AB;for i=1:length(w) P3(i)=1; for k=1:4 P3(i)=P3(i)+Ap(k)*exp(-j*k*w(i); end P3(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,:)=x;km=;for m=2:p+1 mol=0; den=0; for n=m:N mol=mol+(-2)*ef(m-1,n)*eb(m-1,n-1); den=den+(ef(m-1,n)2+(eb(m-1,n-1)2; end km(m-1)=mol/den; for n=m:N ef(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:4 aa=aa;0; bb=aa(length(aa):-1:1); aa=aa+bb*km(i);endfor i=1:length(w) P4(i)=1; for k=2:5 P4(i)=P4(i)+aa(k)*exp(-j*(k-1)*w(i); end P4(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;axis(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:4 P5(i)=P5(i)+aap(k)*exp(-j*k*w(i); end P5(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,10*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; figure(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;hold off;3、计算机练习3：维纳噪声抑制(例6.6的扩展实验)。假定图6.8中所需的信号是一个正弦序列, , 噪声序列和都是AR(1) 过程，分别由如下的一阶差分方程产生：其中是零均值、单位方差的白噪声，与不相关。(a) 试用Matlab程序产生x(n)和的500个样点，画出波形图。(b) 基于x(n)和的500个样点，设计p阶的最优FIR维纳滤波器，由估计，进而估计出，其中阶数p分别取为p=3，6，9，12，试计算各种情况下估计时的平均平方误差(均方误差的样本估计，要叙述估计方案)，并画出对d(n)估计的结果。(c) 有时辅助观测数据中也会漏入一些d(n)信号，即辅助观测信号不仅是，而是 试针对p=12的情况，分别取几个不同的a值(如0.1, 0.5, 1.0)，研究这时的噪声抑制性能。(d) 若只有一路观测的1000个样点，你能想办法近似完成对噪声的有效抑制吗？试解释你的方法的基本原理，叙述你的实现方案。图6.8 有辅观测数据的维纳噪声抑制器的原理图仿真结果：（a）图3.1 的波形图3.2 x(n)的500个样点的波形（b）基于和的500个样点，可以得到、求解Wiener-Hopf方程可以得到最优FIR维纳滤波器。均方误差的样本估计可以用计算得当p=3、6、9、12时，估计时的平均平方误差分别为0.7849、0.2173、0.0747、0.0453。 图3.3 滤波器阶数p=3时的估计值与真实值对比图3 .4滤波器阶数p=6时的估计值与真实值对比图3.5 滤波器阶数p=9时的估计值与真实值对比图3.6 滤波器阶数p=12时的估计值与真实值对比分析以上4个图：红色曲线代表真实值。蓝色曲线代表估计值。由结果来看，随着滤波器阶数的提高，误差越来越小，对的估计也更加精确了，这是因为阶数越高，使用的自相关序列的值的个数就越多，估计的值也就越准确了。（c）图3.7 漏入的参数a=0.1时的估计值与真实值对比图3.8 漏入的参数a=0.5时的估计值与真实值对比图3.9 漏入的参数a=1.0时的估计值与真实值对比漏泄的参数为0.1、0.5、1.0时，估计误差依次为 0.2312、 1.8892、 2.4204，当辅助观测信号受到干扰时，由仿真结果可以看出，干扰的程度越深，即的值越大，估计的性能越差。当漏泄系数时，还可分辨是正弦波形；当漏泄系数时，波形已经基本失真，不能起到分辨效果。（d）原理框图如下：如图，采用“自适应滤波”的方法近似完成对噪声的抑制；假设和都是实值的、零均值过程，且互不相关，另外假设是窄带过程，是宽带过程，对，可以近似认为与自相关序列为0，此时：将观测信号进行延迟产生“参考信号”，将这个参考信号通过自适应滤波器来估计出。因为与不相关，则有另外，由于到自适应滤波器的输入是，因此其输出从而因为和是互不相关，与也是近似不相关的，所以，从而最后的均方误差变为可见，最小化等效于最小化，所以自适应滤波器的输出就是的最小均方估计，即实现了噪声抑制。源程序：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(n);legend(观测信号x(n),正弦信号d(n);grid on;hold off; %d(n),p=3,6,9,12for k=1:13 rv2(k)=0; for n=k:N rv2(k)=rv2(k)+v2(n)*v2(n-(k-1); end rv2(k)=rv2(k)/N;end%Rxv2(k)for k=1:13 rxv2(k)=0; for n=k:N rxv2(k)=rxv2(k)+x(n)*v2(n-(k-1); end rxv2(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(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,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)/N figure(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维纳滤波器阶数p=6输出结果);grid on;hold off;figure(5)plot(0:N-1,d3);