基于均匀圆阵的MUSIC算法

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