求解的题目如下:
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 编辑]