《现代数字信号处理》仿真题(3-4章)

3.17(1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FFT法与直接法求自相关函数的比较%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;f1=0.15;f2=0.17;f3=0.26;n=1:32;d=1;%白噪声的方差取1N=length(n);A1=sqrt((10^3)*2);A2=A1;A3=sqrt((10^2.7)*2);s=A1*cos(2*pi*f1*n+unifrnd(0,2*pi))+A2*cos(2*pi*f2*n+unifrnd(0,2*pi))+A3*cos(2*pi*f3*n+unifrnd(0,2*pi))+normrnd(0,1,1,N);%直接法求自相关函数r1=zeros(1,N);k=1;fork1=0:N-1fork2=k1+1:Nr1(k)=r1(k)+s(k2)*(s(k2-k1))';endr1(k)=r1(k)/N;k=k+1;end%FFT法求自相关函数s=[s,zeros(1,N)];r2=fft(s,2*N);r2=abs(r2).^2/N;r3=ifft(r2);比较结果：两种方法求出的自相关函数相同。(2）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用BT法和周期图法分别进行功率谱估计%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用FFT法求出自相关函数s1=[s,zeros(1,N)];r2=fft(s1,2*N);r3=abs(r2).^2/N;r=ifft(r3);r4=rot90(rot90(r));%用BT法进行功率谱估计M=64;f=-0.5:0.002:0.5;S1=r(1)*ones(1,length(f));fork=2:MS1=S1+r(k)*exp(-j*2*pi*f*(k-1))+r4(k-1)*exp(j*2*pi*f*(k-1));endsubplot(2,1,1),plot(f,S1,'red'),grid;title('BT法功率谱估计得到的曲线');%用周期图法进行功率谱估计U=zeros(1,length(f));fork=1:NU=U+s(k)*exp(-j*2*pi*f*k);endS2=(abs(U).^2)/N;subplot(2,1,2),plot(f,S2,'black'),grid;title('周期图法功率谱估计得到的曲线');仿真结果如下：(3)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用AR模型法进行功率谱估计%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求自相关函数r1=zeros(1,N);r2=fft(s,2*N);r3=abs(r2).^2/N;r=ifft(r3);a=zeros(1,p);a1=zeros(1,p);%用Levinson-Drubin迭代算法求解AR模型系数form=1:pifm==1a(1)=-r(2)/r(1);D=r(1)-abs(r(2))^2/r(1);elsek=0;fori=1:m-1k=k+a(i)*r(m+1-i);endk=-(k+r(m+1))/D;a1=a;a(m)=k;fori=m-1:-1:1a(i)=a1(i)+k*a1(m-i);endD=D*(1-abs(k)^2);endendf=-0.5:0.002:0.5;S=zeros(1,length(f));S1=zeros(1,length(f));fori=1:pS1=S1+a(i)*exp(-j*2*pi*f*i);endS1=abs(S1+1).^2;S=D./S1;plot(f,S),grid;title('16阶AR模型法功率谱估计曲线');仿真结果如下：3.20%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用ROOT-MUSIC法进行功率谱估计%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;M=8;n=1:1000;N=length(n);s=exp(j*0.5*pi*n+j*unifrnd(0,2*pi))+exp(-j*0.3*pi*n+j*unifrnd(0,2*pi))+normrnd(0,1,1,N);r1=fft(s,2*N);%求自相关函数r2=(abs(r1)).^2/N;r=ifft(r2,2*N);R=zeros(M);form=1:M%求自相关矩阵forn=1:Mifn-m0R(m,n)=conj(r(1,m-n+1));elseR(m,n)=r(1,n-m+1);endendend[U,V]=eig(R);G=U(:,1:6);%提取噪声子空间向量V(1,1)symsza1=z.^([0:7]);a2=(z^(-1)).^([0:7]');f=z^7*a1*G*G'*a2;a=sym2poly(f);Y=roots(a);Y1=abs(Y);I=zeros(1,4);fori=1:4%寻找离单位圆最近的四个根（信号源有两个），因为根以共轭形式出现，所以有四个[P2,I(i)]=min(abs(Y1-1));Y1(I(i))=inf;endfori=1:4%求出这四个根的相位，即信号频率的估计w(i)=angle(Y(I(i)))/(2*pi);end运行结果：w=0.25020.2502-0.1497-0.14973.21%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用ESPIRT法进行功率谱估计%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;M=8;n=1:1000;N=length(n);s=exp(j*0.5*pi*n+j*unifrnd(0,2*pi))+exp(-j*0.3*pi*n+j*unifrnd(0,2*pi))+normrnd(0,1,1,N);r1=fft(s,2*N);%求自相关函数r2=(abs(r1)).^2/N;r=ifft(r2,2*N);R=zeros(M+1);form=1:M+1%求自相关矩阵forn=1:M+1ifn-m0R(m,n)=conj(r(1,m-n+1));elseR(m,n)=r(1,n-m+1);endendendRxx=R(1:M,1:M);%Rxx为M-1维的自相关矩阵Rxy=R(2:M+1,1:M);%Rxy为M-1维的互相关矩阵Z=[zeros(1,M-1);eye(M-1)]';Z=[Z;zeros(1,M)];Cxx=Rxx-min(eig(Rxx))*eye(M);Cxy=Rxy-min(eig(Rxx))*Z;a=eig(Cxx,Cxy);a1=abs(abs(a)-1);I=zeros(1,2);fori=1:2%寻找离单位圆最近的两个根（信号源有两个）[P,I(i)]=min(a1);a1(I(i))=inf;endfori=1:2%求出这两个根的相位，即信号频率的估计w(i)=angle(a(I(i)))/(2*pi);end运行结果：w=0.2498-0.15104.18%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%LMS算法进行二阶线性预测%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;a=[1-0.9750.95];%滤波器系数mu=0.005;%初始步长N=512;%样本序列长度M=5;form=1:100;%滤波器独立进行100次实验v=0.0731.*randn(1,N);%滤波器输入为零均值，方差为0.0731的加性白噪声u=filter(1,a,v);%u为输出信号w=zeros(2,N);forn=3:N%LMS算法u1=u(n-1:-1:n-2);y(n)=u1*w(:,n-1);e(m,n)=u(n)-y(n);w(:,n)=w(:,n-1)+mu*u1'*e(m,n);endende=sum(e.^2)/100;plot(e);title('学习曲线');figureplot(w(1,:),'k');holdonplot(w(2,:),'r');title('权向量的变化情况')图1.步长为0.005时的学习曲线和权向量更新情况图2.步长为0.05时的学习曲线和权向量更新情况以上图1和图2比较可知，当步长较大时，收敛速度较快，但失调和剩余均方误差较大，从而使稳态性能变差；而步长较小时，收敛速度虽然较慢，但失调和剩余均方误差减小，从而可以改善稳态性能；