现代数字信号处理仿真作业

第三章仿真作业3.17（1）代码clear;N=32;m=[-N+1:N-1];noise=(randn(1,N)+j*randn(1,N))/sqrt(2);f1=0.15;f2=0.17;f3=0.26;SNR1=30;SNR2=30;SNR3=27;A1=10^(SNR1/20);A2=10^(SNR2/20);A3=10^(SNR3/20);signal1=A1*exp(j*2*pi*f1*(0:N-1));signal2=A2*exp(j*2*pi*f2*(0:N-1));signal3=A3*exp(j*2*pi*f3*(0:N-1));un=signal1+signal2+signal3+noise;uk=fft(un,2*N);sk=(1/N)*abs(uk).^2;r0=ifft(sk);r1=[r0(N+2:2*N),r0(1:N)];r=xcorr(un,N-1,'biased');figuresubplot(2,2,1)stem(m,real(r1));xlabel('m');ylabel('FFT估计r1实部');subplot(2,2,2)stem(m,imag(r1));xlabel('m');ylabel('FFT估计r1虚部');subplot(2,2,3)stem(m,real(r));xlabel('m');ylabel('平均估计r实部');subplot(2,2,4)stem(m,imag(r));xlabel('m');ylabel('平均估计r虚部');仿真结果（2）代码clear;N=256;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);f1=0.15;f2=0.17;f3=0.26;SNR1=30;SNR2=30;SNR3=27;A1=10^(SNR1/20);A2=10^(SNR2/20);A3=10^(SNR3/20);signal1=A1*exp(j*2*pi*f1*(0:N-1));signal2=A2*exp(j*2*pi*f2*(0:N-1));signal3=A3*exp(j*2*pi*f3*(0:N-1));un=signal1+signal2+signal3+noise;-40-2002040-2000020004000mFFT估计r1实部-40-2002040-2000-1000010002000mFFT估计r1虚部-40-2002040-2000020004000m平均估计r实部-40-2002040-2000-1000010002000m平均估计r虚部NF=1024;spr=fftshift((1/NF)*abs(fft(un,NF)).^2);f1=(0:length(spr)-1)*(1/(length(spr)-1))-0.5;M=64;r=xcorr(un,M,'biased');bt=fftshift(abs(fft(r,NF)));f2=(0:length(bt)-1)*(1/(length(bt)-1))-0.5;figuresubplot(1,2,1)plot(f1,10*log10(spr/max(spr)));xlabel('w/2pi');仿真结果（3）代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise;p=8;cx=xcorr(un,p,'biased');rxx=cx(p+1:2*p)';R=toeplitz(rxx);[u,s]=eig(R);-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5-80-70-60-50-40-30-20-100w/2pi-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5-45-40-35-30-25-20-15-10-50w/2pinw=128;ww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k行m列ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果3.20（1）代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise;p=8;cx=xcorr(un,p,'biased');rxx=cx(p+1:2*p)';R=toeplitz(rxx);[u,s]=eig(R);nw=128;-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50.40.60.811.21.41.61.822.22.4x10-4ww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k行m列ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果距离单位圆最近的两个解为-0.2363-0.9717i和0.3747+0.9271i，对应的归一化频率为0.1889和-0.2880（2）代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise;p=8;cx=xcorr(un,p,'biased');rxx=cx(p+1:2*p)';R=toeplitz(rxx);[u,s]=eig(R);nw=128;ww=[-128:128]/128*pi;e=exp(-j*ww'*[0:p-1])%k行m列ev=e*u(:,1:p-2);pw=1./real(diag(ev*ev'));plot(ww/(2*pi),10*log10(pw)/max(pw));仿真结果-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5-2-10123456x10-33.21代码clear;N=1000;fai1=rand(1,1)*2*pi;fai2=rand(1,1)*2*pi;noise=(randn(1,N)+j*randn(1,N))/sqrt(2);un=exp(j*0.5*pi*(0:N-1)+j*fai1)+exp(-j*0.3*pi*(0:N-1)+j*fai2)+noise;p=8;fork=1:N-pxs(:,k)=un(k+p-1:-1:k)';endrxx=xs(:,1:end-1)*xs(:,1:end-1)'/(N-p-1);rxy=xs(:,1:end-1)*xs(:,2:end)'/(N-p-1);[u,e]=svd(rxx);ev=diag(e);emin=ev(end);z=[zeros(p-1,1),eye(p-1);0,zeros(1,p-1)];cxx=rxx-emin*eye(p);cxy=rxy-emin*z;[U,E]=eig(cxx,cxy);Z=diag(E);ph=angle(Z)/(2*pi);err=abs(abs(Z)-1);仿真结果最接近单位圆的两个解分别为0.5867+0.8097i和0.0349-0.9984i，对应的归一化频率为0.1502和-0.2444。第四章仿真题作业4-18（1）产生随机序列代码clear;data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;vn=sqrt(sigma_v_2)*randn(data_len,1,trials);u0=[00];num=1;den=[1a1a2];zi=filtic(num,den,u0);u=filter(num,den,vn,zi);(2)单次实验估计最优权向量代码clear;data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;vn=sqrt(sigma_v_2)*randn(data_len,1,trials);u0=[00];num=1;m=1;den=[1a1a2];zi=filtic(num,den,u0);u=filter(num,den,vn,zi);mu1=0.05;mu2=0.005;w1=zeros(2,data_len);w2=zeros(2,data_len);e1=zeros(data_len,1);e2=zeros(data_len,1);d1=zeros(data_len,1);d2=zeros(data_len,1);%form=1:trialsforn=3:data_len-1w1(:,n+1)=w1(:,n)+mu1*u(n-1:-1:n-2,:,m)*conj(e1(n));w2(:,n+1)=w2(:,n)+mu2*u(n-1:-1:n-2,:,m)*conj(e2(n));d1(n+1)=w1(:,n+1)'*u(n:-1:n-1,:,m);d2(n+1)=w2(:,n+1)'*u(n:-1:n-1,:,m);e1(n+1)=u(n+1,:,m)-d1(n+1);e2(n+1)=u(n+1,:,m)-d2(n+1);endfigure(1)holdon;plot(w1(1,:));plot(w1(2,:));plot(w2(1,:),’g’);plot(w2(2,:),’g’)仿真结果（3）100次实验计算剩余均方误差代码clear;data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;vn=sqrt(sigma_v_2)*randn(data_len,1,trials);0100200300400500600-1.5-1-0.500.511.5u0=[00];num=1;%m=1;den=[1a1a2];zi=filtic(num,den,u0);u=filter(num,den,vn,zi);mu1=0.05;mu2=0.005;w1=zeros(2,data_len);w2=zeros(2,data_len);e1=zeros(data_len,1);e2=zeros(data_len,1);d1=zeros(data_len,1);d2=zeros(data_len,1);jn1=zeros(1,trials);jn2=zeros(1,trials);mse1=zeros(1,trials);mse2=zeros(1,trials);form=1:trialsforn=3:data_len-1w1(:,n+1)=w1(:,n)+mu1*u(n-1:-1:n-2,:,m)*conj(e1(n));w2(:,n+1)=w2(:,n)+mu2*u(n-1:-1:n-2,:,m)*conj(e2(n));d1(n+1)=w1(:,n+1)'*u(n:-1:n-1,:,m);d2(n+1)=w2(:,n+1)'*u(n:-1:n-1,:,m);e1(n+1)=u(n+1,:,m)-d1(n+1);e2(n+1)=u(