clc;
clear; 
sig_num = 1000;%%信號的快拍數(shù)
p=8;%智能天線陣元數(shù)
q=2;%number of signals
x1=randint(100,1,4);%產(chǎn)生sig_sum*1的矩陣，數(shù)值在0-3之間
x2=randint(100,1,4);
y1=dmod(x1,1,0.1,4,'psk',4);
y1=y1(1:sig_num);
y2=dmod(x2,1,0.1,4,'psk',4);
y2=y2(1:sig_num);
x=[y1';y2'];
s=x;
fc=2000*10^6;%*(10^6); %hz
lemda=3*10^8/fc;
d=lemda/2;
bi=1./2;
jh=pi/180;
alpha=[5,25];%%%最小分辨率4度（20與23度只有單峰不能分辨）%%%%%%%%%%%
j=sqrt(-1);
a=[];
for i=1:q
    for n=0:(p-1)
        a(n+1,i)=exp(-j*2*pi*bi*n*sin(alpha(i)*jh));
    end
end
snr=10;
x=a*s;
for i=1:p
   xn(i,:)=awgn(x(i,:),snr,'measured');
end
R=xn*xn'./sig_num;
[U,D,V]=svd(R);
G=U(:,q+1:p);
T=G';
G1=T(:,1:q);
G2=T(:,q+1:p);
gid=zeros(1,p-q);
gid(1)=1;
c=G1'*inv(G2')*gid.';
ccc=fliplr(c.');
cc=[1,ccc];
ra=roots([cc]);
w=angle(ra);
doa=asin(w*lemda./(d*pi*2))./jh