数字信号处理仿真作业.doc

作业1（1）求a = 1 -a1 -a2; 0 -(1+a2) 0; 0 -a1 -1;b = r0 a1*r0 a2*r0;x = inv(a)*b; %x = sigma_v r1 r2;sigma_v=x(1); r1=x(2); r2=x(3);只要更改a1和a2的值，得到a1a2ppr(1)r(2)-0.1950.950.09650.0975+0.9698i0.0975-0.9698i0.1000-0.9305-1.5950.950.03230.7975+0.5604i0.7975-0.5604i0.8179 0.3546-1.91140.950.00380.9557+0.1914i0.9557-0.1914i0.9802 0.9236（2） 产生序列data_len = 1024;trials = 100;n = 1:data_len;a1 = -1.9114;a2 = 0.95;sigma_v_2 =0.003822;v = sqrt(sigma_v_2) * randn(data_len, 1 , trials);u0 = 0 0;num = 1;den = 1 a1 a2;zi = filtic(num,den,u0);u = filter(num,den,v,zi)un=zeros(1,1024)for h=1:1024 un(h)=u(h);end（3） 的估计p = 128;r0 = xcorr(u,p,unbiased);r = r0(p+1:2 * p+1)stem(1:257,r0)（4） sar()、sbt()、sper()sar()图和程序如下：nf=1024;spr=fftshift(1/nf)*abs(fft(un,nf).2)spr=spr/max(spr);spr_unit=10*log10(spr);plot(1:1024,spr_unit); sbt()图和程序如下：m=128;r=xcorr(un,m,biased);nf=1024;bt=fftshift(fft(r,nf);bt=bt/max(bt);bt_unit=10*log10(bt);plot(1:1024,bt_unit); sper()图和程序如下：nf=1024;par=sigma_v_2./fftshift(abs(fft(1,a1,a2,nf).2);par=par/max(par);par_unit=10*log10(par);plot(1:1024,par_unit); 作业2 （1）产生序列close all; clear; clc; n=1024; pv=-10;%-10dbsnr=10;a1=10(snr/20);n=1:n;vn=randn(1,n);un=exp(j*pi*(1.4.*n-1)+exp(j*pi*(1.8.*n-0.9)+a1*vn;figure();plot(n,un); xlabel(n); ylabel(un)title(产生序列un)（2） 的估计for m=0:500 rm_e(m+1)=un(m+1:n)*un(1:n-m)/(n-m); %r(m)估计end figure();plot(1:500,rm_e(1:500); xlabel(m); ylabel(r(m)title(自相关图); （3） 估计四阶样本自相关矩阵m=4;for k=1:n-m xs(:,k)=un(k+m-1:-1:k).;endr=xs*xs/(n-m);（4） sbt()l=256;r=xcorr(un,l,biased);nf=2048;bt=fftshift(fft(r,nf);bt=bt/max(bt);bt_unit=10*log10(bt);plot(1:2048,bt_unit);smvdr()nf=2048; for n=1:nf aq=exp(-j*2*pi*(n-1)/nf*(0:m-1); pmvdr(n)=1/(aq*inv(r)*aq);endpmvdr=pmvdr/max(pmvdr);pmvdr_unit=10*log10(pmvdr);plot(1:2048,pmvdr_unit);sar()首先估计出自相关函数值r(0),r(1).r(16);利用levinson-durbin迭代算法实现ar模型的系数的求出；最后计算16阶ar模型的功率谱。p = 16;r0 = xcorr(un,p,unbiased);r = r0(p+1:2 * p+1)% stem(1:33,r0)a(1,1)=-r(2)/r(1);sigma(1)=r(1)-(abs(r(2)2)/r(1);for m=2:p k(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)/sigma(m-1); a(m,m)=k(m); for i=1:m-1 a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i); end sigma(m)=sigma(m-1)*(1-abs(k(m)2);endnf=1024;par=sigma(p)./fftshift(abs(fft(1,a(p,:),nf).2);par=par/max(par);par_unit=10*log10(par);plot(1:1024,par_unit);smuisc()u,e=svd(r);ev=diag(e);for k=1:m dec=prod(ev(k:m).(1/(m-k+1); nec=mean(ev(k:m); inv=(dec/nec)(m-k+1)*n); aic(k)=-2*log(inv)+2*(k-1)*(2*m-k+1)endamin,k=min(aic);n1=k-1;en=u(:,n1+1:m);nf=2048;for n=1:nf aq=exp(-j*2*pi*(n-1)/nf*(0:m-1); pmusic(n)=1/(aq*en*en*aq);endpmusic=pmusic/max(pmusic);pmusic_unit=10*log10(pmusic);plot(1:2048,pmusic_unit);首先根据aic准则估计出信号源个数n1，再计算出smuisc()作业3u=0.05;vt=0.0732;a1=-0.975;a2=0.95;k=100;w=0;n1=1000;n2=512;r=0;j=0;jn2=0;for k=1:k v=wgn(n1,1,10*log10(vt); x1=zeros(1,n1); x1(1)=0; x1(2)=0; for m=1:n1-2 x1(m+2)=-a1*x1(m+1)-a2*x1(m)+v(m); end; %产生n1个采样值 x=x1(n1-n2+1:n1); %取n1中后n2个值作为稳定的采样值 clear x1 v; %清楚多余变量 x(:,1)=0;0; x(:,2)=x(1);0; for n=3:1:n2 x(:,n)=x(n-1);x(n-2); end; %构造输入向量 w(:,1)=0;0; for n=1:1:n2 y(n)=w(:,n)*x(:,n); e(n)=x(n)-y(n); w(:,n+1)=w(:,n)+u*x(:,n)*conj(e(n); end; w=w+w; jn2=jn2+abs(e(n2)2; j=j+abs(e).2;end;jn2=jn2/k;j=j/k;w=w/k; %对w的估计平均w=w(:,n2+1) %输出w经512次迭代后的值eig1=1.5;eig2=0.5;m=u*(eig1+eig2)/2-u*(eig1+eig2) % m的理论值 jexn2=jn2-vt me=jexn2/vt %m的估计值jex=m*vt %剩余均方误差的理论值n=1:n2;plot(n,j);title(u=0.05,a1=-0.975,a2=0.95);xlabel(n);ylabel(j(n);仿真图如下 由图和表可以看出：u比较大时，收敛较快，但是失调较大；u比较小时，收敛较慢，在经过512次迭代之后w仍未收敛到w0，但是失调比较小。作业4n1=1000;n2=512;vt=0.0731;a1=-0.975;a2=0.95;n=100;w=0;j=0;for k=1:1:n v=wgn(n1,1,10*log10(vt); x1=zeros(1,n1); x1(1)=0; x1(2)=0; for m=1:n1-2 x1(m+2)=-a1*x1(m+1)-a2*x1(m)+v(m); end; %产生n1个采样值 x=x1(n1-n2+1:n1); %取n1中后n2个值作为稳定的采样值 clear x1 v; %清楚多余变量 x(:,1)=0;0; x(:,2)=x(1);0; for n=3:1:n2 x(:,n)=x(n-1);x(n-2); end; %构造输入向量 jmin=0.005; c=0.01; w(:,1)=0;0; k=c.*eye(2); %对k赋初值 for n=1:1:n2 an=x(:,n)*k*x(:,n)+jmin; %新息矢量的相关矩阵 gn=k*x(:,n)*inv(an); %kalman增益 an(n)=x(n)-x(:,n)*w(:,n); %新息过程 w(:,n+1)=w(:,n)+gn*an(n); %一步预测值 k=k-gn*x(:,n)*k; %更新预测误差相关矩阵 end; j=j+(abs(an).2)/n; end;n=1:1:n2;plot(n,j);title(a1=-0.975,a2=0.95,v2=0.0731);xlabel(n);ylabel(j(n); %画出kalman的学习曲线仿真图如下 和作业三的图形相比较，kalman算法估计的学习曲线收敛明显比lms算法快。作业5（1）lsw0=1.2*pi;n=25;m=18;v=wgn(n,1,-10);for n=1:n u(n)=exp(j*pi*n-j*pi)+exp(j*w0*n-j*0.79*pi)+v(n);end; %产生fblp的输入a=zeros(m,2*(n-m); %fblp输入数据矩阵b=zeros(1,2*(n-m);af=zeros(m,(n-m); %flp输入矩阵的转置共轭矩阵bf=zeros(1,(n-m);ab=zeros(m,(n-m); %blp输入矩阵的转置矩阵 bb=zeros(1,(n-m);for m=1:m for n=1:(n-m) af(m,n)=u(m-m+n); ab(m,n)=conj(u(m+n); bf(n)=u(m+n); bb(n)=conj(u(n); end;end; %构造af,ab,bf,bba=af ab;b=bf bb;y v x=svd(a); %对矩阵a进行奇异值分解theta=a*b;w=x(:,1)*x(:,1)*theta/v(1,1)2+x(:,2)*x(:,2)*theta/v(2,2)2; %解出ls权向量w=0:2*pi/999:2*pi;a=zeros(m,1000);for i=1:m a(i,:)=exp(-j*i*w);end;par=1./abs(1-w*a).2; %计算ar(m) psdplot(w,par)xlabel(w);ylabel(par(w);title(ar(m) psd for m=18 );仿真图如下 （2） rlsw0=1.2*pi;n=200;m=8;bc=0.05;j1=0;w1=0; %给lms算法赋初值l=1;c=0.004;j2=0;w2=0; %给rls算法赋初值for k=1:100 v=wgn(n,1,-10); for n=1:n u(n)=exp(j*pi*n-j*pi)+exp(j*w0*n-j*0.79*pi)+v(n); end; %产生输入数据 a=zeros(m,n); b=zeros(m-1),1); a(:,1)=zeros(m,1); for m=1:n-1 b=a(1,m);a(2,m);a(3,m);a(4,m);a(5,m);a(6,m);a(7,m); a(:,m+1)=u(m);b; end; %产生输入数据矩阵 w1=zeros(m,n); for n=1:n d(n)=w1(:,n)*a(:,n); e1(n)=u(n)-d(n); w1(:,n+1)=w1(:,n)+bc.*a(:,n)*conj(e1(n); end; %用lms算法进行迭代 j1=j1+abs(e1).2; w1=w1+w1; p=(1/c).*eye(m); k=(1/l).*p*a(:,1)./1+(1/l)*a(:,1)*p*a(:,1); b(1)=u(1); w2(:,1)=k*conj(b(1); e2(1)=u(1)-a(:,1)*w2(:,1); p=(1/l)*p-(1/l)*k*a(:,1)*p; for n=2:1:n k=1/l*p*a(:,n)/1+1/l*a(:,n)*p*a(:,n); b(n)=u(n)-a(:,n).*conj(w2(:,n-1);