現(xiàn)代數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告

1、現(xiàn)代數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告1、估計(jì)隨機(jī)信號(hào)的樣本自相關(guān)序列。先以白噪聲為例。(a) 產(chǎn)生零均值單位方差高斯白噪聲的1000個(gè)樣點(diǎn)。(b) 用公式： 估計(jì)的前100個(gè)自相關(guān)序列值。與真實(shí)的自相關(guān)序列相比較，討論你的估計(jì)的精確性。(c) 將樣本數(shù)據(jù)分成10段，每段100個(gè)樣點(diǎn)，將所有子段的樣本自相關(guān)的平均值作為自相關(guān)的估值，即： 與(b)的結(jié)果相比，該估計(jì)值有什么變化？它更接近真實(shí)自相關(guān)序列嗎？(d) 再將1000點(diǎn)的白噪聲通過濾波器產(chǎn)生1000點(diǎn)的y(n)，試重復(fù)(b)的工作，估計(jì)y(n)的前100個(gè)自相關(guān)序列值，并與真實(shí)的自相關(guān)序列相比較，討論你的估計(jì)的精確性。仿真結(jié)果：(a)圖1.1 零均值單

2、位方差高斯白噪聲的1000個(gè)樣本點(diǎn)分析圖1.1：這1000個(gè)樣本點(diǎn)是均值近似為0，方差為1的高斯白噪聲。(b) 圖1.2 的前100個(gè)自相關(guān)序列值分析上圖可知：當(dāng)k=0時(shí)取得峰值，且峰值大小比較接近于1，而當(dāng)k0時(shí)估計(jì)的自相關(guān)值在0附近有小幅度的波動(dòng)，這與真實(shí)自相關(guān)序列rx(k)=(k)比較接近，k0時(shí)估計(jì)值非常接近0，說明了估計(jì)的結(jié)果是比較精確的。(c) 圖1.3基于Bartlett法的前100個(gè)自相關(guān)序列值與(b)的結(jié)果相比，同樣在k=0時(shí)達(dá)到峰值，k0時(shí)0值附近上下波動(dòng)；估計(jì)值的方差比較小，隨著k的增大波動(dòng)幅度逐漸變小，在k較大時(shí)它更接近真實(shí)自相關(guān)序列。即采用分段方法得到的自相關(guān)序列的估

3、計(jì)值更加接近rx(k)=(k)。分析仿真圖也可以看出：將樣本數(shù)據(jù)分段，將所有子段的樣本自相關(guān)的平均值作為自相關(guān)的估值時(shí)，可以有效的降低自相關(guān)估計(jì)的方差，而分段樣本估計(jì)的優(yōu)點(diǎn)在于，估計(jì)自相關(guān)序列與實(shí)際自相關(guān)序列的方差減小，且當(dāng)分段數(shù)越大，估計(jì)值越趨向于無偏估計(jì)。(d) 圖1.4 y(n)的前100個(gè)自相關(guān)序列值與真實(shí)值的對(duì)比從圖中可以看出在k=0時(shí)估計(jì)與真實(shí)的自相關(guān)序列之間有較小的誤差，隨著k的增大，估計(jì)得到的值有較大的波動(dòng)，存在一定誤差。源程序clcclear%產(chǎn)生1000個(gè)高斯白噪聲的樣本點(diǎn)x=randn(1,1000);K=1000;figure(1);k=0:K-1; stem(k,x,

4、'.'); %繪制1000個(gè)高斯白噪聲title('零均值單位方差高斯寶噪聲，1000個(gè)樣本點(diǎn)');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('根據(jù)樣本點(diǎn)估計(jì)出的前100自相關(guān)序列值')

5、;xlabel('k');ylabel('r_essk');hold on; realvalue=1,zeros(1,99); stem(k,realvalue,'r','.'); legend('根據(jù)樣本點(diǎn)估計(jì)出的前100自相關(guān)序列值','真實(shí)的自相關(guān)序列'); 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

6、,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法估計(jì)功率譜方法得出的前100個(gè)自相關(guān)序列值');xlabel('k');ylabel('r_ess2k');legend(&

7、#39;Bartlett法估計(jì)功率譜方法得出的前100個(gè)自相關(guān)序列值','真實(shí)的自相關(guān)序列'); 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)/100

8、0;endfigure(4);stem(r_ess3,'.');title('yn前100個(gè)自相關(guān)序列估計(jì)值');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.'

9、;,'.');legend('yn前100個(gè)自相關(guān)序列估計(jì)值','yn的真實(shí)自相關(guān)序列');2、計(jì)算機(jī)練習(xí)2：AR過程的線性建模與功率譜估計(jì)?？紤]AR過程：是單位方差白噪聲。(a) 取b(0)=1, a(1)=0.-7348, a(2)=-1.8820, a(3)=-0.7057, a(4)=-0.8851，產(chǎn)生x(n)的N=256個(gè)樣點(diǎn)。(b) 計(jì)算其自相關(guān)序列的估計(jì)，并與真實(shí)的自相關(guān)序列值相比較。(c) 將的DTFT作為x(n)的功率譜估計(jì)，即：。(d) 利用所估計(jì)的自相關(guān)值和Yule-Walker法(自相關(guān)法)，估計(jì)和的值，并討論估計(jì)的精度

10、。(e) 用(d)中所估計(jì)的和來估計(jì)功率譜為：。(f) 將(c)和(e)的兩種功率譜估計(jì)與實(shí)際的功率譜進(jìn)行比較，畫出它們的重疊波形。(g) 重復(fù)上面的(d)(f)，只是估計(jì)AR參數(shù)分別采用如下方法：(1) 協(xié)方差法；(2) Burg方法；(3) 修正協(xié)方差法。試比較它們的功率譜估計(jì)精度。仿真結(jié)果：（a）圖2.1 x(n)的N=256個(gè)樣點(diǎn)（b）圖2.2自相關(guān)序列的估計(jì)值與真實(shí)的對(duì)比圖2.2中估計(jì)的自相關(guān)序列與真實(shí)的自相關(guān)序列差異較大；在k>100后，真實(shí)的自相關(guān)序列就波動(dòng)得很小，而估計(jì)的自相關(guān)序列則仍有較大的波動(dòng)，但總體上來言兩者都隨著k的增大而逐漸衰減，當(dāng)k>150時(shí)，真實(shí)自相關(guān)

11、序列逐漸趨于0，而估計(jì)出的自相關(guān)序列卻仍有較大的波動(dòng)，這是因?yàn)楣烙?jì)的點(diǎn)數(shù)較少，使得估計(jì)精度不夠。（c）圖2.3 估計(jì)的功率譜與真實(shí)功率譜對(duì)比（d）Yule-Walker法(自相關(guān)法)估計(jì)的參數(shù)值為：b(0)= 1.1537a(1) a(2) a(3) a(4)=- 0.7174 -1.7952 -0.6387 -0.8167圖2.4 Yule-Walker法估計(jì)的功率譜與真實(shí)功率譜對(duì)比Yule-Walker法(自相關(guān)法)估計(jì)的參數(shù)，其系數(shù)的符號(hào)正確，但數(shù)值大小相差較大，并且多次實(shí)驗(yàn)的相差值也很大，參數(shù)估計(jì)的精度遠(yuǎn)遠(yuǎn)不夠。因此從圖2.4中也能看出，該方法得到功率譜與真實(shí)的譜相差很大（e）協(xié)方差法

12、圖2.5 協(xié)方差法估計(jì)的功率譜與真實(shí)功率譜對(duì)比采用協(xié)方差法估計(jì)的參數(shù)，其系數(shù)與真實(shí)的系數(shù)非常接近，該方法能夠較精確的估計(jì)出功率譜。（f）修正協(xié)方差圖2.6 修正的協(xié)方差法估計(jì)的功率譜與真實(shí)功率譜對(duì)比采用修正的協(xié)方差法估計(jì)的參數(shù)，其系數(shù)雖然沒有協(xié)方差法和burg法那么精確，但是估計(jì)誤差也很小。從圖2.6中也能看出，該方法能夠較精確的估計(jì)出功率譜。（g）Burg算法圖2.7 burg法估計(jì)的功率譜與真實(shí)功率譜對(duì)比采用burg估計(jì)的參數(shù)，其系數(shù)幾乎等于真實(shí)的系數(shù)，分析圖2.7中也能看出，該方法非常精確的估計(jì)出功率譜，幾乎與真實(shí)的功率譜曲線重合。源程序：clc;clear;N=256;NN=20000

13、;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)+

14、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)&#

15、39;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

16、-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)*

17、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=d

18、en+(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

19、(P4(i)2;endPP4=PP4+P4;endfigure(1)stem(1:N,x,'.'); title('xn的256個(gè)樣本點(diǎn)');xlabel('n');ylabel('xn');figure(2)plot(0:N-1,r/50); hold on;plot(0:N-1,r1/50,'r');title('xn的256個(gè)樣本點(diǎn)估計(jì)出的前256個(gè)自相關(guān)序列和真實(shí)值');ylabel('Rx(k)');xlabel('k');legend('估計(jì)值&#

20、39;,'真實(shí)值');grid on;axis(0 256 -40 40);hold off; figure(3)plot(w/pi,10*log10(P/50); hold on;title('功率譜估計(jì)');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,1

21、0*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); ho

22、ld on;title('功率譜估計(jì)比較');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('Rx(k)的DTFT','真實(shí)頻譜');grid on;hold off; figure(5)plot(w/pi,10*log10(PP1/50); hold on;title('Yule-Walker法功率譜估計(jì)比較');ylabel('P(dB)');xlabel('w/pi&#

23、39;);plot(w/pi,10*log10(P5),'r');legend('Yule-Walke法','真實(shí)頻譜');grid on;hold off; figure(6)plot(w/pi,10*log10(PP2/50); hold on;title('協(xié)方差法功率譜估計(jì)比較');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('協(xié)方差法','真實(shí)頻譜');g

24、rid on;hold off; figure(7)plot(w/pi,10*log10(PP3/50); hold on;title('修正協(xié)方差法功率譜估計(jì)比較');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('修正協(xié)方差法','真實(shí)頻譜');grid on;hold off; figure(8)plot(w/pi,10*log10(PP4/50); hold on;title('Burg法功率譜估計(jì)

25、比較');ylabel('P(dB)');xlabel('w/pi');plot(w/pi,10*log10(P5),'r');legend('Burg法','真實(shí)譜');grid on;hold off;3、計(jì)算機(jī)練習(xí)3：維納噪聲抑制(例6.6的擴(kuò)展實(shí)驗(yàn))。假定圖6.8中所需的信號(hào)是一個(gè)正弦序列, , 噪聲序列和都是AR(1) 過程，分別由如下的一階差分方程產(chǎn)生：其中是零均值、單位方差的白噪聲，與不相關(guān)。(a) 試用Matlab程序產(chǎn)生x(n)和的500個(gè)樣點(diǎn)，畫出波形圖。(b) 基于x(n)和的500個(gè)

26、樣點(diǎn)，設(shè)計(jì)p階的最優(yōu)FIR維納濾波器，由估計(jì)，進(jìn)而估計(jì)出，其中階數(shù)p分別取為p=3，6，9，12，試計(jì)算各種情況下估計(jì)時(shí)的平均平方誤差(均方誤差的樣本估計(jì)，要敘述估計(jì)方案)，并畫出對(duì)d(n)估計(jì)的結(jié)果。(c) 有時(shí)輔助觀測(cè)數(shù)據(jù)中也會(huì)漏入一些d(n)信號(hào)，即輔助觀測(cè)信號(hào)不僅是，而是 試針對(duì)p=12的情況，分別取幾個(gè)不同的a值(如0.1, 0.5, 1.0)，研究這時(shí)的噪聲抑制性能。(d) 若只有一路觀測(cè)的1000個(gè)樣點(diǎn)，你能想辦法近似完成對(duì)噪聲的有效抑制嗎？試解釋你的方法的基本原理，敘述你的實(shí)現(xiàn)方案。圖6.8 有輔觀測(cè)數(shù)據(jù)的維納噪聲抑制器的原理圖仿真結(jié)果：（a）圖3.1 的波形圖3.2 x(n)

27、的500個(gè)樣點(diǎn)的波形（b）基于和的500個(gè)樣點(diǎn)，可以得到、求解Wiener-Hopf方程可以得到最優(yōu)FIR維納濾波器。均方誤差的樣本估計(jì)可以用計(jì)算得當(dāng)p=3、6、9、12時(shí)，估計(jì)時(shí)的平均平方誤差分別為0.7849、0.2173、0.0747、0.0453。 圖3.3 濾波器階數(shù)p=3時(shí)的估計(jì)值與真實(shí)值對(duì)比圖3 .4濾波器階數(shù)p=6時(shí)的估計(jì)值與真實(shí)值對(duì)比圖3.5 濾波器階數(shù)p=9時(shí)的估計(jì)值與真實(shí)值對(duì)比圖3.6 濾波器階數(shù)p=12時(shí)的估計(jì)值與真實(shí)值對(duì)比分析以上4個(gè)圖：紅色曲線代表真實(shí)值。藍(lán)色曲線代表估計(jì)值。由結(jié)果來看，隨著濾波器階數(shù)的提高，誤差越來越小，對(duì)的估計(jì)也更加精確了，這是因?yàn)殡A數(shù)越高，使用

28、的自相關(guān)序列的值的個(gè)數(shù)就越多，估計(jì)的值也就越準(zhǔn)確了。（c）圖3.7 漏入的參數(shù)a=0.1時(shí)的估計(jì)值與真實(shí)值對(duì)比圖3.8 漏入的參數(shù)a=0.5時(shí)的估計(jì)值與真實(shí)值對(duì)比圖3.9 漏入的參數(shù)a=1.0時(shí)的估計(jì)值與真實(shí)值對(duì)比漏泄的參數(shù)為0.1、0.5、1.0時(shí)，估計(jì)誤差依次為 0.2312、 1.8892、 2.4204，當(dāng)輔助觀測(cè)信號(hào)受到干擾時(shí)，由仿真結(jié)果可以看出，干擾的程度越深，即的值越大，估計(jì)的性能越差。當(dāng)漏泄系數(shù)時(shí)，還可分辨是正弦波形；當(dāng)漏泄系數(shù)時(shí)，波形已經(jīng)基本失真，不能起到分辨效果。（d）原理框圖如下：如圖，采用“自適應(yīng)濾波”的方法近似完成對(duì)噪聲的抑制；假設(shè)和都是實(shí)值的、零均值過程，且互不相關(guān)

29、，另外假設(shè)是窄帶過程，是寬帶過程，對(duì)，可以近似認(rèn)為與自相關(guān)序列為0，此時(shí)：將觀測(cè)信號(hào)進(jìn)行延遲產(chǎn)生“參考信號(hào)”，將這個(gè)參考信號(hào)通過自適應(yīng)濾波器來估計(jì)出。因?yàn)榕c不相關(guān)，則有另外，由于到自適應(yīng)濾波器的輸入是，因此其輸出從而因?yàn)楹褪腔ゲ幌嚓P(guān)，與也是近似不相關(guān)的，所以，從而最后的均方誤差變?yōu)榭梢?，最小化等效于最小化，所以自適應(yīng)濾波器的輸出就是的最小均方估計(jì)，即實(shí)現(xiàn)了噪聲抑制。源程序：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(

30、'輔助噪聲觀測(cè)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('觀測(cè)信號(hào)x(n)和正弦信號(hào)d(n)');ylabel('x(n)');xlabel('n');legend('觀測(cè)信號(hào)x(n)','正弦信號(hào)d(n)&

31、#39;);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

32、(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=fil

33、ter(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維納濾波器階數(shù)p=3輸出結(jié)果');grid on;hold off;figure(4)plot(0:N-1,d2); hold on;plot(0:N-1,d,'r');title('FIR維納濾波器階數(shù)p=6輸出結(jié)果');grid on;hold off;figure(5)plot(0:N-1,d3); hold on;plot(0:N-1,d,'r');title('FIR維納濾波器階數(shù)p=9輸出結(jié)果');grid on;