维纳滤波器噪声消除

1 %Problem 1 Since rv1,v2 = E[v1(n)v2 ∗(n − k)] rx,v2 = E[x(n)v2 ∗(n − k)] x(n)=d(n)+v1(n) Given in the problem that d(n) is uncorrelated with v2(n) v1(n)v2 ∗(n − k) = x(n)v2 ∗(n − k) So rv1,v2 = rx,v2 2 3 4 Ra = Columns 1 through 3 1.5413 -0.9164 0.5289 Columns 4 through 6 -0.2953 0.1487 -0.0962 Columns 7 through 9 0.0683 -0.0794 0.0552 Columns 10 through 12 -0.0241 0.0116 0.0082 Columns 13 through 15 -0.0077 -0.0082 -0.0009 Columns 16 through 18 -0.0124 -0.0001 0.0106 Columns 19 through 20 -0.0321 0.0332 5 Rc = Columns 1 through 3 1.5461 -0.9186 0.5308 Columns 4 through 6 -0.2974 0.1502 -0.0972 Columns 7 through 9 0.0695 -0.0808 0.0550 Columns 10 through 12 -0.0247 0.0120 0.0074 Columns 13 through 15 -0.0071 -0.0086 -0.0005 Columns 16 through 18 -0.0124 0.0002 0.0102 Columns 19 through 20 -0.0316 0.0332 6 M=1 Autocorrelation values R1=1.5651 Cross-correlation values P1=0.6598 Filter coefficients -1.7170 7 8 M=6 Autocorrelation mR7 = 1.5736 -0.9400 0.5635 -0.3434 0.2239 -0.1371 -0.9400 1.5736 -0.9400 0.5635 -0.3434 0.2239 0.5635 -0.9400 1.5736 -0.9400 0.5635 -0.3434 -0.3434 0.5635 -0.9400 1.5736 -0.9400 0.5635 0.2239 -0.3434 0.5635 -0.9400 1.5736 -0.9400 -0.1371 0.2239 -0.3434 0.5635 -0.9400 1.5736 Cross-correlation P7 = 0.6876 0.5541 0.4425 0.3489 0.2966 0.2348 Filter coefficients w6 = 0.7248 -0.0430 0.3213 0.0627 0.1309 0.1652 9 10 M=12 Autocorrelation R8 =1.4806 -0.8504 0.4904 -0.2770 0.1580 -0.0684 0.0120 0.0001 -0.0175 0.0178 - 0.0165 0.0134 mR8 = Columns 1 through 10 1.4806 -0.8504 0.4904 -0.2770 0.1580 -0.0684 0.0120 0.0001 -0.0175 0.0178 -0.8504 1.4806 -0.8504 0.4904 -0.2770 0.1580 -0.0684 0.0120 0.0001 -0.0175 0.4904 -0.8504 1.4806 -0.8504 0.4904 -0.2770 0.1580 -0.0684 0.0120 0.0001 -0.2770 0.4904 -0.8504 1.4806 -0.8504 0.4904 -0.2770 0.1580 -0.0684 0.0120 0.1580 -0.2770 0.4904 -0.8504 1.4806 -0.8504 0.4904 -0.2770 0.1580 -0.0684 -0.0684 0.1580 -0.2770 0.4904 -0.8504 1.4806 -0.8504 0.4904 -0.2770 0.1580 0.0120 -0.0684 0.1580 -0.2770 0.4904 -0.8504 1.4806 -0.8504 0.4904 -0.2770 0.0001 0.0120 -0.0684 0.1580 -0.2770 0.4904 -0.8504 1.4806 -0.8504 0.4904 -0.0175 0.0001 0.0120 -0.0684 0.1580 -0.2770 0.4904 -0.8504 1.4806 -0.8504 0.0178 -0.0175 0.0001 0.0120 -0.0684 0.1580 -0.2770 0.4904 -0.8504 1.4806 -0.0165 0.0178 -0.0175 0.0001 0.0120 -0.0684 0.1580 -0.2770 0.4904 -0.8504 0.0134 -0.0165 0.0178 -0.0175 0.0001 0.0120 -0.0684 0.1580 -0.2770 0.4904 Columns 11 through 12 -0.0165 0.0134 0.0178 -0.0165 -0.0175 0.0178 0.0001 -0.0175 0.0120 0.0001 -0.0684 0.0120 0.1580 -0.0684 -0.2770 0.1580 0.4904 -0.2770 -0.8504 0.4904 1.4806 -0.8504 -0.8504 1.4806 Cross-correlation P8 = 0.6725 0.5767 0.4425 0.3717 0.2900 0.2586 0.1784 0.1495 0.1022 0.0889 0.0646 0.0549 11 Filter coefficients w12 = 0.6490 0.0369 0.2679 0.1182 0.1137 0.1448 0.0268 0.0992 -0.0051 0.0682 -0.0125 0.0572 12 There must be something wrong for all three plots from problem 6 to problem 8 don’t seem to get better as expected, but as the order grow larger, the estimate of v1hat should be closer to the real value, so that e(n) will get closer to d(n) 13 % problem #2 %d(n)=sin(n*w0+fi) w0=0.05*pi; fi=2*pi*rand(1)-pi; d(1)=sin(w0+fi); v(1)=randn(1); v1(1)=v(1); v2(1)=v(1); x(1)=d(1)+v1(1); for n=2:10000 d(n)=sin(n*w0+fi); v(n)=randn(1); v1(n)=0.8*v1(n-1)+v(n); v2(n)=-0.6*v2(n-1)+v(n); x(n)=d(n)+v1(n); end %problem no.3 figure(1) dp=plot (d(1:200)); set(dp,'color','red'); xlabel('n'); title('d(n) and x(n)'); hold all xp=plot (x(1:200)); set(xp,'color','blue'); figure(3) v2p=plot (v2(1:200)); xlabel('n'); title('v2(n)'); %problem no.4 v2(n)autocorrelation for ka=0:19 suma=0; 14 for ma=1:10000-ka suma=suma+v2(ma)*v2(ma+ka); end Ra(ka+1)=suma; end Ra=Ra/10000; %problem no.5 x(n) and v2(n) cross-correlation for kc=0:19 sumc=0; for mc=1:10000-kc sumc=sumc+v2(mc)*x(mc+kc); end Rc(kc+1)=sumc; end Rc=Rc/10000; %problem no.6 M=1 sum6=0; for n=1:10000 sum6=sum6+v2(n)*v2(n); end R1=sum6; R1=R1/10000 sum62=0; for n=1:10000 sum62=sum62+v2(n)*x(n); end P1=sum62; P1=P1/10000; w1=R1*P1 ;%inverse R1=R1 w1=v2(1); v1hat(1)=w1*v2(1); for n=1:10000 sum63=0; sum63=sum63+w1*v2(n); v1hat(n)=sum63; end e1=x-v1hat; figure(1); d=plot (d(1:200)); set(d,'color','red'); hold all e1=plot(e1(1:200)); set(e1,'color','blue') 15 title('d(n) and e(n) for first 200 samples') xlabel('n'); legend('d','e1') %problem no.7 M=6 for k7=0:5 sum7=0; for n=1:10000-k7 sum7=sum7+v2(n)*v2(n+k7); end R7(k7+1)=sum7; end R7=R7/10000; mR7=[R7(1) R7(2) R7(3) R7(4) R7(5) R7(6), R7(2) R7(1) R7(2) R7(3) R7(4) R7(5), R7(3) R7(2) R7(1) R7(2) R7(3) R7(4), R7(4) R7(3) R7(2) R7(1) R7(2) R7(3), R7(5) R7(4) R7(3) R7(2) R7(1) R7(2), R7(6) R7(5) R7(4) R7(3) R7(2) R7(1)]; %since v2(n) only has real part, inverse matrix R7 is equal to R7 for k72=0:5 sum7=0; for n=1:10000-k72 sum7=sum7+v2(n)*x(n+k72); end P7(k72+1)=sum7; end P7=P7/10000; w6=mR7*P7'; for n=1:10000 sum72=0; for k=0:5 if n-k>=1 sum72=sum72+w6(k+1)*v2(n-k); end end v1hat6(n)=sum72; end e6=x-v1hat6; figure(1); d=plot (d(1:200)); set(d,'color','red'); 16 hold all e6=plot(e6(1:200)); set(e6,'color','blue') title('d(n) and e(n) for first 200 samples') xlabel('n'); legend('d','e6') %problem no.8 M=12 for k8=0:11 sum8=0; for n=1:10000-k8 sum8=sum8+v2(n)*v2(n+k8); end R8(k8+1)=sum8; end R8=R8/10000 mR8=[R8(1) R8(2) R8(3) R8(4) R8(5) R8(6) R8(7) R8(8) R8(9) R8(10) R8(11) R8(12), R8(2) R8(1) R8(2) R8(3) R8(4) R8(5) R8(6) R8(7) R8(8) R8(9) R8(10) R8(11), R8(3) R8(2) R8(1) R8(2) R8(3) R8(4) R8(5) R8(6) R8(7) R8(8) R8(9) R8(10), R8(4) R8(3) R8(2) R8(1) R8(2) R8(3) R8(4) R8(5) R8(6) R8(7) R8(8) R8(9), R8(5) R8(4) R8(3) R8(2) R8(1) R8(2) R8(3) R8(4) R8(5) R8(6) R8(7) R8(8), R8(6) R8(5) R8(4) R8(3) R8(2) R8(1) R8(2) R8(3) R8(4) R8(5) R8(6) R8(7), R8(7) R8(6) R8(5) R8(4) R8(3) R8(2) R8(1) R8(2) R8(3) R8(4) R8(5) R8(6), R8(8) R8(7) R8(6) R8(5) R8(4) R8(3) R8(2) R8(1) R8(2) R8(3) R8(4) R8(5), R8(9) R8(8) R8(7) R8(6) R8(5) R8(4) R8(3) R8(2) R8(1) R8(2) R8(3) R8(4), R8(10) R8(9) R8(8) R8(7) R8(6) R8(5) R8(4) R8(3) R8(2) R8(1) R8(2) R8(3), R8(11) R8(10) R8(9) R8(8) R8(7) R8(6) R8(5) R8(4) R8(3) R8(2) R8(1) R8(2), R8(12) R8(11) R8(10) R8(9) R8(8) R8(7) R8(6) R8(5) R8(4) R8(3) R8(2) R8(1)] %since v2(n) only has real part, inverse matrix R8 is equal to R7 for k82=0:11 sum8=0; for n=1:10000-k82 sum8=sum8+v2(n)*x(n+k82); end P8(k82+1)=sum8; end P8=P8/10000 w12=mR8*P8' for n=2:10000 sum82=0; for k=1:12 if n-k>=1 sum82=sum82+w12(k)*v2(n-k); 17 end end v1hat12(n)=sum82; %e12(n)=x(n)-v1hat12(n); end e12=x-v1hat12; figure(2); d=plot (d(1:200)); set(d,'color','red'); hold all e12=plot(e12(1:200)); set(e12,'color','blue') title('d(n) and e(n) for first 200 samples') xlabel('n'); legend('d','e12')