大连理工大学信号第11章作业编程题答案

4、利用MATLAB编程实现：
已知AR(4)模型为：
，() 2.7607(1) 3.8106(2) 2.6535(3)0.9238(4)()−−+−−−+−= x n x n x n x n x n w n
其中()
w n是零均值、单位方差的高斯白噪声。

w n；
（1）产生长度为1024点的高斯分布白噪声()
（2）让该白噪声通过这个AR模型；
x n；
（3）得到输出信号()
x n的AR模型的参数；
（4）再估计()
（5）比较估计的结果与原来给定的AR参数。

程序清单：
%******利用AR模型得到观测数据******%
clc
clear %清除数据
P=4; %AR模型阶数
N=1024; %生成数据长度
w=wgn(N,1,0); %方差为1高斯白噪声
x=zeros(N+P,1); %观测数据向量
a=[2.7607, -3.8106, 2.6535, -0.9238]; %权向量
for n=P+1:1:N+P %求解观测向量X
x(n)=a(1:P)*x(n-1:-1:n-P)+w(n-P);
end
x=x(P+1:end); %得到观测向量X
x=x-mean(x); %观测向量去均值
r=zeros(P+1,1); %求自相关系数
for m=0:1:P
r(m+1)=x(1:end-m)'*x(m+1:end);
end
r=r./N; %有偏的求平均
%****利用高斯消去法求解AR模型系数****%
R=zeros(P,P);
for m=1:P-1
R(m,m:end)=r(1:P+1-m);
R(m:end,m)=r(1:P+1-m);
end
R(P,P)=r(1);
b=r(2:end)*-1;
a=R\b;
a=[1;a]; %求解AR模型的系数
sigma=r'*a; %得到白噪声方差
%****利用levinson算法求解AR模型系数*****%
A(1,1)=-r(2)/r(1);
for m=2:1:P %首先求解A(m,m)
temp1=A(m-1,1:m-1)*r(2:m);
temp2=A(m-1,1:m-1)*r(m:-1:2);
A(m,m)=-(r(m+1)+temp2)/(r(1)+temp1);
for n=1:1:m-1
A(m,n)=A(m-1,n)+A(m,m)*A(m-1,m-n);
end %迭代求解其他的系数
end
%******利用Burg算法求解AR模型系数*****%
Error_f=zeros(P,N);
Error_b=zeros(P,N);
Error_f(1,:)=x; %初始化前向预测误差
Error_b(1,:)=x; %初始化后向预测误差
B=zeros(P,P);
temp1=Error_f(1,2:end)*Error_b(1,1:end-1)';
temp2=sum(Error_f(1,2:end).^2+Error_b(1,1:end-1).^2);
B(1,1)=-2*temp1/temp2;
for m=2:1:P %前后向预测误差的迭代求解
Error_f(m,m:end)=Error_f(m-1,m:end)+B(m-1,m-1)*Error_b(m-1,m-1:end-1); Error_b(m,m:end)=Error_f(m-1,m:end)*B(m-1,m-1)+Error_b(m-1,m-1:end-1);
temp1=Error_f(m,m+1:end)*Error_b(m,m:end-1)';
temp2=sum(Error_f(m,m+1:end).^2+Error_b(m,m:end-1).^2);
B(m,m)=-2*temp1/temp2; %求解系数B(m,m)
for n=1:1:m-1 %利用levinson求解其他系数
B(m,n)=B(m-1,n)+B(m,m)*B(m-1,m-n);
end
end
结果与分析：
分析：以上是两次单独运行程序得到的结果，对比两次结果容易发现：利用高斯消去法和Levinson算法求解得到的AR模型的参数是相等的，而利用Burg算法求解得到的AR模型参数较前述的两种算法精度更高，更加接近AR模型的实际参数，而且通过对比容易发现，利用高斯消去法和Levinson算法两次求解的AR模型参数相差较大，而Burg算法的参数波动较小，为了更好的说明算法求解的不同，将500次的求解叠加平均，结果如下所列：
分析：以上结果是运行500次后叠加平均的结果，从表中容易发现，利用高斯消去法和Levinson算法求解得到的AR模型参数与实际参数相对误差大，而利用Burg算法得到的AR 模型参数与实际参数相对误差小，这就说明了利用Burg算法求解AR模型参数的优越性。