AR建模仿真自相关算法

%AR（3）产生数据，用3阶的LD算法计算，输出估计的功率谱
clear;
p=100;          %AR计算的阶数

u=wgn（100011）;   %u（n）是一个均值为零，功率为一的白噪声
x（1）=u（1）;
x（2）=u（2）;
x（3）=u（3）;
for kk=4:1000
   x（kk）=14/24*x（kk-1）+9/24*x（kk-2）-1/24*x（kk-3）+u（kk）;
end           %产生指定的AR（3）模型的x（n）
figure 
pburg（x3）
R = zeros（1 p+1）;
   sample_num = length（x）;
   for n = 1:p+2
       R（n） = mean（ x（1:sample_num-n+1） .* x（n:sample_num） ）;%自相关函数
   end
%for  n= 2:p+2
   % R（n）=R（n）/R（1）;
%end
    ssigma = zeros（1 p）;   % 噪声方差矩阵
    D = zeros（1 p）;        % Dk矩阵
    gama = zeros（1 p）;     % gama矩阵
    a = zeros（p p）;       
   
    a0 = 1;
    D0 = a0 * R（1+1）;
    ssigma0 = R（1）;
    gama（1） = D0 / （ssigma0 + eps）;
    a（11） = -gama（1）;
    
    for k = 1:p
        %公式中的 R（k） is R（k+1） in this program
        ssigma（k） = R（1） + sum（ a（k 1:k） .* R（2:k+1） ）;
        D（k） = a0*R（k+2） + sum（ a（k 1:k） .* fliplr（R（2 : k+1