% 離散隨機信號處理 Project I
% 功率譜估計
% 馮偉  2005310367  fengw@wireless.mdc.tsinghua.edu.cn
% 功能: 自相關法進行譜估計   2007年5月18日

function S=PSD_zxg(x,p,M)     
N=length(x);
pi=3.1415926;
x_Np = zeros(1,N+p)';
for ii=1:N
    x_Np(ii) = conj(x(ii));
end
A_H_Np = zeros(p,N+p);
S0 = zeros(1,M+1);
for ii=1:p
    for jj=1:N+p
        if (jj-ii) > 0 && (jj-ii) < N+1
            A_H_Np(ii,jj)=x(jj-ii);
        else
            A_H_Np(ii,jj) = 0;
        end
    end
end

% AR模型參數(shù)估計
a=-(A_H_Np*A_H_Np')^-1*A_H_Np*x_Np;
Pv=1/(N+p)*(x_Np'*x_Np+a'*A_H_Np*x_Np);

% 功率譜估計
for n=1:M+1
    k=1:p;
    S0(n)=10*log10(Pv/abs(1+a'*exp(-j*pi/M*(n-1)*k).')^2);
end
S = S0;