clc;
clear all;

Pd=8000;
s1=1/1000*exp(j*pi/3*[0:Pd])';
s2=1/1000*exp(j*pi/4*[0:Pd])';
s3=1/1000*exp(j*pi/5*[0:Pd])';
save sig3 s1 s2 s3;

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

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)));
A(1,t)=1;
A(2,t)=exp(j*pi*cos(st(t)*degrad)*cos(dt(t)*degrad));
A(3,t)=exp(j*2*pi*cos(st(t)*degrad)*cos(dt(t)*degrad));
A(4,t)=exp(j*pi*sin(st(t)*degrad)*cos(dt(t)*degrad));
A(5,t)=exp(j*pi*cos(st(t)*degrad)*cos(dt(t)*degrad)+j*pi*sin(st(t)*degrad)*cos(dt(t)*degrad));
A(6,t)=exp(j*2*pi*cos(st(t)*degrad)*cos(dt(t)*degrad)+j*pi*sin(st(t)*degrad)*cos(dt(t)*degrad));
%構造A陣
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)));
AA=[1;exp(j*pi*cos(st1*degrad)*cos(dt1*degrad));exp(j*2*pi*cos(st1*degrad)*cos(dt1*degrad));
    exp(j*pi*sin(st1*degrad)*cos(dt1*degrad));exp(j*pi*cos(st1*degrad)*cos(dt1*degrad)+j*pi*sin(st1*degrad)*cos(dt1*degrad));
    exp(j*2*pi*cos(st1*degrad)*cos(dt1*degrad)+j*pi*sin(st1*degrad)*cos(dt1*degrad))];
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;