clear;clc;
N_x=8192;
N=6;
l=10;
d=0.5*l;
M=2;
w=[24*pi*(10.^7) 12*pi*(10.^7)]';
R=6;
x1=6;x2=4;x3=-4;x4=-6;x5=-4;x6=4;
y1=0;y2=4.47;y3=4.47;y4=0;y5=-4.47;y6=-4.47;
fy1=10;fw1=10;fy2=40;fw2=40;
B=[exp(-j*2*pi*(x1*cos(fy1/180*pi)*cos(fw1/180*pi)+y1*sin(fy1/180*pi)*cos(fw1/180*pi))/l)  exp(-j*2*pi*(x2*cos(fy1/180*pi)*cos(fw1/180*pi)+y2*sin(fy1/180*pi)*cos(fw1/180*pi))/l)  exp(-j*2*pi*(x3*cos(fy1/180*pi)*cos(fw1/180*pi)+y3*sin(fy1/180*pi)*cos(fw1/180*pi))/l) exp(-j*2*pi*(x4*cos(fy1/180*pi)*cos(fw1/180*pi)+y4*sin(fy1/180*pi)*cos(fw1/180*pi))/l) exp(-j*2*pi*(x5*cos(fy1/180*pi)*cos(fw1/180*pi)+y5*sin(fy1/180*pi)*cos(fw1/180*pi))/l) exp(-j*2*pi*(x6*cos(fy1/180*pi)*cos(fw1/180*pi)+y6*sin(fy1/180*pi)*cos(fw1/180*pi))/l);
   exp(-j*2*pi*(x1*cos(fy2/180*pi)*cos(fw2/180*pi)+y1*sin(fy2/180*pi)*cos(fw2/180*pi))/l)  exp(-j*2*pi*(x2*cos(fy2/180*pi)*cos(fw2/180*pi)+y2*sin(fy2/180*pi)*cos(fw2/180*pi))/l)  exp(-j*2*pi*(x3*cos(fy2/180*pi)*cos(fw2/180*pi)+y3*sin(fy2/180*pi)*cos(fw2/180*pi))/l) exp(-j*2*pi*(x4*cos(fy2/180*pi)*cos(fw2/180*pi)+y4*sin(fy2/180*pi)*cos(fw2/180*pi))/l) exp(-j*2*pi*(x5*cos(fy2/180*pi)*cos(fw2/180*pi)+y5*sin(fy2/180*pi)*cos(fw2/180*pi))/l) exp(-j*2*pi*(x6*cos(fy2/180*pi)*cos(fw2/180*pi)+y6*sin(fy2/180*pi)*cos(fw2/180*pi))/l);]'
   
xxx=2*exp(j*w*[0:N_x-1]);
x=B*xxx;
R=x*x';
[V D]=eig(R);
[lambda,index]=sort((diag(D)));
UU=V(:,index(1:2));
%ULA估計(jì)角度變化的范圍和頻率選擇
 fy=0:90;
 fw=0:90;
 for i=1:length(fy)
 for jj=1:length(fw)
  
   AA=[exp(-j*2*pi*(x1*cos(fy(i)/180*pi)*cos(fw(jj)/180*pi)+y1*sin(fy(i)/180*pi)*cos(fw(jj)/180*pi))/l)  exp(-j*2*pi*(x2*cos(fy(i)/180*pi)*cos(fw(jj)/180*pi)+y2*sin(fy(i)/180*pi)*cos(fw(jj)/180*pi))/l)  exp(-j*2*pi*(x3*cos(fy(i)/180*pi)*cos(fw(jj)/180*pi)+y3*sin(fy(i)/180*pi)*cos(fw(jj)/180*pi))/l) exp(-j*2*pi*(x4*cos(fy(i)/180*pi)*cos(fw(jj)/180*pi)+y4*sin(fy(i)/180*pi)*cos(fw(jj)/180*pi))/l) exp(-j*2*pi*(x5*cos(fy(i)/180*pi)*cos(fw(jj)/180*pi)+y5*sin(fy(i)/180*pi)*cos(fw(jj)/180*pi))/l) exp(-j*2*pi*(x6*cos(fy(i)/180*pi)*cos(fw(jj)/180*pi)+y6*sin(fy(i)/180*pi)*cos(fw(jj)/180*pi))/l);];
   WW=AA*UU*UU'*AA';
   Pmusic(i,jj)=abs(1./WW);%角譜
 % PCap(i) = abs(1/(AA'*inv(R)*AA));
 
end
end
mesh(fy,fw,Pmusic);
grid on
xlabel('rad')