clear?all;
close?all;
tic;
f=0.15*10^8;
c=3.0*10^8;
lmda=c/f;
r=0.6533*lmda;???%圓陣半徑
snr=10;
D=2;???%信號源數
kp=1024;
N=16;???%陣元數
%%%%%%%%%%陣元位置%%%%%%%%%%%
for?n=1:N
????px（n）=r*cos（2*pi*（n-1）/N）;
????py（n）=r*sin（2*pi*（n-1）/N）;
????pz（n）=0;
end
rad=pi/180;
sita=[50?75];??%輸入信號俯仰角
fei=[85?230];???%輸入信號方位角
%%%%%%輸入信號波數%%%%%%%%%%%%
for?n=1:D
????kx（n）=-2*pi/lmda*sin（sita（n）*rad）*cos（fei（n）*rad）;
????ky（n）=-2*pi/lmda*sin（sita（n）*rad）*sin（fei（n）*rad）;
????kz（n）=-2*pi/lmda*cos（sita（n）*rad）;
end
%%%%%%%%陣列流型矩陣%%%%%%%%%%%%
for?m=1:N
????for?n=1:D
????a（mn）=exp（1j*（px（m）*kx（n）+py（m）*ky（n）+pz（m）*kz（n）））;
????end
end
S=randn（Dkp）;
X=a*S;???
X1=awgn（Xsnr‘measured‘）;
R=X1*X1‘/kp;%協方差矩陣
J=zeros（NN）;
for?m=1:N
????for?n=1:N
????????if?m+n==N+1
????????????J（mn）=1;
????????else?J（mn）=0;
????????end
????end
end
for?m=1:N
????for?n=1:N
????????if?m<=N/2
????????????if?m==n
????????????????U（mn）=1;
????????????else?U（mn）=0;
????????????end
????????elseif?m>N/2
????????????if?m+n==3*N/2+1
????????????????U（mn）=1;
????????????else?U（mn）=0;
????????????end
????????end
????end
end
R1=U*R*U‘;
R2=J*conj（R1）*J;
R3=（R1+R2）/2;
RF=J*R3‘*J-R3;
Q=zeros（N:N）;
%%%%描述Q矩陣%%????圓陣與線陣的Q是不相同的
for?m=1:N
????for?n=1:N
????????if?m<=N/2
?