您好,欢迎访问三七文档
当前位置:首页 > 电子/通信 > 综合/其它 > 现代数字信号处理仿真作业
3.17N=32;noise=(randn(1,N)+j*rand(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');subplot(2,2,1);m=-31:31;stem(m,real(r1));title('实部');axis([-31,31,-1500,2500]);subplot(2,2,2);m=-31:31;stem(m,imag(r1));title('虚部');axis([-31,31,-2000,2000]);subplot(2,2,3);m=-31:31;stem(m,real(r));title('实部');axis([-31,31,-1500,2500]);subplot(2,2,4);m=-31:31;stem(m,imag(r));title('虚部');axis([-31,31,-2000,2000]);N=256;noise=(randn(1,N)+j*rand(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;NF=1024;Pw_p=(abs(fft(un,NF)).^2)/NF;Pw_p=10*log10(Pw_p);Pw_p=Pw_p-max(Pw_p);figure();dw=2*pi/NF;df=dw/(2*pi);plot(-pi/(2*pi):df:pi/(2*pi)-df,fftshift(Pw_p));title('周期图法');xlabel('w/2π');ylabel('归一化功率谱/dB');M=64;r=xcorr(un,M,'biased');NF=1024;Pw_bt=abs(fft(r,NF));Pw_bt=10*log10(Pw_bt);Pw_bt=Pw_bt-max(Pw_bt);figure();dw=2*pi/length(Pw_bt);df=dw/(2*pi);plot(-pi/(2*pi):df:pi/(2*pi)-df,fftshift(Pw_bt));title('BT法');xlabel('w/2π');ylabel('归一化功率谱/dB');p=16;r0=xcorr(un,p,'biased');r=r0(p+1:2*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-1a(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;Pw_LD=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);Pw_LD=10*log10(Pw_LD);Pw_LD=Pw_LD-max(Pw_LD);dw=2*pi/length(Pw_LD);df=dw/(2*pi);figure();plot(-pi/(2*pi):df:pi/(2*pi)-df,Pw_LD);title('16阶AR模型的功率谱估计');xlabel('w/2π');ylabel('归一化功率谱/dB');3.20N=1000;M=8;K=2;v=sqrt(1/2)*randn(1,N)+j*sqrt(1/2)*randn(1,N);n=1:N;fi=2*pi*rand(1,2);u(n)=exp(j*0.5*pi*n+j*fi(1))+exp(-j*0.3*pi*n+j*fi(2))+v(n);fork=1:N-M;x(:,k)=u(k+M-1:-1:k).';endR=x*x'/(N-M);NF=2048;forn=1:NFAq=exp(-j*2*pi*(n-1)/NF*(0:M-1)');Pmvdr(n)=1/(Aq'*inv(R)*Aq);%MVDR谱endPmvdr=10*log10(Pmvdr/max(Pmvdr));dw=2*pi/length(Pmvdr);df=dw/(2*pi);figure;plot(-pi/(2*pi):df:pi/(2*pi)-df,Pmvdr);title('MVDR谱估计');xlabel('w/2π');ylabel('归一化功率谱/dB');NF=2048;[G,D]=eigs(R,M-K,'sm');forn=1:NF;aw=exp(-j*2*pi*(n-1)/NF*(0:M-1)');P_music(n)=1/(aw'*G*G'*aw);endp1=abs(P_music);dw=2*pi/length(P_music);df=dw/(2*pi);figure;plot(-pi/(2*pi):df:pi/(2*pi)-df,fftshift(10*log10(p1/max(p1))));title('MUSIC谱估计');xlabel('w/2π');ylabel('归一化功率谱/dB');Gr=G*G';a=zeros(2*M-1,1);form=1:Ma(m:m+M-1)=a(m:m+M-1)+Gr(M:-1:1,m);endx=roots(a);a2=x(abs(x)1);[lmdaI]=sort(abs(abs(a2)-1));f=a2(I(1:2));so1=(angle(f(1)))/(2*pi)so2=(angle(f(2)))/(2*pi)结果So1=0.25So2=-0.14983.21N=1000;M=8;K=2;v=sqrt(1/2)*randn(1,N)+j*sqrt(1/2)*randn(1,N);n=1:N;fi=2*pi*rand(1,2);u(n)=exp(j*0.5*pi*n+j*fi(1))+exp(-j*0.3*pi*n+j*fi(2))+v(n);fork=1:N-M;x(:,k)=u(k+M-1:-1:k).';endRxx=x(:,1:end-1)*x(:,1:end-1)'/(N-M-1);Rxy=x(:,1:end-1)*x(:,2:end)'/(N-M-1);[U,E]=svd(Rxx);ev=diag(E);emin=ev(end);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);err=abs(abs(z)-1);[lmdaI]=sort(abs(abs(z)-1));f=ph(I(1:2));so1=f(1)so2=f(2)结果So1=0.2495So2=-0.15034.18data_len=512;trials=100;n=1:data_len;a1=-0.975;a2=0.95;sigma_v_2=0.0731;error1=zeros(data_len,1);error2=zeros(data_len,1);weight1=zeros(2,data_len);weight2=zeros(2,data_len);forloop=1:100v=sqrt(sigma_v_2)*randn(data_len,1,trials);u0=[00];num=1;den=[1a1a2];Zi=filtic(num,den,u0);u=filter(num,den,v,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);forn=3:data_len-1w1(:,n+1)=w1(:,n)+mu1*u(n-1:-1:n-2,:,1)*conj(e1(n));w2(:,n+1)=w2(:,n)+mu2*u(n-1:-1:n-2,:,1)*conj(e2(n));d1(n+1)=w1(:,n+1)'*u(n:-1:n-1,:,1);d2(n+1)=w2(:,n+1)'*u(n:-1:n-1,:,1);e1(n+1)=u(n+1,1)-d1(n+1);e2(n+1)=u(n+1,1)-d2(n+1);enderror1=error1+e1.^2;error2=error2+e2.^2;weight1=weight1+w1;weight2=weight2+w2;enderror1=error1./100;error2=error2./100;weight1=weight1./100;weight2=weight2./100;figure(1);plot(w1(1,:));holdon;plot(w1(2,:));holdon;plot(weight1(1,:));holdon;plot(weight1(2,:));title('步长为0.05时权向量收敛曲线');xlabel('迭代次数');ylabel('权向量');figure(2);plot(w2(1,:));holdon;plot(w2(2,:));holdon;plot(weight2(1,:));holdon;plot(weight2(2,:));title('步长为0.005时权向量收敛曲线');xlabel('迭代次数');ylabel('权向量');N=1:data_len;figure(3);plot(N,error1,'*',N,error2);title('步长分别为0.05和0.005时100次独立实验的学习曲线');legend('\mu1=0.05','\mu2=0.005');xlabel('迭代次数');ylabel('均方误差');wopt=zeros(2,trials);Jmin=zeros(1,trials);sum_eig=zeros(trials,1);form=1:trialsrm=xcorr(u(:,:,m),'biased');R=[rm(512),rm(513);rm(511),rm(512)];p=[rm(511);rm(510)];wopt(:,m)=R\p;[v,d]=eig(R);Jmin(m)=rm(512)-p'*wopt(:,m);sum_eig(m)=d(1,1)+d(2,2);ende1_100trials_ave=0.01*sum(sum_eig);e2_1
本文标题:现代数字信号处理仿真作业
链接地址:https://www.777doc.com/doc-1521688 .html