求解的题目如下：

Write a small MATLAB program that implements the pthorder 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 L-D (total number of computations for a pthorder L-D).

列文森-杜宾算法我已经实现了，用AR(2) process (b0=1,a1=0, a2=0.81) 来验证程序也是正确的，现在的关键问题是对于MA(2)过程的功率谱问题。我做的方法如下：

clear all; close all; clc;

N = 1000;

p_AR = 2;

p_MA = 10;

%平稳随机信号x2是由高斯白噪声通过一个线性定常系统产生，信号模型

%为：x(n)=u(n)+u(n-1)+u(n-2)，u(n)为高斯白噪声

a21 = 1; b20 = 1; b21 = 1; b22 = 1;

a2 = [a21];

b2 = [b20,b21,b22];

x2 = filter(b2,a2,u); % x(n)=u(n)+u(n-1)+u(n-2)

%---------------------计算信号的自相关矩阵-----------------

R2 = xcorr(x2,'coeff') ;

R2 = R2(N:end);

%-------------------调用Levinson-Durbin算法---------------

[A2,E2]=Levinson_Durbin_Algo(R2,p_MA);

%------利用Levinson-Durbin算法得到的模型参数计算功率谱------

M=1000;%在0~2*pi的范围内等角距的取1000个点，用于计算傅里叶变换

for i=1:M

sum2=0;

for k=2:p_MA+1

sum2=sum2+A2(k)*exp(-j*2*pi/M*(k-1)*(i-1));

end

S2(i)=E2/(abs(1+sum2)*abs(1+sum2));

end

%----------------画出功率谱的图形---------------------------------------

W = linspace(-pi,pi,M);

figure;

subplot(2,1,2);

plot(W,(abs(S2)))

xlabel('\omega');

title('Case2:x(n)=u(n)+u(n-1)+u(n-2)下用L-D算法求得的功率谱');

最后得到的结果与实际结果相差很大，请版上各位高手指点迷津！

[本帖最后由 mingcheng 于 2009-12-5 20:53 编辑]