%%%%%%可能是程序寫的不好，算法的仿真結果不行，分辨力較差@@@@@@@@@
clc;
clear ;
close all;
tic
M=20;%陣元數目
N=3;%信源數目
b=5;%形成的波束個數
snap=1000;%快拍數目
C=3e8;
f0=10e6;
lamda=C/f0;
d=0.5*lamda;
% k=d/lamda;
theta0=50;
theta1=30;
theta2=20;
fs=1000;
ts=1/fs;
t=(0:snap-1)*ts;
a=[0:M-1]';%陣列矢量
u0=5;
u1=10;
u2=20;


s0=exp(j*2*pi*(f0*t+0.5*u0*t.^2));
s1=exp(j*2*pi*(f0*t+0.5*u1*t.^2));
s2=exp(j*2*pi*(f0*t+0.5*u2*t.^2));

%陣列流行矢量

a_theta0=exp(j*2*pi*d/lamda*a*sin(theta0/180*pi));
a_theta1=exp(j*2*pi*d/lamda*a*sin(theta1/180*pi));
a_theta2=exp(j*2*pi*d/lamda*a*sin(theta2/180*pi));

A=[a_theta0 a_theta1 a_theta2];
S=[s0;s1;s2];
X0=A*S;
SNR=10;
randn('state',0);
real_noise=randn(size(X0));
randn('state',3);
imag_noise=randn(size(X0));
noise0=(real_noise+j*imag_noise)/2^0.5;
noise=10^(-SNR/20)*noise0;
X=X0+noise;
Rx=X*X'/length(t);

jj=1;
for n=0:M-1
    temp=0;
    for ii=1:M-n
        temp=temp+Rx(ii,ii+n);
        
    end
  
    rt(jj)=(temp/(M-n));
    rt(2*M-jj)=conj(rt(jj));
    jj=jj+1;
end



Rt=zeros(M,M);
for iii=1:M
    for jjj=1:M
        Rt(iii,jjj)=rt(iii-jjj+M);
    end
end
Rt

[V,D]=eig(Rt);
D=abs(real((sum(D,1))));%按列相加再求實部
U=V; %所有的特征向量   

%＝＝＝＝＝＝＝尋找噪聲子空間，找到P個信源對應的大特征值并把它清空再尋找次大的
for i=1:N
    [temp,index]=max(D);%找出特征向量中的max
    D(index)=[];        %將其置空
    U(:,index)=[];      %在特征矩陣中把最大特征向量對應的列特征向量置空
end



theta=-90:0.1:90;
for iiii=1:length(theta)
    a_theta=exp(j*2*pi*d/lamda*a*sin(theta(iiii)/180*pi));
    P(iiii)=10*log10(1/sum((abs(a_theta'*U)).^2));
end

figure(1);
plot(theta,P);grid on;
toc