function [a, e] = AR2(x, p)
% AR自相關(guān)算法譜估計，輸入采樣數(shù)據(jù)x，預(yù)測階數(shù)p，輸出模型傳輸函數(shù)系數(shù)a
% 返回H(z)=e/(1+A(z))中的A，e是sigma^2
N = length(x);
A = zeros(p, N+p);
for n=1:p
    for m=1:N+p
        if m<=n || m-n>=N+1
            A(n, m) = 0;
        else 
            A(n, m) = x(m-n);  %A'
        end
    end
end
xf = zeros(N+p, 1);
xf(1:N) = x';
a = -1*pinv(A')*xf;
e = (xf'*xf+a'*A*xf)/(N+p);
%test
% clear;
% close all;
% a = [1, -1.3817, 1.5632, -0.8843, 0.4096];
% b = [1, 0.3544, 0.3508, 0.1736, 0.2401];
% x = xx(a, b, 200, 256);
% [coa, coe] = AR2(x, 4);
% coa = [1;coa];
% cob = sqrt(coe);
% [w,P] = spectrum(coa, cob, 128);
% plot(w, P);