clear all
n=1:128;
N1=128;
m=30;
%p=2;
N=50;
%設(shè)置觀測數(shù)據(jù)
x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(1,N1);
%建立自相關(guān)矩陣
Rxx=zeros(m);
for k=1:N;
        for i=1:m;
             xm(i)=x(i+k-1);
        end
    rxx=xm'*xm;
    Rxx=Rxx+rxx;
end
Rx=Rxx/N;
%做特征值分解
[V,D]=eig(Rx);
p=4;
%對(duì)特征值和特征向量進(jìn)行處理
for i=1:m;
    d(i)=D(m+1-i,m+1-i);
end
V1=zeros(m);
for i=1:m;
    V1(1:m,m+1-i)=V(1:m,i);
end
%取信號(hào)子空間
% l=m;
% for i=1:m-1;
%    dd=d(i)-d(i+1);
%    if dd>l,
%       l=i;
%        break
%    end
% end
s=zeros(m,p);
for i=1:p,
   s(1:m,i)=V1(1:m,i);
end
%用信號(hào)子空間公式求功率譜
df=0.001;
I=eye(m);
f=0.001:0.001:0.5;
for i=1:length(f);
    wi=2*pi*df*i;
    for t=1:m;
        aw(t)=exp(-j*(t-1)*wi);
    end
    p(i)=1/(aw*(I-s*s')*aw');
end
pw=10*log10(abs(p));
plot(f,pw);