clc;
clear all;

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

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

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;

% 圓半徑為半波長(zhǎng)
l=m/2;
%F=pi*sin
%jn(F)=
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;