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

hold off;

%d(n),p=3,6,9,12 for 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=A\\B;

A2=toeplitz(rv2(1:p(2)+1)',rv2(1:p(2)+1)); B2=rxv2(1:p(2)+1)'; w2=A2\\B2;

A3=toeplitz(rv2(1:p(3)+1)',rv2(1:p(3)+1)); B3=rxv2(1:p(3)+1)'; w3=A3\\B3;

A4=toeplitz(rv2(1:p(4)+1)',rv2(1:p(4)+1)); B4=rxv2(1:p(4)+1)'; w4=A4\\B4; %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); 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; w1' w2' w3' w4'

N=500;

g=normrnd(0,1,1,N); v1=filter([1],[1,-0.8],g); v2=filter([1],[1,0.6],g); d=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:13 rv2(k)=0; for n=k:N

rv2(k)=rv2(k)+v0(n)*v0(n-(k-1)); end

rv2(k)=rv2(k)/N; end %

for k=1:13

rxv2(k)=0; for n=k:N

rxv2(k)=rxv2(k)+x(n)*v0(n-(k-1)); end

rxv2(k)=rxv2(k)/N; end %

pp=12;

A=toeplitz(rv2(1:pp+1)',rv2(1:pp+1)); B=rxv2(1:pp+1)'; w1=A\\B; %

v11=filter(w1',[1],v0);

d1=[d1 x-v11]; (v1-v11)*(v1-v11)'/N end

figure(7)

plot(0:N-1,d1(1:N)); hold on; plot(0:N-1,d,'r');

title('辅助观测数据中漏入d[n]的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('辅助观测数据中漏入d[n]的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('辅助观测数据中漏入d[n]的1.0'); legend('估计结果','期望值'); axis([0 N -4 4]); grid on; hold off;