clc
clear
N=128;  %采樣點%
n=1:N;
M=110;
QE=110;
PE=100;
Xn=(20.^0.5)*sin(2*pi*0.2*n)+(2.^0.5)*sin(2*pi*0.213*n)+randn(size(n));
Rx1=xcorr(Xn,'unbiased');
for i=1:M;
    for j=0:PE-1;
        R1(i,j+1)=Rx1(QE+i-j+1);
    end
end
r1=Rx1(QE+1:QE+M);
r=conj(r1)';
R=[r,R1];
%求SVD%
[U,S,V] = svd (R);
%求階數p%
for k=1:M;
    a=S(k,k)/S(1,1);
    if(a<0.01)
        break
    end
end
p=k;
%求S（p）矩陣%
Sp=zeros(p+1,p+1);
for j=1:p;
    for i=1:(PE+1-p);
    Sp=(S(j,j).^2)*V(i:i+p,j)*(conj(V(i:i+p,j))')+Sp;
    end
end
%求參數%
Spni=inv(Sp);
for i=1:p;
    x(1,i)=Spni(i+1,1)/Spni(1,1);
end
%求功率譜%
freqz(1,[1,x],1024,1);