《现代数字信号处理及其应用》仿真作业

第三章仿真作业3.17(1)N=32;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(1i*2*pi*f1*(0:N-1));signal2=A2*exp(1i*2*pi*f2*(0:N-1));signal3=A3*exp(1i*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');r11=real(r1);r12=imag(r1);r1=real(r);r2=imag(r);m=1-N:N-1;subplot(2,2,1);stem(m,r11,'o');xlabel('m');ylabel('实部');title('基于FFT的自相关函数快速计算');subplot(2,2,2);stem(m,r12,'o');xlabel('m');ylabel('虚部');subplot(2,2,3);stem(m,r1);xlabel('m');ylabel('实部');title('教材式（3.1.2）估计的自相关函数');subplot(2,2,4);stem(m,r2);xlabel('m');ylabel('虚部');（2）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(1i*2*pi*f1*(0:N-1));signal2=A2*exp(1i*2*pi*f2*(0:N-1));signal3=A3*exp(1i*2*pi*f3*(0:N-1));un=signal1+signal2+signal3+noise;NF=1024;Spr=fftshift((1/NF)*abs(fft(un,NF)).^2);Sprmax=max(Spr);Spr=Spr/Sprmax;f=(-(NF/2)+1:(NF/2))/NF;plot(f,20*log(Spr));xlabel('f');ylabel('归一化功率谱/dB');title('周期图法');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(1i*2*pi*f1*(0:N-1));signal2=A2*exp(1i*2*pi*f2*(0:N-1));signal3=A3*exp(1i*2*pi*f3*(0:N-1));un=signal1+signal2+signal3+noise;M=64;r=xcorr(un,N-1,'biased');NF=1024;BT=fftshift(fft(r,NF));BTmax=max(BT);BT=BT/BTmax;%归一化幅度f=(-(NF/2)+1:(NF/2))/NF;plot(f,20*log(BT));xlabel('f');ylabel('归一化功率谱/dB');title('BT法');（3）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(1i*2*pi*f1*(0:N-1));signal2=A2*exp(1i*2*pi*f2*(0:N-1));signal3=A3*exp(1i*2*pi*f3*(0:N-1));un=signal1+signal2+signal3+noise;p=16;r0=xcorr(un,p,'biased');r=r0(p+1:2*p+1);%从p+1开始取a(1,1)=-r(2)/r(1);sigma(1)=r(1)-(abs(r(2)).^2)/r(1);form=2:pk(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)))/sigma(m-1);a(m,m)=k(m);fori=1:m-1;a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i));endsigma(m)=sigma(m-1)*(1-abs(k(m).^2));endNF=1024;Par=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);Parmax=max(Par);Par=Par/Parmax;%归一化幅度f=(-(NF/2)+1:(NF/2))/NF;plot(f,20*log(Par));xlabel('f');ylabel('归一化功率谱/dB');title('16阶AR模型的功率谱估计');3.20（1）clc;clearall;closeall;N=1000;noise=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal1=exp(j*0.5*pi*(0:N-1)+j*unifrnd(0,2*pi,1,1));signal2=exp(-j*0.3*pi*(0:N-1)+j*unifrnd(0,2*pi,1,1));un=signal1+signal2+noise;M=8;xs=zeros(M,N-M);fork=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endR=(xs*xs')/(N-M);[U,E,V]=svd(R);G=U(:,3:M);Gr=G*G';co=zeros(2*M-1,1);form=1:Mco(m:m+M-1)=co(m:m+M-1)+Gr(M:-1:1,m);endz=roots(co);ph=angle(z)/(2*pi);z1=abs(z);z2=abs(z1-1);[estallv,estain]=sort(z2,'ascend');所以单Root-MUSIC算法中最近单位圆的两个根为0.0033-0.9977i0.5856+0.8074i上述根的相位对应的归一化频率为0.2495-0.1501（2）clc;clearall;closeall;N=1000;noise=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal1=exp(j*0.5*pi*(0:N-1)+j*2*pi*rand);signal2=exp(-j*0.3*pi*(0:N-1)+j*2*pi*rand);un=signal1+signal2+noise;M=8;xs=zeros(M,N-M);fork=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endR=(xs*xs')/(N-M);[U,E,V]=svd(R);G=U(:,3:M);f=[-0.5:1/999:0.5];forff=1:length(f)w=f(ff)*2*pi;forl=1:Maw(l)=exp(-j*w*(l-1));%计算a(w)endWW=aw*G*G'*aw';Pmusic(ff)=abs(1./WW);%谱扫描函数endPmusic=10*log10(Pmusic);f=[-0.5:1/999:0.5];plot(f,Pmusic);holdon%绘图输出xlabel('ω/2π')ylabel('归一化Music谱/dB');3.21clc;clearall;closeall;N=1000;noise=(randn(1,N)+1i*randn(1,N))/sqrt(2);signal1=exp(j*0.5*pi*(0:N-1)+j*unifrnd(0,2*pi,1,1));signal2=exp(-j*0.3*pi*(0:N-1)+j*unifrnd(0,2*pi,1,1));un=signal1+signal2+noise;M=8;xs=zeros(M,N-M);fork=1:N-Mxs(:,k)=un(k+M-1:-1:k).';endRxx=xs(:,1:length(xs)-1)*xs(:,1:length(xs)-1)'/(N-M-1);Rxy=xs(:,1:length(xs)-1)*xs(:,2:length(xs))'/(N-M-1);[U,E]=svd(Rxx);emin=min(ev);Z=[zeros(M-1,1),eye(M-1);0,zeros(1,M-1)];Cxx=Rxx-emin*eye(M);Cxy=Rxy-emin*Z;[U,E]=eig(Cxx,Cxy);z=diag(E);ph=angle(z)/(2*pi);z1=abs(z);z2=abs(z1-1);[estallv,estain]=sort(z2,'ascend');所以单次ESPRIT算法中最近单位圆的两个特征值为0.5845+0.8117i0.0036-0.9994i上述特征值的相位对应的归一化频率为-0.15070.2494第四章仿真作业4.18%（1）产生N=512点的样本序列data_len=512;%样本序列的长度trials=100;%随机试验的次数A=zeros(data_len,2);EA=zeros(data_len,1);B=zeros(data_len,2);EB=zeros(data_len,1);form=1:trialsa1=-0.975;a2=0.95;sigma_v_2=0.0731;v=sqrt(sigma_v_2)*randn(data_len,1,trials);%产生v(n)u0=[00];num=1;den=[1a1a2];Zi=filtic(num,den,u0);%滤波器的初始条件u=filter(num,den,v,Zi);%产生样本序列u(n)%（2）用LMS滤波器来估计w1和w2mu1=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);%LMS迭代过程forn=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);endA=A+conj(w1)';EA=EA+e1.^2;B=B+conj(w2)';EB=EB+e2.^2;end%剩余均方误差和失调参数wopt=zeros(2,trials);Jmin=zeros(1,trials);sum_eig=zeros