现代信号处理及其应用仿真.doc

仿真作业周子超 2008200030433.17(1)clear allclose allf1=0.15;f2=0.17;f3=0.26;v=randn(1,256);A1=10(30/20);A2=10(30/20);A3=10(30/27);t=0:255;x1=A1*exp(j*2*pi*f1*t);x2=A2*exp(j*2*pi*f2*t);x3=A3*exp(j*2*pi*f3*t);u=v+x1+x2+x3;u=u(:,1:32);u1=u,zeros(1,32);uk=fft(u1);for i=1:64uuk(i)=(1/32)*uk(i)*conj(uk(i);endr0=ifft(uuk);rr1=r0(1,34:64),r0(1,1:32);r=0;r10=0;r1=zeros(1,63);u=u,zeros(1,32);for n=1:32r=(1/32)*u(n)*conj(u(n);r10=r+r10;endfor m=1:31 for n=m+1:32 r=(1/32)*u(n)*conj(u(n-m); r1(1,m)=r1(1,m)+r; endendfor m=32:62 for n=1:31 r=(1/32)*u(n)*conj(u(n+m-31); r1(1,m)=r1(1,m)+r; endendrr2=r1(1,62:-1:32),r10,r1(1,1:31);通过上面的计算可以发现基于FFT的自相关函数快速算法估计的自相关函数和式(3.1.2)估计出的自相关函数相等。(2)clear allclose allf1=0.15;f2=0.17;f3=0.26;v=randn(1,256);A1=10(30/20);A2=10(30/20);A3=10(27/20);t=0:255;x1=A1*exp(j*2*pi*f1*t);x2=A2*exp(j*2*pi*f2*t);x3=A3*exp(j*2*pi*f3*t);u=v+x1+x2+x3;U=fft(u);for i=1:256 SPER(i)=abs(U(i)*conj(U(i)/256;endSPER=10*log(SPER);n=0:length(u)-1;k=n/(length(u)-1);figureplot(k,SPER)xlabel(归一化频率)ylabel(功率谱(dB)title(周期图法)u=u(:,1:256);r=0;r10=0;r1=zeros(1,511);u=u,zeros(1,256);for n=1:256r=(1/256)*u(n)*conj(u(n);r10=r+r10;endfor m=1:255 for n=m+1:256 r=(1/256)*u(n)*conj(u(n-m); r1(1,m)=r1(1,m)+r; endendfor m=256:510 for n=1:255 r=(1/256)*u(n)*conj(u(n+m-255); r1(1,m)=r1(1,m)+r; endendrr=r1(1,510:-1:256),r10,r1(1,1:255);rr=rr(1,193:319);k=-63:63;w=pi*k/63+pi;Sbtw=rr*(exp(-j).(k*w);SBTW=abs(Sbtw);P=10*log(SBTW);figureplot(w/(2*pi),P)xlabel(归一化频率)ylabel(功率谱(dB)title(BT法) (3)clear allclose allf1=0.15;f2=0.17;f3=0.26;v=randn(1,256);A1=10(30/20);A2=10(30/20);A3=10(27/20);t=0:255;x1=A1*exp(j*2*pi*f1*t);x2=A2*exp(j*2*pi*f2*t);x3=A3*exp(j*2*pi*f3*t);u=v+x1+x2+x3;r=0;r10=0;r1=zeros(1,511);u=u,zeros(1,256);for n=1:256r=(1/256)*u(n)*conj(u(n);r10=r+r10;endfor m=1:255 for n=m+1:256 r=(1/256)*u(n)*conj(u(n-m); r1(1,m)=r1(1,m)+r; endendr00=0;r11=0;a=zeros(16,16);delta=zeros(16);r=r10,r1(1,1:16);a(1,1)=-r(2)/r(1);delta(1)=r(1)-r(2)*conj(r(2)/r(1);for m=2:16 r11=0; for i=1:m-1 r00=a(i,m-1)*r(m+1-i); r11=r11+r00; end k=-(r(m+1)+r11)/delta(m-1); a(m,m)=k; for i=1:m-1 a(i,m)=a(i,m-1)+k*conj(a(m-i,m-1); end delta(m)=delta(m-1)*(1-abs(k)2);endap=a(:,16).;delta=delta(16);ap=1,ap,zeros(1,1007);i=1;for w=0:2*pi/1023:2*pi e=exp(-j*w.*(0:1023); AP=ap*e.; SARW(i)=delta/(1+AP)*conj(1+AP); i=i+1;endP=abs(10*log10(SARW);w=0:1/1023:1;figureplot(w,P)xlabel(归一化频率)ylabel(功率谱(dB)title(Levinson-Durbin迭代算法)3.20clear allclose allM=8;N=1000;n=1:N;x=exp(j*2*pi*0.25*n)+exp(j*2*pi*0.6*n)+randn(size(n);Rxx=zeros(M,M); for n=M:N X=zeros(1,M); X=x(n:-1:n-M+1); Rxx=Rxx+X*X; endRxx=(1/(N-M+1)*Rxx; V D=eig(Rxx,nobalance);weights=diag(D);weights=weights(1:6);for k=1:6; noise_eigenvects(:,k)=V(:,k);end% root-musicP=0;for i=1:length(weights), P=P+conv(noise_eigenvects(:,i),conj(flipud(noise_eigenvects(:,i);endroots_P=roots(P);roots_P1=roots_P(abs(roots_P)0) f1rootmusic=angle(sorted_roots(1)/(2*pi);else f1rootmusic=(angle(sorted_roots(1)+2*pi)/(2*pi);endif(angle(sorted_roots(2)0) f2rootmusic=angle(sorted_roots(2)/(2*pi);else f2rootmusic=(angle(sorted_roots(2)+2*pi)/(2*pi);endfrootmusic=f1rootmusic f2rootmusic;% musici1=1;for w=0:2*pi/1023:2*piaw=1 exp(j*w) exp(j*2*w) exp(j*3*w) exp(j*4*w) exp(j*5*w) exp(j*6*w) exp(j*7*w).;Pmusic(i1)=1/abs(aw*noise_eigenvects*noise_eigenvects*aw); Pmusic(i1)=10*log10(Pmusic(i1);i1=i1+1;endw=0:1/1023:1;plot(w,Pmusic)xlabel(归一化频率)ylabel(功率谱(dB)title(Music算法)通过上面的计算，Root-Music算法估计的频率为f1=0.25，f2=0.5999，而Music算法的估计结果为f1=0.248，f2=0.63，所以Root-Music算法估计的频率较Music算法估计的准确。3.21clearcloseN=1000;M=8;delta=0; I=eye(M,M);Z=zeros(M,M); for i=1:M-1 Z(i,i+1)=1;endw=randn(1,1000) ; n=1:1000; x(n)=exp(j*0.5*pi*n)+exp(j*1.2*pi*n+j*4/pi)+w(1,n); Rxx=zeros(M,M); Rxy=zeros(M,M);for n=1:N-M+1 X=zeros(1,M); X=x(n:n+M-1); Rxx=Rxx+X*X; endfor n=1:N-M X=zeros(1,M); Y=zeros(1,M); X=x(n:n+M-1); Y=x(n+1:n+M); Rxy=Rxy+X*Y; endRxx=(1/(N-M+1)*Rxx; Rxy=(1/(N-M)*Rxy; Cxx=zeros(M,M); Cxy=zeros(M,M);V1,D1=eig(Rxx);V2,D2=eig(Rxy);delta=min(diag(D1); Cxx=Rxx-delta*I; Cxy=Rxy-delta*Z; V,D=eig(Cxx,Cxy);lamda=diag(D);lamda1=abs(lamda);error=abs(diag(I)-lamda1);if(-angle(D(1,1)0) f2=-angle(D(2,2)/(2*pi)else f2=(-angle(D(2,2)+2*pi)/(2*pi)end通过上面的计算，ESPRIT算法估计的信号频率为f1=0.2503，f2=0.6112.4.18(1)clear allclose allN=512; p=0.0731; a=sqrt(p);v=a*randn(N,1); num=1 0 0;den=1 -0.975 0.95;Gz=tf(num,den,-1);u=lsim(Gz,v); (2)clear allclose allN=512; p=0.0731; a=sqrt(p);v=a*randn(N,1); num=1 0 0;den=1 -0.975 0.95;Gz=tf(num,den,-1);u=lsim(Gz,v); d=u; mu=0.05;% mu=0.005;w=zeros(2,1);for n=1:N if(n=1) x=0;0; else if(n=2) x=u(1);0; else x=u(n-1:-1:n-2); end end e(n)=d(n)-w*x; w=w+mu*e(n)*x; Recordedw(1:2,n)=w; ee(1,n)=e(n)2;end 通过上面的计算，=0.05时，w1=0.9885，w2=-0.9507；=0.005时，w1=0.7665，w2=-0.7361(采样点数太少，w1和w2还没有收敛)。(3)、(4)clear allclose allN=512;g=100;a1=-0.975;a2=0.95;mu=0.05;% mu=0.005;pp=zeros(g,N);qq=zeros(g,N);ppp=zeros(1,N);qqq=zeros(1,N);ee=zeros(g,N);eee=zeros(1,N);for q=1:gp=0.0731; a=sqrt(p);v=a*randn(N,1); num=1 0 0;den=1 -0.975 0.95;Gz=tf(num,den,-1);u=lsim(Gz,v); d=u;w=zeros(2,1);for n=1:N if(n=1) x=0;0; else if(n=2) x=u(1);0; else x=u(n-1:-1:n-2); end end e(n)=d(n)-w*x; w=w+mu*e(n)*x; Recordedw(1:2,n)=w; e(1,n)=e(n)2;endfor n=1:N ee(q,n)=e(1,n); pp(q,n)=Recordedw(1,n); qq(q,n)=Recordedw(2,n);endfor q=1:N eee(1,q)=sum(ee(:,q)/g; ppp(1,q)=sum(pp(:,q)/g; qqq(1,q)=sum(qq(:,q)/g; endendfiguresemilogy(eee,g);title(平均100次的均方误差曲线);xlabel(Samples);figurehold onplot(ppp,r);plot(qqq,b);legend(w1(n)的变化曲线,w2(n)的变化曲线,1);title(平均100次的变化曲线);xlabel(Samples);ylabel(Weights value);hold offP=zeros(2,1);Rxx=zeros(2,2);w=ppp(1,512);qqq(1,512);for n=3:N x=u(n-1:-1:n-2); Rxx=Rxx+x*x; P=P+x*u(n);endRxx=Rxx/(N-3+1);P=P/(N-3+1);daltau=(1+a2)/(1-a2)*(1+a2)2-a12);Jmin=daltau-P*w;V D=eig(Rxx);Jex=mu*Jmin*(D(1,1)+D(2,2)/(2-mu*(D(1,1)+D(2,2);M=Jex/Jmin;=0.05通过上面的计算，此时剩余均方误差Jex=0.5902，失调参数M=0.0458.=0.005通过上面的计算，此时剩余均方误差Jex=0.0614，失调参数M=0.0047.但是通过上图可知由于512采样点数太少，此时w1和w2还没有收敛，故计算结果不正确。6.13clear allclose allf1=0.1;f2=0.25;f3=0.27;N=1000;n=1:N;v=randn(1,N);A1=10(30/20);A2=10(30/20);A3=10(27/20);x1=A1*exp(j*2*pi*f1*n);x2=A2*exp(j*2*pi*f2*n);x3=A3*exp(j*2*pi*f3*n);u=v+x1+x2+x3;for i=4:N AH(:,i-3)=u(i:-1:i-4+1);endA=AH;U,S,V=svd(A);P=0;m=1;s=zeros(1,4);ss=zeros(1,4); for w=0:2*pi/1023:2*pi aw=1 exp(-j*w) exp(-j*2*w) exp(-j*3*w).;for i=1:4 s(1,i)=aw*V(:,i)/S(i,i); ss(1,i)=s(1,i)*conj(s(1,i); P=P+ss(1,i);end Pmvdr(1,m)=1/P; m=m+1; P=0;endw=0:1/1023:1;plot(w,Pmvdr)xlabel(归一化频率)ylabel(功率谱) 通过上图可知，f1=0.1可以分辨出来，但是对于抽头个数为4太少，故f2=0.25和f3=0.27不能够分辨出来。6.15clear allclose all% Number of signal pointsm=2000;g=100;%Take only N points for training N=100 ;%forgetting factorlamda=0.98; delta=1/0.02; pp=zeros(g,N);qq=zeros(g,N);ppp=zeros(1,N);qqq=zeros(1,N);ee=zeros(g,N);eee=zeros(1,N);for q=1:gp=0.995; a=sqrt(p);u=a*randn(m,1); num=1 0;den=1 0.99;Gz=tf(num,den,-1);out=lsim(Gz,u); d=out;%initial P matrixP=delta*eye(2);w=zeros(2,1);for n=1:N if(n=1) x=0;0; else if(n=2) x=out(1);0; else x=out(n-1:-1:n-2); end end phi=x*P; k=phi/(lamda+phi*x); y(n)=w*x; e(n)=d(n)-y(n); w=w+k*e(n); P=(P-k*phi)/lamda; % Just for plotting Recordedw(1:2,n)=w; e(1,n)=e(n)2;endfor n=1:N ee(q,n)=e(1,n); pp(q,n)=Recordedw(1,n); qq(q,n)=Recordedw(2,n);endfor q=1:N eee(1,q)=sum(ee(:,q)/g; ppp(1,q)=sum(pp(:,q)/g; qqq(1,q)=sum(qq(:,q)/g;endendhold onplot(eee,g);plot(ppp,r);plot(qqq,b);legend(e(n)的变化曲线,w1(n)的变化曲线,w2(n)的变化曲线,1);title(平均100次的变化曲线);xlabel(Samples);ylabel(Weights value);hold off7.14clear allclose all% Number of signal pointsm=2000;v=0.0332; %a1=sqrt(v);u=a1*randn(m,1); %num=1 0 0 0 0;den=1 -1.595 0.95 -0.84 1.211;Gz=tf(num,den,-1);out=lsim(Gz,u); %d=out; %Take only N points for training N=300 ; %forgetting factor% lamda=0.9; %initial P matrix% delta =50; p=eye(4);w=1;0;0;0;for n=1:N if(n=1) x=0;0;0;0; else if(n=2) x=out(1);0;0;0; else if(n=3) x=out(2);out(1);0;0; else if(n=4) x=out(3);out(2);out(1);0; else x=out(n-1);out(n-2);out(n-3);out(n-4); end end end end d1=x.*w; a=d(n)-d1; A=x.*p*x+v; k=p*x*A(-1); p=p-k*x.*p; w=w+k*a; Recordedw(1:4,n)=w;end figurehold onplot(Recordedw(1,1:N),r);plot(Recordedw(2,1:N),b);plot(Recordedw(3,1:N),g);plot(Recordedw(4,1:N),y);legend(w1(n)的变化曲线,w2(n) w1(n)的变化曲线,w3(n) w1(n)的变化曲线,w4(n) w1(n)的变化曲线,1);title(单次运算的变化曲线);xlabel(Samples);ylabel(Weights value);hold off通过上图可以看出，用卡尔曼滤波算法估计的最优权值和原模型参数接近，但是通过多次运行发现，该方法估计的权值波动很大。8.16(1)clear; clcclose all;k=1:100;f1=0.1; f2=0.25;theta=-90:0.01:90;s1=sqrt(10)*sin(2*pi*f1*k); s2=10*sin(2*pi*f2*k);S=s1; w1=pi*sin(2*pi/9); w2=pi*sin(-pi/18);a1=1 exp(j*w1) exp(j*2*w1) exp(j*3*w1) exp(j*4*w1) exp(j*5*w1) exp(j*6*w1) exp(j*7*w1) exp(j*8*w1) exp(j*9*w1).; A1=a1;a2=1 exp(j*w2) exp(j*2*w2) exp(j*3*w2) exp(j*4*w2) exp(j*5*w2) exp(j*6*w2) exp(j*7*w2) exp(j*8*w2) exp(j*9*w2).;A2=a2;noise=zeros(10,100);for i=1:10noise(i,:)=randn(1,100);endU=zeros(10,100);for i=1:100 U(:,i)=A1*s1(:,i)+A2*s2(:,i)+noise(:,i);endcovmatrix=zeros(10,10);for i=1:100temp=U(:,i)*(U(:,i);covmatrix=(covmatrix+temp)/100;endV D=eig(covmatrix,nobalance);for k=1:6;g(:,k)=V(:,k);endi1=1;for theta1=-90:0.01:90aw=1 exp(j*pi*sin(theta1/180*pi) exp(j*pi*2*sin(theta1/180*pi) exp(j*pi*3*sin(theta1/180*pi) exp(j*pi*4*sin(theta1/180*pi) exp(j*pi*5*sin(theta1/180*pi) exp(j*pi*6*sin(theta1/180*pi) exp(j*pi*7*sin(theta1/180*pi) exp(j*pi*8*sin(theta1/180*pi) exp(j*pi*9*sin(theta1/180*pi).;Pmusic(i1)=1/abs(aw*g*g*aw); Pmusic(i1)=10*log10(Pmusic(i1);i1=i1+1;endplot(theta,Pmusic)xlabel(空间角度)ylabel(Music谱(dB)(2)clear;clcclose allk=1:100;f1=0.1; f2=0.25;theta=-90:0.01:90;s1=sqrt(10)*exp(j*2*pi*f1*k); s2=10*exp(j*2*pi*f2*k);S=s1; w1=pi*sin(2*pi/9); w2=pi*sin(-pi/18);a1=1 exp(j*w1) exp(j*2*w1) exp(j*3*w1) exp(j*4*w1) exp(j*5*w1) exp(j*6*w1) exp(j*7*w1) exp(j*8*w1) exp(j*9*w1).; A1=a1;a2=1 exp(j*w2) exp(j*2*w2) exp(j*3*w2) exp(j*4*w2) exp(j*5*w2) exp(j*6*w2) exp(j*7*w2) exp(j*8*w2) exp(j*9*w2).;A2=a2; noise=zeros(10,100);for i=1:10noise(i,:)=randn(1,100);endU=zeros(10,100);for i=1:100U(:,i)=A1*s1(:,i)+A2*s2(:,i)+noise(:,i);endcovmatrix=zeros(10,10);for i=1:100temp=U(:,i)*(U(:,i);covmatrix=(covmatrix+temp)/100;endV D=eig(covmatrix,nobalance);weights=diag(D);weights=weights(1:8);for k=1:8; noise_eigenvects(:,k)=V(:,k);endP=0;for i=1:length(weights), P=P+conv(noise_eigenvects(:,i),conj(flipud(noise_eigenvects(:,i);endroots_P=roots(P);roots_P1=roots_P(abs(roots_P)1);not_used,indx=sort(abs(abs(roots_P1)-1); sorted_roots=roots_P1(indx);fai=angle(sorted_roots(1:2)/pi;theta=asin(fai)*180/pi;通过上面的计算，两个不相干的信号源的波达方向为-10.0688和39.2571。(3)clear; clcclose all;k=1:100;f1=0.1; f2=0.25;theta=-90:0.01:90;s1=sqrt(10)*exp(j*2*pi*f1*k); s2=10*exp(j*2*pi*f2*k);w1=pi*sin(2*pi/9);w2=pi*sin(-pi/18);a1=1 exp(j*w1) exp(j*2*w1) exp(j*3*w1) exp(j*4*w1) exp(j*5*w1) exp(j*6*w1) exp(j*7*w1) exp(j*8*w1) exp(j*9*w1); %A1=a1.;a2=1 exp(j*w2) exp(j*2*w2) exp(j*3*w2) exp(j*4*w2) exp(j*5*w2) exp(j*6*w2) exp(j*7*w2) exp(j*8*w2) exp(j*9*w2);A2=a2.;noise=zeros(10,100);for i=1:10noise(i,:)=randn(1,100);endU=zeros(10,100);for i=1:100 U(:,i)=A1*s1(1,i)+A2*s2(1,i)+noise(:,i);endcovmatrix=zeros(10,10);for i=1:100temp=U(:,i)*(U(:,i);covmatrix=(covmatrix+temp)/100;endV D=eig(covmatrix,nobalance);for k=9:10;s(:,k-8)=V(:,k);ends1=s(1:9,:);s2=s(2:10,:);s12=s1 s2;ss=s12*s12;U,S,V=svd(ss);u11=U(1:2,1:2);u12=U(1:2,3:4);u21=U(3:4,1:2);u22=U(3:4,3:4);fai=(-1)*u12*i