1、现代数字信号处理及其应用仿真作业第三章仿真作业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+
2、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);s
3、tem(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
4、*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.
5、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/BTma
6、x;%归一化幅度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); s
7、ignal3=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);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
8、sigma(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;clear all;close all;N=1000;noise=(randn(1,N)+1i*randn(1,N)/sqrt(2);s
9、ignal1=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);for k=1:N-M xs(:,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); for m=1:M co(m:m+M-1)=co(m:m+M-1)+Gr(M:-1:1,m)
10、; end z=roots(co); ph=angle(z)/(2*pi); z1=abs(z); z2=abs(z1-1);estallv,estain = sort(z2,ascend);所以单Root-MUSIC算法中最近单位圆的两个根为0.0033 - 0.9977i 0.5856 + 0.8074i上述根的相位对应的归一化频率为0.2495 -0.1501（2）clc;clear all;close all;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);
11、signal2=exp(-j*0.3*pi*(0:N-1)+j*2*pi*rand);un=signal1+signal2+noise;M=8;xs=zeros(M,N-M);for k=1:N-M xs(:,k)=un(k+M-1:-1:k).;end R=(xs*xs)/(N-M); U,E,V=svd(R); G=U(:,3:M); f=-0.5:1/999:0.5;for ff=1:length(f) w=f(ff)*2*pi; for l=1:M aw(l)=exp(-j*w*(l-1); %计算a(w) end WW=aw*G*G*aw; Pmusic(ff)=abs(1./WW)
12、; %谱扫描函数endPmusic=10*log10(Pmusic);f=-0.5:1/999:0.5;plot(f,Pmusic); hold on %绘图输出xlabel(/2)ylabel(归一化Music谱/dB);3.21clc;clear all;close all;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
13、+signal2+noise;M=8;xs=zeros(M,N-M);for k=1:N-M xs(:,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=di
14、ag(E);ph=angle(z)/(2*pi);z1=abs(z); z2=abs(z1-1);estallv,estain = sort(z2,ascend);所以单次ESPRIT算法中最近单位圆的两个特征值为 0.5845 + 0.8117i 0.0036 - 0.9994i上述特征值的相位对应的归一化频率为 -0.1507 0.2494第四章仿真作业4.18%（1）产生N=512点的样本序列data_len = 512; %样本序列的长度trials = 100; %随机试验的次数A=zeros(data_len,2);EA=zeros(data_len,1);B=zeros(data
15、_len,2);EB=zeros(data_len,1);for m = 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 = 0 0;num = 1;den = 1 a1 a2;Zi = filtic(num, den, u0); %滤波器的初始条件u = filter(num, den, v, Zi); %产生样本序列u(n)%（2）用LMS滤波器来估计w1和w2mu1 = 0.05;mu2 = 0.005;w1 = zero
16、s(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迭代过程for n =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(
17、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(tria
18、ls,1);for m=1:trials;rm=xcorr(u(:,:,m),biased);R=rm(512),rm(513);rm(511),rm(512);p=rm(511);rm(510);wopt(:,m)=Rp;v,d=eig(R);Jmin(m)=rm(512)-p*wopt(:,m);sum_eig(m)=d(1,1)+d(2,2);endsJmin=sum(Jmin)/trials;e1_100trials_ave=sum(e1)/trials;e2_100trials_ave=sum(e2)/trials;Jex1=e1_100trials_ave-sJmin;Jex2=e
19、2_100trials_ave-sJmin;sum_eig_100trials=sum(sum_eig)/100;Jexfin=mu1*sJmin*(sum_eig_100trials/(2-mu1*sum_eig_100trials);Jexfin2=mu2*sJmin*(sum_eig_100trials/(2-mu2*sum_eig_100trials);M1=Jexfin/sJminM2=Jexfin2/sJminfigure(1);plot(A/trials);hold on;plot(conj(w1);xlabel(迭代次数);ylabel(权向量);title(步长为0.05权向
20、量收敛曲线);figure(2);plot(B/trials);hold on;plot(conj(w2);xlabel(迭代次数);ylabel(权向量);title(步长为0.005权向量收敛曲线);figure(3);plot(EA/trials,*);hold on;plot(EB/trials,-);xlabel(迭代次数);ylabel(均方误差);title(步长分别为0.05和0.005学习曲线);仿真结果：剩余均方误差和失调参数结果分别为：M1 =0.0502M2 =0.0048第五章仿真作业5.10解：（1）求解M=2时的相关矩阵。根据u（n）表达式可知a1=-0.99,2
21、=0.93627,构造Yule-Walker方程 由于u(n)为是信号，故根据自相关函数的共轭对称性，有r(1)=r(-1),求解方程可得 故滤波器权系数个数为M=2时的相关系数为 （2）根据（1）中同样的方法，可得滤波器权系数个数为M=3时的相关系数为 （3）对相关矩阵进行特征分解 R2=47.0487 46.5783;46.5783 47.0487;R3=47.0487 46.5783 46.1125;46.5783 47.0487 46.5783;46.1125 46.5783 47.0487;eig_value2=eig(R2)eig_value3=eig(R3)eig_spread2
22、=max(eig_value2)/min(eig_value2)eig_spread3=max(eig_value3)/min(eig_value3)（4）根据LMS算法均方误差收敛特性，可知步长因子应满足 当M=2时，步长因子应满足，M=3时，有，所以题目中给的不符合要求，故在试验中M=2时，令。clc;clear all;M=2;%抽头个数data_len=10000;N=500;n=1:data_len;a1=-0.99;sigma_v=0.936271;v=sqrt(sigma_v)*randn(data_len,1,N);u0=0;num=1;den=1 a1;Zi=filtic(n
23、um ,den,u0);u=filter(num,den,v,Zi);u_len1=0.0006;w1=zeros(M,data_len);e1=zeros(data_len,1);d1=zeros(data_len,1);ww1=zeros(2*N,data_len);%求学习曲线for i=1:Nfor n=M+1:data_len-1 w1(:,n+1)=w1(:,n)+u_len1*u(n-1:-1:n-M,:,i)*conj(e1(n); d1(n+1)=w1(:,n+1)*u(n:-1:n-M+1,:,i); e1(n+1)=u(n+1,:,i)-d1(n+1); ee1(n)=e
24、1(n)*e1(n);endwe1(i,:)=ee1;endEwe1=mean(we1);plot(Ewe1)hold onxlabel(迭代次数);ylabel(均方误差);title(步长为0.001的500次独立实验的学习曲线);第六章仿真作业6.13M=8;%阶数N=1000;noise=(randn(1,N)+j*randn(1,N)/sqrt(2);f1=0.1;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
25、*(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;A=zeros(M,N-M+1);for n=1:N-M+1 A(:,n)=un(M+n-1:-1:n);endU,S,V=svd(A);invphi=V*inv(S*S)*V;P=1024;f=linspace(-0.5,0.5,P);omega=2*pi*f;a=zeros(M,P);for k=1:P for m=1:M a(m,k)=exp(-j*omega(k)*
26、(m-1); endendPmvdr=zeros(size(omega);for k=1:P Pmvdr(k)=1/(a(:,k)*invphi*a(:,k);endPmvdr=abs(Pmvdr/max(Pmvdr);Pmvdr=10*log10(Pmvdr);f=(-(P/2)+1:(P/2)/P;plot(f,Pmvdr);xlabel(f);ylabel(归一化功率谱/dB);title(M=8时MVDR谱);6.15clc;clear all;M=2;%抽头个数data_len=1000;N=500;a1=0.99;sigma_v=0.995;W=zeros(data_len,2);
27、for i=1:Nv=sqrt(sigma_v)*randn(data_len,1);num=1;den=1 a1;u0=zeros(length(den) -1,1);Zi=filtic(num ,den,u0);u=filter(num,den,v,Zi);n0=1;b=u(n0+1:data_len);L=length(b);u1=zeros(M-1,1).,u.;A=zeros(M,L);for k=1:L A(:,k)=u1(M-1+k:-1:k);end delta=0.004;lambda=0.98;w=zeros(M,L+1);epsilon=zeros(L,1);P1=eye
28、(M)/delta; for k=1:L-1 PIn=P1*A(:,k+1); denok=lambda+A(:,k+1).*PIn; kn=PIn./denok; epsilon(k+1)=b(k+1)-w(:,k).*A(:,k+1); w(:,k+1)=w(:,k)+kn.*conj(epsilon(k+1); P1=P1/lambda-kn*A(:,k+1).*P1/lambda; epsilon2(k+1)=epsilon(k+1)*epsilon(k+1);endMSE(i,:)=epsilon2;%也可以求和在求平均W =W + conj(w);endfigure(1);MSE=
29、mean(MSE);plot(MSE);xlabel(迭代次数);ylabel(MSE);title(500次重复独立实验的学习曲线);figure(2);plot(W/500);hold onplot(conj(w);xlabel(迭代次数);ylabel(权值);title(权向量收敛曲线);第七章仿真作业7.14Clear；close all；N=3000;sigma_v=0.0332;v=sqrt(sigma_v)*randn(1,N);%根据给定的AR模型产生u（n）序列a1=1.6;a2=-1.46;a3=0.616;a4=-0.1525;u(1:4)=0;for i=1:(N-4) u(i+4)=a1*u(i+
