clc
clear all
signal_number=2;
n1=1:0.1:1024;
signal_1=100*cos(2*pi*200/4000*n1);
signal_2=100*cos(2*pi*200/4000*n1+pi/2);
angle_1=pi*45/180;
angle_2=pi*100/180;

%采樣
sample_number=1024;%采樣點數 
n2=1:sample_number; 
signal_1=signal_1(n2);
signal_2=signal_2(n2);
S=[signal_1;signal_2];

%確定導向矢量以及接收信號
array_number=4;
t=0.34;%距離除以波長
i=sqrt(-1); 
j=i;
delay_1=[j*2*pi*t*cos(angle_1);
       j*2*pi*t*cos(angle_1-pi/2);
       j*2*pi*t*cos(angle_1-pi);
       j*2*pi*t*cos(angle_1-3*pi/2)];
delay_2=[j*2*pi*t*cos(angle_2);
       j*2*pi*t*cos(angle_2-pi/2);
       j*2*pi*t*cos(angle_2-pi);
       j*2*pi*t*cos(angle_2-3*pi/2)]; 
delay=[delay_1,delay_2];   
A=exp(delay); %導向矢量
X=A*S; %接收信號 

%確定協方差矩陣以及對協方差矩陣進行分解 
cov=X*X'/sample_number;%協方差矩陣 
[V,D]=eig(cov); 
signal_space=V(:,array_number-signal_number+1:array_number);%確定信號子空間

angle=1:0.1:360;
%設置最后要求的方向矢量函數
A_value=[j*2*pi*t*cos(angle*pi/180);
       j*2*pi*t*cos(angle*pi/180-pi/2);
       j*2*pi*t*cos(angle*pi/180-pi);
       j*2*pi*t*cos(angle*pi/180-3*pi/2)];
A_value=exp(A_value);  

%確定適應度函數
for n=1:length(angle);
PA=A_value(:,n)*inv(A_value(:,n)'*A_value(:,n))*A_value(:,n)';
tr_matrix=PA*signal_space*signal_space';
[M N]=size(tr_matrix);
fitness(n)=0;
for n3=1:M
    fitness(n)=fitness(n)+tr_matrix(n3,n3);
end

fitness(n)=abs(fitness(n));
end