clc;
%信源
%--------------------------------------------
Pd=2000;
Fd=1;
Fs=4*Fd;
R=0.5;
Delay=5;
M=4;
x1=randint(Pd,1,M); 
x2=randint(Pd,1,M);
y1=modmap(x1,Fd,Fs,'qask',M);
y2=modmap(x2,Fd,Fs,'qask',M);
[rcv_a1,ti]=rcosflt(y1,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a2,ti]=rcosflt(y2,Fd,Fs,'fir/sqrt/Fs',R,Delay);

s1=amodce(rcv_a1,10,'qam');
% s2=amodce(rcv_a2,10,'qam');
s2=s1; %假設信號s1和s2相干
save sig3 s1 s2
%--------------------------------------------

%全局變量定義
%--------------------------------------------
clear
m=8; %陣元數
angle1=-10;
angle2=10;
th=[angle1;angle2];
nn=1024;
SN1=10;
SN2=10;
SN3=12;
sn=[SN1;SN2];
degrad=pi/180;
%--------------------------------------------

%生成離散信號
%--------------------------------------------
load sig3 
tt=1:nn;
S=[s1(tt).';s2(tt).'];
nr=randn(m,nn);
ni=randn(m,nn);
U=nr+j*ni; %每路噪聲均值為2.
Ps=S*S'/nn;
ps=diag(Ps);
refp=2*10.^(sn/10);
tmp=sqrt(refp./ps);
S2=diag(tmp)*S;
%--------------------------------------------

%協方差矩陣及前后向平滑
%--------------------------------------------
tmp2=[0:m-1]'; %陣元位置
A=exp(-i*pi*tmp2*sin(th'*degrad)); %方向矩陣
X=A*S2+U;
Rxx=X*X'/nn;
a=eye(8);
b=a(:,8:-1:1);%構造交換矩陣
RXX=(Rxx+b*Rxx.'*b)/2;
%--------------------------------------------

%奇異值分解
%--------------------------------------------
[Q,SS,W]=svd(RXX); 
Vu=Q(:,3:8 ); %小特征值對應的特征向量，即噪聲子空間
%--------------------------------------------

%做圖
%--------------------------------------------
th2=[-90:90]';
tmp2=[0:m-1]';
a2=-i*tmp2*pi*sin(th2'*degrad);
A2=exp(a2);
num=diag(A2'*A2);
Ena=Vu'*A2;
den=diag(Ena'*Ena);
doa=num./den;
semilogy(th2,abs(doa));
axis([-90 90 0.1 1e5]);
grid on;
hold on;
%--------------------------------------------