clc;
clear;
sig_num = 2000;%%信號的快拍數
fc=2000; %mhz
p=8;%智能天線陣元數
q=2;%number of signals
snr=10;
rou=0;
alp=[10,30].';
alpha=alp*pi./180;
rou=0;
x1=randint(sig_num,1,4);%產生sig_sum*1的矩陣，數值在0-3之間
x2=randint(sig_num,1,4);
x3=randint(sig_num,1,4);
x222=rou*x1+sqrt(1-(abs(rou))^2)*x3;%產生相關的兩路信號
x2=round(x222);
for i=1:sig_num
    if x2(i)>=4 
        x2(i)=3;
    elseif x2(i)<=0
        x2(i)=0;
    end
end
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;
for i=1:q
    for m=1:p
        a(m,i)=exp(-j*2*pi*0.5*(m-1)*sin(alpha(i)));
    end
end
a1=a(1:7,:);
a2=a(2:8,:);
x1=a1*s;
x2=a2*s;
for i=1:p-1
   xn1(i,:)=awgn(x1(i,:),snr,'measured');
   xn2(i,:)=awgn(x2(i,:),snr,'measured');
end
r1=xn1*xn1'./sig_num;
r2=xn1*xn2'./sig_num;
 [u1,s1,v1]=svd(r1);
cigma=min(diag(s1));
cxx=r1-cigma*eye(7,7);
z=[zeros(1,7);eye(6,6) zeros(6,1)];
cxy=r2-cigma*z;
% E1=eig(cxx,cxy);%基本的ESPRIT方法
% [ei1,index1]=sort(abs(abs(E1)-1));
% doa1=-asin(angle(E1(index1(1:2)))./pi)./pi*180
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TLS-ESPRIT方法
[U,S,V]=svd(cxx);
U1=U(:,1:2);%2為接受的信號數
V1=V(:,1:2);
S1=S(1:2,1:2);
cxy1=U1'*cxy*V1;
E=eig(S1,cxy1);
[ei,index]=sort(abs(abs(E)-1));
doa=-asin(angle(E(index(1:2)))./pi)./pi*180
% doaa=sort(doa);
% doa1=doaa(1);
% doa2=doaa(2);