clear all;

Pd=2000;
Fd=1;
Fs=4*Fd;
R=0.5;
Delay=5;
M=4;
x1=randint(Pd,1,M); 
x2=randint(Pd,1,M);
x3=randint(Pd,1,M);
y1=modmap(x1,Fd,Fs,'qask',M);
y2=modmap(x2,Fd,Fs,'qask',M);
y3=modmap(x3,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);
[rcv_a3,ti]=rcosflt(y3,Fd,Fs,'fir/sqrt/Fs',R,Delay);
s1=amodce(rcv_a1,10,'qam'); 
s2=amodce(rcv_a2,10,'qam');
s3=amodce(rcv_a3,10,'qam');
save sig3 s1 s2 s3;

m=8; % 陣元數
p=3; % 信號數
%俯仰角為入射方向與圓陣軸線之間的夾角,方位角為x軸到入射線在圓陣上投影的夾角
st1=10; % 俯仰角
dt1=80; % 方位角
st2=20;
dt2=58;
st3=35;
dt3=69;

st=[st1;st2;st3];
dt=[dt1;dt2;dt3];

nn=4096; 
SN1=20; 
SN2=20;
SN3=20;
sn=[SN1;SN2;SN3];
degrad=pi/180;

load sig3
tt=1:nn;
S=[s1(tt),s2(tt),s3(tt)].';
noise=randn(m,nn)+j*randn(m,nn); 
Ps=S*S'/nn;
ps=diag(Ps);
refp=2*10.^(sn/10);
tmp=sqrt(refp./ps);
S2=diag(tmp)*S;

% 圓半徑為半波長
A=zeros(m,p);
k=[0:m-1]';
for t=1:p
A(:,t)=exp(j*pi*(sin(st(t)*degrad)*cos(dt(t)*degrad)*cos(2*pi*k/m)+sin(st(t)*degrad)*sin(dt(t)*degrad)*sin(2*pi*k/m)));
end
X=A*S2+noise;

Rxx=X*X'/nn;
[U,s,v]=svd(Rxx);
Vs=U(:,1:p);
Vu=U(:,p+1:m);

for st1=0:90
for dt1=0:90
k=[0:m-1]';
AA=exp(j*pi*(sin(st1*degrad)*cos(dt1*degrad)*cos(2*pi*k/8)+sin(st1*degrad)*sin(dt1*degrad)*sin(2*pi*k/8)));
WW=AA'*Vu*Vu'*AA;
Pmusic(st1+1,dt1+1)=abs(8./WW);%角譜
end 
end 
sst=0:90;
ddt=0:90;
mesh(sst,ddt,Pmusic);
title('MUSIC空間譜');
xlabel('方位角'); 
ylabel('俯仰角');
zlabel('角譜');
grid on;