數(shù)字信號處理仿真作業(yè)

1、作業(yè)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.1914i

2、0.9557-0.1914i0.9802 0.9236（2） 產(chǎn)生序列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） 的估

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=1

4、0*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); 作業(yè)2 （1）產(chǎn)生序列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

5、-0.9)+a1*vn;figure();plot(n,un); xlabel(n); ylabel(un)title(產(chǎn)生序列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,

6、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()首先估計出自相關函數(shù)值r(0),r(1).r(16);利用levinso

7、n-durbin迭代算法實現(xiàn)ar模型的系數(shù)的求出；最后計算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)=sig

8、ma(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=

9、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);首先根據(jù)aic準則估計出信號源個數(shù)n1，再計算出smuisc()作業(yè)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

10、;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; %產(chǎn)生n1個采樣值 x=x1(n1-n2+1:n1); %取n1中后n2個值作為穩(wěn)定的采樣值 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

11、)*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經(jīng)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(

12、u=0.05,a1=-0.975,a2=0.95);xlabel(n);ylabel(j(n);仿真圖如下 由圖和表可以看出：u比較大時，收斂較快，但是失調較大；u比較小時，收斂較慢，在經(jīng)過512次迭代之后w仍未收斂到w0，但是失調比較小。作業(yè)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); e

13、nd; %產(chǎn)生n1個采樣值 x=x1(n1-n2+1:n1); %取n1中后n2個值作為穩(wěn)定的采樣值 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(:,

14、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的學習曲線仿真圖如下 和作業(yè)三的圖形相比較，kalman算法估計的學習曲線收斂明顯比lms算法快。作業(yè)5（1）lsw0=1.2*pi;n=25;m=18;v=wgn(n,1,-10);for n=1:n u(n)=exp(j*p

15、i*n-j*pi)+exp(j*w0*n-j*0.79*pi)+v(n);end; %產(chǎn)生fblp的輸入a=zeros(m,2*(n-m); %fblp輸入數(shù)據(jù)矩陣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

16、; %構造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 );仿真圖

17、如下 （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; %產(chǎn)生輸入數(shù)據(jù) 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

18、(7,m); a(:,m+1)=u(m);b; end; %產(chǎn)生輸入數(shù)據(jù)矩陣 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);