现代数字信号处理实验报告-2009.docVIP

现代数字信号处理实验报告 院系：电子与信息工程系 班级：互联网4班 姓名/学号：/ 姓名/学号：/ 实验一： 1.题目 Write a small MATLAB program that implements the pth order Levinson-Durbin (L-D). Run/Test the program using a AR(2) process (b0=1,a1=0, a2=0.81) and an MA(2) (bn=1,1,1) process - about 1000 samples. Use L-D with p=2 (for the AR) and 10 (for the MA). Plot the AR spectra produced in the two cases with L-D. List the direct form and the reflection coefficients in a table. Profile the LD (total number of computations for a pth order L-D). 1.1算法说明 首先由给定的已知模型产生出信号x(n) ，这里已知模型为AR(2)(参数b0=1,a1=0, a2=0.81)和MA(2)(参数bn=1,1,1)。然后AR(2) 由得到的信号x(n) 用AR(2) 模型去估计，而MA(2) 由得到的信号x(n) 用AR(10) 模型去估计。 由Wold分解定理可以得到，任何广义平衡随机过程可分解成一个完全随机的部分和一个确定的部分。它的一个推论是：如果功率谱完全是连续的，那么任何ARMA 或者AR过程，可以用一个无限阶的MA过程表示。 Kolmogorov定理也有类似结论，任何ARMA 或者MA过程可以用一个无限阶的AR过程来表示。 求解AR模型的参数只需要求解AR模型的Yule-Walker方程，常用解法，如高斯消元法需要的运算量数量级为p3而Levinson-Durbin算法的运算量数量级为p2。这是一种按阶次进行递推的算法，即首先以AR(0)和AR(1)模型参数作为初始条件，计算AR(2)模型参数，然后根据这些参数计算AR(3)模型参数，等等，一直到计算出AR(p)模型参数为止。 由k阶模型参数求k+1阶模型参数的计算公式： ， ，； 对于AR(p)模型，递推计算直到k + 1 = p 为止。 1.2 MATLAB仿真代码 1.2.1对AR模型的估计 clear all close all %由已知信号模型产生信号 N=1000; x=zeros(1,N); u=randn(1,N); for n=3:N x(n)=-0.81*x(n-2)+u(n); end % L-D算法 p=2; D=zeros(1,p); sigma=zeros(1,p+1); ak(1)=1; [X,R]=corrmtx(x,p); sigma(1)=R(1,1); ak(2)=-R(1,2)/R(1,1); sigma(2)=(R(1,1)^2-R(1,2)^2)/R(1,1); for k=2:p [X,R]=corrmtx(x,k); D(k)=R(k+1,:)*[ak(1:k) 0]; r(k+1)=D(k)/sigma(k); sigma(k+1)=(1-r(k+1)^2)*sigma(k); ak1=[ak(1:k) 0]-r(k+1)*[0 ak(k:-1:1)]; ak=ak1; end %计算功率谱密度 w=linspace(-pi,pi,N); sum=0; for n=2:p+1 sum=sum+ak(n)*exp(-(n-1)*w*1j); end Sw=sigma(p+1)./abs((1+sum).^2); an=[1 0 0.81]; for n=1:N t(n)=an(2:3)*exp(-j*w(n)*[1:2]); SxxW(n)=1/abs(1+t(n))^2; end %画图 figure1： plot(x); title(信号x(n) ); xlabel(n); ylabel(x(n)); subplot(3,1,2); plot(w/pi,SxxW); title(x(n)的功率谱密度); xlabel(PI); ylabel(Sxx(w)); subplot(3,1,3); plot(w/pi,Sw); title(AR(2)估计后的功率谱估计); xlabel(PI); ylabel(S(w)); clear all close all %由已知信号模型产生信号 N=1000; x=zeros(1,N); u=randn(1,N); for n=3:N x