function [ap_burg, en_burg] = psd_burg(xn, p);
%[ap_burg, en_burg] = psd_burg(xn, p) use the burg
%method to estimate the psd of signal xn with a p phase AR model
%ap_burg is the AR coiefficients
%en_burg is the deviation estimated

Phase=p;
N = length(xn);
efn = xn; %forward prediction error
ebn = xn; %backard prediction error
q = (sum(xn.^2))/N; %deviation
ap = zeros(Phase,1);
for p=1:Phase
    num = 0;
    don = 0;
    ap0 = ap;
    efn0 = efn;
    ebn0 = ebn;
    for n=p:(N-1)
        temp1 = efn0(n+1)*ebn0(n);
        temp2 = ((efn0(n+1))^2)+((ebn0(n+1))^2);
        num = num + temp1;
        don = don + temp2;
    end
    kp = -2*num/don;
    q  = (1 - (kp^2))*q;
    ap(p)=kp;
    for i=1:(p-1)
        ap(i)=ap0(i)+kp*ap0(p-i);
    end
    for n=(p+1):(N-1)
        efn(n+1)=efn0(n+1)+kp*ebn0(n);
        ebn(n+1)=ebn0(n)+kp*efn0(n+1);
    end
end
ap_burg = ap;
en_burg = q;