DOA——MUSIC算法

一、均匀圆阵(UCA, Uniform Circular Array)的MUSIC算法

假设一个半径为R的M元均匀圆阵的所有阵元均位于坐标系X-Y平面内，第k-1个阵元坐标为，第i个窄带信号波长为，来波方向为，如图1，则第k-1个阵元到圆心(即原点)的波程差为：

均匀圆阵

存在P个入射信号均匀圆阵的接收模型可表示为：

其他步骤与基于ULA的MUSIC算法一致。

令任意两阵元间的波程差为：

当时，即产生相位模糊。将均匀圆阵各阵元投影到入射方向，得到一个随入射方向变动的非均匀线阵。需要保证在任意入射方向上投影出的非均匀线阵，其最小间隔总是小于信号波长，模糊谱峰对测向结果影响较小，即：

在方向角相同时，水平入射()信号的波程差最长，且投影出的非均匀线阵随方向角不同周期变化，因此只需要讨论水平入射信号对应投影线阵的不同情况。

在MUSIC算法中，阵元的最小间隔越大模糊谱峰峰值就越大。但在均匀圆阵中，阵元间隔随着入射波方向变化，因此算法性能受到最小间隔最大值的影响。根据来波方向不同，入射方向上的第k个阵元投影间隔分别为：

当M为奇数时，对着阵元入射，投影点重合为${ m{(M + 1)/2}}$个；当M为偶数时，对着相邻阵元连线中点入射，投影点重合为M/2个，此时投影线阵的非零最小间隔的值最大，且取得该最大值k=1时。进一步可求得半径的选取关系：

选取半径时，按上式等比例缩放，即能使对应的奇数阵与偶数阵有近似的抗相位模糊的性能。

每一路接收的结构图：

二、MUSIC

1、

clear all;
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsym=500;%符号个数
Fsym=1;%符号速率
M=3;%一个符号对应的比特数
Fbit=M*Fsym;%比特速率
Nsour=3;%信源数
Angle=[5,15,35];%信源的来波方向
Fc=10;%载波频率
Fs=100;%抽样频率
R=0.5;%滚降因子
Del=5;%群延迟因子
% Nsamp=50;%采样点数或者快拍数
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%构造相干信源--信源1、信源2与信源3
Rcos2=Rcos11+Rcos21+Rcos31;%构造相干信源--信源1、信源2与信源3
Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相干信源--信源1、信源2与信源3
save xyc3 Rcos1 Rcos2 Rcos3
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsamp=1024;%采样点数或者快拍数 
i=sqrt(-1);
j=i;
Ntx=8;%阵列数
SNR=[2,2,2];%三信源的信噪比
% sn=10;                        %----单信号源
Lamda=2;%信号波长
D=Lamda/2;%线性阵列的距离
p=3;%子阵个数
L=Ntx-p+1;%子阵阵元数
nr=randn(Ntx,Nsamp);
ni=randn(Ntx,Nsamp);
n=nr+j*ni;%产生背景噪声
load xyc3;
t=1:Nsamp;
% s1=[Rcos1(t).'];%接收信号的采样点数%----单信号源
s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩阵维数=信源数*抽样点数
ps=diag((s1*s1')/Nsamp);%无噪声信号功率--%矩阵维数=信源数*1
delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1
% delta1=ps./(2*10.^(sn/10));        %----单信号源
delta2=diag(delta1);%矩阵维数=1*1
delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1
Rev_s1=(1./delta')*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数
%计算各信源SNR比条件下，阵列接收到的信号幅度%
Pn=zeros(Nsamp,1);
pn=zeros(Ntx,Nsamp);
Pn=diag(n'*n);
for h=1:Nsamp
    pn(:,h)=n(:,h)./sqrt(Pn(h,:));
end
Rev_n=pn;
%计算各阵列接收到的背景噪声下的信号幅度%
tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数
% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda;   %----单信号源
tmp1=[0:Ntx-1]';%矩阵维数=阵元数*1
tmp4=[0:L-1]';%子矩阵维数=子矩阵阵元数*1
a1=tmp1*tmp;%矩阵维数=阵元数*信源数
A=exp(a1);%方向矩阵--%矩阵维数=阵元数*信源数
X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数
Rxx=(X*X')/Nsamp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sub_FRxx=zeros(L,L);
sub_BRxx=zeros(L,L);
for i=1:p
    sub_FR=zeros(L,Nsamp);
    sub_BR=zeros(L,Nsamp);
    sub_FR=X(i:1:i+L-1,:);
    last=Ntx+1-i;
    first=last-L+1;
    sub_BR=conj(X(last:-1:first,:));
    sub_FRxx=sub_FRxx+((sub_FR*sub_FR')./Nsamp);
    sub_BRxx=sub_BRxx+((sub_BR*sub_BR')./Nsamp);
end
sub_FRxx=sub_FRxx./p;
sub_BRxx=sub_BRxx./p;
sub_Rxx=(sub_FRxx+sub_BRxx)./2;
[VFB,HFB]=eig(sub_Rxx);
[HFB,IFB]=sort(diag(HFB),1);
VFB=VFB(:,IFB);
VnFB=VFB(:,1:L-Nsour);
VsFB=VFB(:,L-Nsour+1:L);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ScanAng=[-90:1:90];
for i=1:length(ScanAng)
    tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;
    tmp3=tmp2*tmp1;
    tmp5=tmp2*tmp4;
    A_Sita=exp(tmp3);
    Sub_Sita=exp(tmp5);
    Sub_FBsita(i)=(Sub_Sita'*Sub_Sita)/(Sub_Sita'*VnFB*VnFB'*Sub_Sita);
end 
figure(1);
semilogy(ScanAng,real(Sub_FBsita),'bo-');
axis([-60 60 0.1 1e7]);
xlabel('M_Angle(deg)');
ylabel('M_Spectrum');
grid on

　　

2、

clear all;
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsym=500;%符号个数
Fsym=1;%符号速率
M=3;%一个符号对应的比特数
Fbit=M*Fsym;%比特速率
Nsour=3;%信源数
Angle=[10,40,80];%信源的来波方向
Fc=10;%载波频率
Fs=100;%抽样频率
R=0.5;%滚降因子
Del=5;%群延迟因子
% Nsamp=50;%采样点数或者快拍数
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=Rcos11;%构造相干信源--信源1、信源2与信源3
Rcos2=Rcos21;%构造相干信源--信源1、信源2与信源3
Rcos3=Rcos31;%构造相干信源--信源1、信源2与信源3
save xyc3 Rcos1 Rcos2 Rcos3
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsamp=1024;%采样点数或者快拍数 
i=sqrt(-1);
j=i;
Ntx=8;%阵列数
SNR=[10,10,10];%三信源的信噪比
% sn=10;                        %----单信号源
Lamda=2;%信号波长
D=Lamda/2;%线性阵列的距离
p=3;%子阵个数
L=Ntx-p+1;%子阵阵元数
nr=randn(Ntx,Nsamp);
ni=randn(Ntx,Nsamp);
n=nr+j*ni;%产生背景噪声
load xyc3;
t=1:Nsamp;
% s1=[Rcos1(t).'];%接收信号的采样点数%----单信号源
s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩阵维数=信源数*抽样点数
ps=diag((s1*s1')/Nsamp);%无噪声信号功率--%矩阵维数=信源数*1
delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1
% delta1=ps./(2*10.^(sn/10));        %----单信号源
delta2=diag(delta1);%矩阵维数=1*1
delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1
Rev_s1=(1./delta')*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数
%计算各信源SNR比条件下，阵列接收到的信号幅度%
Pn=zeros(Nsamp,1);
pn=zeros(Ntx,Nsamp);
Pn=diag(n'*n);
for h=1:Nsamp
    pn(:,h)=n(:,h)./sqrt(Pn(h,:));
end
Rev_n=pn;
%计算各阵列接收到的背景噪声下的信号幅度%
tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数
% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda;   %----单信号源
tmp1=[0:Ntx-1]';%矩阵维数=阵元数*1
tmp4=[0:L-1]';%子矩阵维数=子矩阵阵元数*1
a1=tmp1*tmp;%矩阵维数=阵元数*信源数
A=exp(a1);%方向矩阵--%矩阵维数=阵元数*信源数
X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数
Rxx=(X*X')/Nsamp;
[V,H]=eig(Rxx);%MUSIC算法---MUltiSIgnal Classification
[H,I]=sort(diag(H),1);%特征值按照升序排列
V=V(:,I);%特征值对应的特征向量也按照相应特征值的升序排列
Vn=V(:,1:Ntx-Nsour);%噪声子空间---协方差的特征向量--最小特征值对应的特征向量
Vs=V(:,Ntx-Nsour+1:Ntx);%信号子空间---协方差的特征向量--最大特征值对应的特征向量
ScanAng=[-90:1:90];
for i=1:length(ScanAng)
    tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;
    tmp3=tmp2*tmp1;
    tmp5=tmp2*tmp4;
    A_Sita=exp(tmp3);
MUSIC_Sita(i)=(A_Sita'*A_Sita)/(A_Sita'*Vn*Vn'*A_Sita);
end 
figure(1);
semilogy(ScanAng,real(MUSIC_Sita),'g*-');
axis([-90 90 0.1 1e7]);
xlabel('M_Angle(deg)');
ylabel('M_Spectrum');
grid on

　　

3、

clear all;
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsym=500;%符号个数
Fsym=1;%符号速率
M=3;%一个符号对应的比特数
Fbit=M*Fsym;%比特速率
Nsour=3;%信源数
Angle=[5,8,35];%信源的来波方向
Fc=10;%载波频率
Fs=100;%抽样频率
R=0.5;%滚降因子
Del=5;%群延迟因子
% Nsamp=50;%采样点数或者快拍数
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%构造相干信源--信源1、信源2与信源3
Rcos2=Rcos11+Rcos21+Rcos31;%构造相干信源--信源1、信源2与信源3
Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%构造相干信源--信源1、信源2与信源3
save xyc3 Rcos1 Rcos2 Rcos3
%产生三信源，角度分别为-40°、30°、45°，采用8PSK调制，滚降系数为0.5的平方根升余弦滤波
Nsamp=512;%采样点数或者快拍数 
i=sqrt(-1);
j=i;
Ntx=8;%阵列数
SNR=[2,2,2];%三信源的信噪比
% sn=10;                        %----单信号源
Lamda=2;%信号波长
D=Lamda/2;%线性阵列的距离
p=3;%子阵个数
L=Ntx-p+1;%子阵阵元数
nr=randn(Ntx,Nsamp);
ni=randn(Ntx,Nsamp);
n=nr+j*ni;%产生背景噪声
load xyc3;
t=1:Nsamp;
% s1=[Rcos1(t).'];%接收信号的采样点数%----单信号源
s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩阵维数=信源数*抽样点数
ps=diag((s1*s1')/Nsamp);%无噪声信号功率--%矩阵维数=信源数*1
delta1=(1./(2*10.^(SNR/10)))*ps;%矩阵维数=1*1
% delta1=ps./(2*10.^(sn/10));        %----单信号源
delta2=diag(delta1);%矩阵维数=1*1
delta=sqrt(delta2);%噪声幅度值--%矩阵维数=1*1
Rev_s1=(1./delta')*s1;%SNR条件下的信号幅度--%矩阵维数=信源数*抽样点数
%计算各信源SNR比条件下，阵列接收到的信号幅度%
Pn=zeros(Nsamp,1);
pn=zeros(Ntx,Nsamp);
Pn=diag(n'*n);
for h=1:Nsamp
    pn(:,h)=n(:,h)./sqrt(Pn(h,:));
end
Rev_n=pn;
%计算各阵列接收到的背景噪声下的信号幅度%
tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩阵维数=1*信源数
% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda;   %----单信号源
tmp1=[0:Ntx-1]';%矩阵维数=阵元数*1
tmp4=[0:L-1]';%子矩阵维数=子矩阵阵元数*1
a1=tmp1*tmp;%矩阵维数=阵元数*信源数
A=exp(a1);%方向矩阵--%矩阵维数=阵元数*信源数
X=A*Rev_s1+Rev_n;%阵列接收到的信号幅度--%矩阵维数=阵元数*抽样点数
Rxx=(X*X')/Nsamp;
Rxx_fb=zeros(L,L);
Rxx_f=zeros(L,L);
Rxx_b=zeros(L,L);
J=fliplr(eye(L));
for m=1:p
    for k=1:p
        Rxx_f=Rxx_f+Rxx(m:1:m+L-1,k:1:k+L-1)*Rxx(k:1:k+L-1,m:1:m+L-1);
        Rxx_b=Rxx_b+J*conj(Rxx(m:1:m+L-1,k:1:k+L-1))*conj(Rxx(k:1:k+L-1,m:1:m+L-1))*J;
    end
end
Rxx_f=Rxx_f./p;
Rxx_b=Rxx_b./p;
Rxx_fb=(Rxx_f+Rxx_b)./p;
[V_fb,H_fb]=eig(Rxx_fb);%特征分解---MUltiSIgnal Classification
[H_fb,I_fb]=sort(diag(H_fb),1);%特征值按照升序排列
V_fb=V_fb(:,I_fb);%特征值对应的特征向量也按照相应特征值的升序排列
Vn_fb=V_fb(:,1:L-Nsour);
Vs_fb=V_fb(:,L-Nsour+1:L);
ScanAng=[-90:1:90];
for i=1:length(ScanAng)
    tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;
    tmp3=tmp2*tmp1;
    tmp5=tmp2*tmp4;
    A_Sita=exp(tmp3);
    Sub_Sita=exp(tmp5);
    fb_sita(i)=(Sub_Sita'*Sub_Sita)/(Sub_Sita'*Vn_fb*Vn_fb'*Sub_Sita);
end 
figure(1);
semilogy(ScanAng,real(fb_sita),'r-');
axis([-60 60 0.1 1e7]);
xlabel('M_Angle(deg)');
ylabel('M_Spectrum');
grid on

　　

4、

%=========================================================================
%      UCA_multi_in_2D
%       
%=========================================================================

clc;
clear all;
close all;

%------------------------常数表-------------------------------
c = 3e8;
namda = c/18e9;
est_num = 1;
iteration = 100;
sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5]; 
phi = 60;

%% ----------------入射信号模型-------------------------------
N_x = 2^5;                 %快拍点数
F0 = 18e9;                   %中心频率
B = 20e6;                   %带宽
Fs = 2*B;                   %采样频率
Ts = 1/Fs;                  %采样时间
T = (N_x-1)*Ts;             %快拍持续时间
u = B/T;                    %频率变化率
t = -T/2:Ts:T/2;            %时间轴点
l = c/18e9;
st = exp(1j*2*pi*(F0*t+.5*u*t.^2));

dir7 =(46:.25:56)*pi/180;  %(-50:.25:-40)*pi/180;
dir8 =(62.5:.25:72.5)*pi/180;  
dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;
dir10=(49:.25:59)*pi/180;  
ang=(50:.25:70)*pi/180;

e_dir7 = zeros(1,length(sr_array));
e_dir8 = zeros(1,length(sr_array));
e_dir9 = zeros(1,length(sr_array));
e_dir10= zeros(1,length(sr_array));


for ss = 1:length(sr_array)
    snr  = sr_array(ss);
    %--------------------7阵元---------------------------------------
    sensor_num = 7;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    R =  namda/(1-cos(d_angle(2)));
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = d_angle(2)*180/pi;
    for it = 1:iteration
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));            
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        xt = A*(sqrt(10^(snr/10))*st)+n;                
    % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir7),length(ang));
        for i=1:length(dir7)
            for k=1:length(ang)
               a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i)));
               a_theta = exp(-1j*2*pi*a_tao/l);
               Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end
        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir7(a)*180/pi;
        e_dir7(ss) = e_dir7(ss)+(aa-theta)^2;
    end
    
    %--------------------8阵元-------------------------------------
    sensor_num = 8;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = 1.5*d_angle(2)*180/pi;
    for it = 1:iteration         
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
        xt = A*(sqrt(10^(snr/10))*st)+n;
        % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir8),length(ang));
        for i=1:length(dir8)
            for k=1:length(ang)
                a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i)));
                a_theta = exp(-1j*2*pi*a_tao/l);
                Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end

        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir8(a)*180/pi;
        e_dir8(ss) = e_dir8(ss)+(aa-theta)^2;
    end
    
    
    
    %---------------9阵元---------------------------------------------------
    sensor_num = 9;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    R =  namda/(1-cos(d_angle(2)));
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = d_angle(2)*180/pi;
    for it = 1:iteration
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));            
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        xt = A*(sqrt(10^(snr/10))*st)+n;                
    % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir9),length(ang));
        for i=1:length(dir9)
            for k=1:length(ang)
               a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i)));
               a_theta = exp(-1j*2*pi*a_tao/l);
               Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end
        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir9(a)*180/pi;
        e_dir9(ss) = e_dir9(ss)+(aa-theta)^2;
    end
    
    %---------------10阵元---------------------------------------------------
    sensor_num = 10;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = 1.5*d_angle(2)*180/pi;
    for it = 1:iteration         
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
        xt = A*(sqrt(10^(snr/10))*st)+n;
        % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir10),length(ang));
        for i=1:length(dir10)
            for k=1:length(ang)
                a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i)));
                a_theta = exp(-1j*2*pi*a_tao/l);
                Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end

        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir10(a)*180/pi;
        e_dir10(ss) = e_dir10(ss)+(aa-theta)^2;
    end
    
end
figure;
% subplot(121);
% plot(sr_array+45,e_dir10(1,:)/iteration,'-^k',sr_array+45,e_dir9(1,:)/iteration,'-*k');
% legend('8元UCA','9元UCA');
% grid on;%axis([-10,20,-.01,.45]);
% colormap gray;
% xlabel('信噪比/dB');ylabel('均方误差/°');
% title('半径相同精度试验');
% 
% 
% subplot(122);
plot(sr_array+45,e_dir7/iteration,'-^k',sr_array+45,e_dir8/iteration,'-*k');hold on;
plot(sr_array+45,e_dir9/iteration,'-sk',sr_array+45,e_dir10/iteration,'-dk');
legend('7阵元','8阵元','9阵元','10阵元');
grid on;axis([-5.5,18,-.02,.65]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方误差/°');
title('不同阵元最大半径测向试验');

　　

5、

%=========================================================================
%      UCA_multi_in_2D
%       
%=========================================================================

clc;
clear all;
close all;

%------------------------常数表-------------------------------
c = 3e8;
namda = c/18e9;
est_num = 1;
iteration = 100;
sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5]; 
phi = 60;

%% ----------------入射信号模型-------------------------------
N_x = 2^5;                 %快拍点数
F0 = 18e9;                   %中心频率
B = 20e6;                   %带宽
Fs = 2*B;                   %采样频率
Ts = 1/Fs;                  %采样时间
T = (N_x-1)*Ts;             %快拍持续时间
u = B/T;                    %频率变化率
t = -T/2:Ts:T/2;            %时间轴点
l = c/18e9;
st = exp(1j*2*pi*(F0*t+.5*u*t.^2));
R = 0.01;
dir9(1,:)=(-50:.25:-40)*pi/180;
dir9(2,:)=(75:.25:85)*pi/180;  
dir10(1,:)=(-51:.25:-41)*pi/180;
dir10(2,:)=(107:.25:117)*pi/180;  
ang=(50:.25:70)*pi/180;

e_dir9= zeros(2,length(sr_array));
e_ang9= zeros(2,length(sr_array));
e_dir10= zeros(2,length(sr_array));
e_ang10= zeros(2,length(sr_array));


for ss = 1:length(sr_array)
    snr  = sr_array(ss);
    %---------------7阵元---------------------------------------------------
    sensor_num = 9;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    x = R*cos(d_angle);y = R*sin(d_angle);
    theta_array = [-9*d_angle(2)/8,d_angle(3)]*180/pi;

    for it = 1:iteration
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
        for dd = 1:2               
            theta = theta_array(dd);
            tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
            A = exp(-1j*2*pi*tao./l);
            
            xt = A*(sqrt(10^(snr/10))*st)+n;                
        % -------------------2D-MUSIC算法-----------------------
            Rxx = xt*xt'/N_x;
            [U,S] = svd(Rxx);
            [~,index] = sort(diag(S));
            Un = U(:,index(1:sensor_num-est_num));
            Gn = Un*Un';           
            Pmusic = zeros(length(dir9(dd,:)),length(ang));
            for i=1:length(dir9(dd,:))
                for k=1:length(ang)
                   a_tao = sin(ang(k))*(x*cos(dir9(dd,i))+y*sin(dir9(dd,i)));
                   a_theta = exp(-1j*2*pi*a_tao/l);
                   Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
                end
            end

            [a,b]=find(Pmusic==max(max(Pmusic)));
            aa = dir9(dd,a)*180/pi;
            bb = ang(b)*180/pi;
            e_dir9(dd,ss) = e_dir9(dd,ss)+(aa-theta)^2;
            e_ang9(dd,ss) = e_ang9(dd,ss)+(bb-phi)^2;
        end
    end
       %---------------8阵元---------------------------------------------------
    sensor_num = 8;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    x = R*cos(d_angle);y = R*sin(d_angle);
    theta_array = [-9*d_angle(2)/8,2.5*d_angle(2)]*180/pi;
    for it = 1:iteration
        
        for dd = 1:2            
            theta = theta_array(dd);
            tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
            A = exp(-1j*2*pi*tao./l);
            n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
            xt = A*(sqrt(10^(snr/10))*st)+n;
            % -------------------2D-MUSIC算法-----------------------
            Rxx = xt*xt'/N_x;
            [U,S] = svd(Rxx);
            [~,index] = sort(diag(S));
            Un = U(:,index(1:sensor_num-est_num));
            Gn = Un*Un';           
            Pmusic = zeros(length(dir10(dd,:)),length(ang));
            for i=1:length(dir10(dd,:))
                for k=1:length(ang)
                    a_tao = sin(ang(k))*(x*cos(dir10(dd,i))+y*sin(dir10(dd,i)));
                    a_theta = exp(-1j*2*pi*a_tao/l);
                    Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
                end
            end

            [a,b]=find(Pmusic==max(max(Pmusic)));
            aa = dir10(dd,a)*180/pi;
            bb = ang(b)*180/pi;
            e_dir10(dd,ss) = e_dir10(dd,ss)+(aa-theta)^2;
            e_ang10(dd,ss) = e_ang10(dd,ss)+(bb-phi)^2;
        end
    end
end
figure;
subplot(121);
plot(sr_array+45,e_dir9(1,:)/iteration,'-^k',sr_array+45,e_dir9(2,:)/iteration,'-*k');
legend('-9/8 360/M','3 360/M');
grid on;%axis([-10,20,-.01,.45]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方误差/°');
title('9元阵方向角误差');
% 
% subplot(222);
% plot(sr_array+45,e_ang9(1,:)/iteration,'-^k',sr_array+45,e_ang9(2,:)/iteration,'-*k');
% legend('-9/8 360/M','3 360/M');
% grid on;axis([-10,20,-.1,1.45]);
% colormap gray;
% xlabel('信噪比/dB');ylabel('均方误差/°');
% title('7元阵俯仰角误差');


subplot(122);
plot(sr_array+45,e_dir10(1,:)/iteration,'-^k',sr_array+45,e_dir10(2,:)/iteration,'-*k');
legend('-5/4 360/M','5/2 360/M');
grid on;%axis([-10,20,-.01,.45]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方误差/°');
title('8元阵方向角误差');
% 
% 
% subplot(224);
% plot(sr_array+45,e_ang10(1,:)/iteration,'-^k',sr_array+45,e_ang10(2,:)/iteration,'-*k');
% legend('-5/4 360/M','3 360/M');
% grid on;axis([-10,20,-.1,1.45]);
% colormap gray;
% xlabel('信噪比/dB');ylabel('均方误差/°');
% title('8元阵俯仰角误差');

　　

6、

%=========================================================================
%      UCA_multi_in_2D
%       
%=========================================================================

clc;
clear all;
close all;

%------------------------常数表-------------------------------
c = 3e8;
namda = c/18e9;
est_num = 1;
iteration = 100;
sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5]; 
phi = 60;

%% ----------------入射信号模型-------------------------------
N_x = 2^5;                 %快拍点数
F0 = 18e9;                   %中心频率
B = 20e6;                   %带宽
Fs = 2*B;                   %采样频率
Ts = 1/Fs;                  %采样时间
T = (N_x-1)*Ts;             %快拍持续时间
u = B/T;                    %频率变化率
t = -T/2:Ts:T/2;            %时间轴点
l = c/18e9;
st = exp(1j*2*pi*(F0*t+.5*u*t.^2));

dir7 =(46:.25:56)*pi/180;  %(-50:.25:-40)*pi/180;
dir8 =(62.5:.25:72.5)*pi/180;  
dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;
dir10=(49:.25:59)*pi/180;  
ang=(50:.25:70)*pi/180;

e_dir7 = zeros(1,length(sr_array));
e_dir8 = zeros(1,length(sr_array));
e_dir9 = zeros(1,length(sr_array));
e_dir10= zeros(1,length(sr_array));
R = 0.1;

for ss = 1:length(sr_array)
    snr  = sr_array(ss);
    %--------------------7阵元---------------------------------------
    sensor_num = 7;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = d_angle(2)*180/pi;
    for it = 1:iteration
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));            
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        xt = A*(sqrt(10^(snr/10))*st)+n;                
    % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir7),length(ang));
        for i=1:length(dir7)
            for k=1:length(ang)
               a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i)));
               a_theta = exp(-1j*2*pi*a_tao/l);
               Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end
        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir7(a)*180/pi;
        e_dir7(ss) = e_dir7(ss)+(aa-theta)^2;
    end
    
    %--------------------8阵元-------------------------------------
    sensor_num = 8;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = 1.5*d_angle(2)*180/pi;
    for it = 1:iteration         
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
        xt = A*(sqrt(10^(snr/10))*st)+n;
        % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir8),length(ang));
        for i=1:length(dir8)
            for k=1:length(ang)
                a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i)));
                a_theta = exp(-1j*2*pi*a_tao/l);
                Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end

        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir8(a)*180/pi;
        e_dir8(ss) = e_dir8(ss)+(aa-theta)^2;
    end
    
    
    
    %---------------9阵元---------------------------------------------------
    sensor_num = 9;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = d_angle(2)*180/pi;
    for it = 1:iteration
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));            
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        xt = A*(sqrt(10^(snr/10))*st)+n;                
    % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir9),length(ang));
        for i=1:length(dir9)
            for k=1:length(ang)
               a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i)));
               a_theta = exp(-1j*2*pi*a_tao/l);
               Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end
        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir9(a)*180/pi;
        e_dir9(ss) = e_dir9(ss)+(aa-theta)^2;
    end
    
    %---------------10阵元---------------------------------------------------
    sensor_num = 10;
    d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
    x = R*cos(d_angle);y = R*sin(d_angle);
    
    theta = 1.5*d_angle(2)*180/pi;
    for it = 1:iteration         
        tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
        A = exp(-1j*2*pi*tao./l);
        n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
        xt = A*(sqrt(10^(snr/10))*st)+n;
        % -------------------2D-MUSIC算法-----------------------
        Rxx = xt*xt'/N_x;
        [U,S] = svd(Rxx);
        [~,index] = sort(diag(S));
        Un = U(:,index(1:sensor_num-est_num));
        Gn = Un*Un';           
        Pmusic = zeros(length(dir10),length(ang));
        for i=1:length(dir10)
            for k=1:length(ang)
                a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i)));
                a_theta = exp(-1j*2*pi*a_tao/l);
                Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
            end
        end

        [a,b]=find(Pmusic==max(max(Pmusic)));
        aa = dir10(a)*180/pi;
        e_dir10(ss) = e_dir10(ss)+(aa-theta)^2;
    end
    
end
figure;

plot(sr_array+45,e_dir7/iteration,'-^k',sr_array+45,e_dir8/iteration,'-*k');hold on;
plot(sr_array+45,e_dir9/iteration,'-sk',sr_array+45,e_dir10/iteration,'-dk');
legend('7阵元','8阵元','9阵元','10阵元');
grid on;axis([-5.5,18,-.02,.25]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方误差/°');
title('不同阵元相同半径测向试验');

　　

7、

%=========================================================================
%      Circular Array Classical-Music
%       
%=========================================================================

clc;
clear all;
close all;

c = 3e8;
phi = 60;
namda = c/18e9;
R = 7.5/100;
R = 10/100;

snr = -35;                   %信噪比
N_x = 2^5;                  %快拍点数
F0 = 18e9;                  %中心频率
B = 20e6;                   %带宽
Fs = 40e6;                  %采样频率
Ts = 1/Fs;                  %采样时间
T = (N_x-1)*Ts;             %快拍持续时间
u = B/T;                    %频率变化率
t = -T/2:Ts:T/2;            %时间轴点
l = c/F0;

st = sqrt(10^(snr/10))*exp(1j*2*pi*(F0*t+.5*u*t.^2));
%% -----------------9------------------------
sensor_num = 9;
R = 7.1239/100;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
theta =d_angle(2)*180/pi;
x = R*cos(d_angle);y = R*sin(d_angle);   

tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));
A = exp(-1j*2*pi*tao/l);
n = randn(sensor_num,N_x)*sqrt(10^(-45/10));
xt = A*st+n;
% -------------------2D-MUSIC算法-----------------------
Rx = xt*xt'/N_x;
[U,S] = eig(Rx);
est_sour = 1;
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]);
Gn = Un*Un';

dir=(-180:.25:179.8)*pi/180;  
ang=(20:.25:91)*pi/180;
Pmusic9 = zeros(length(dir),length(ang));

for i=1:length(dir)
    for k=1:length(ang)
        a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i)));
        a_theta = exp(-1j*2*pi*a_tao/l);
        Pmusic9(i,k)=1./abs((a_theta)'*Gn*a_theta);
    end
end
P_music9 = 10*log10(Pmusic9/min(min(Pmusic9)));

%% -----------------9------------------------
sensor_num = 10;
R = 4.5879/100;
%R = 10/100;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
theta =1.6*d_angle(2)*180/pi;
x = R*cos(d_angle);y = R*sin(d_angle);   

tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));
A = exp(-1j*2*pi*tao/l);
n = randn(sensor_num,N_x)*sqrt(10^(-45/10));
xt = A*st+n;
% -------------------2D-MUSIC算法-----------------------
Rx = xt*xt'/N_x;
[U,S] = eig(Rx);
disp(est_sour);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]);
Gn = Un*Un';

dir=(-180:.25:179.8)*pi/180;  
ang=(20:.25:91)*pi/180;
Pmusic10 = zeros(length(dir),length(ang));

for i=1:length(dir)
    for k=1:length(ang)
        a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i)));
        a_theta = exp(-1j*2*pi*a_tao/l);
        Pmusic10(i,k)=1./abs((a_theta)'*Gn*a_theta);
    end
end
P_music10 = 10*log10(Pmusic10/min(min(Pmusic10)));


figure;
% subplot(221);
% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
% mesh(xx,yy,P_music9);
% title('9元阵二维空间谱');
% xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空间谱/dB');
% axis([20,91,-180,180,0,24]);%colormap gray;

subplot(121);
[xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
mesh(xx,yy,P_music9);
title('9元阵方向角空间谱');
xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空间谱/dB');
axis([20,91,-180,180,0,24]);%colormap gray;

% subplot(222);
% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
% mesh(xx,yy,P_music10);
% title('10元阵二维空间谱');
% xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空间谱/dB');
% axis([20,91,-180,180,0,24]);%colormap gray;

subplot(122);
[xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
mesh(xx,yy,P_music10);
title('10元阵方向角空间谱');
xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空间谱/dB');
axis([20,91,-180,180,0,24]);%colormap gray;

　　

8、

% 波程差图
clear all;
close all;
clc;


sensor_num = 9;
c = 3e8;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
phi = 90;

namda = c/18e9;
    if mod(sensor_num,2);
        R =  namda/(1-cos(d_angle(2)));%相邻最小间距最大值
        %R = namda/(cos(floor(sensor_num/4)*d_angle(2))+cos(d_angle(2)*floor((sensor_num+1)/4)-.5));%相邻最大间距最大值
        %R = namda/(1+cos(d_angle(2)/2));%任意阵元间距最大
    else
        R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));%相邻最小间距最大值        
        %R = namda/(2*sin(d_angle(2)/2));%相邻最大间距最大值
        %R = namda/2;%任意阵元间距最大值
    end
    
R = 1;

dir = (0:.4:d_angle(2)*180/pi)*pi/180;
%idx=nchoosek(1:sensor_num,2);  % 16取2的组合
idx1 = [1,2;2,9;9,3;3,8;8,4;4,7;7,5;5,6];
idx2 = [2,1;1,3;3,9;9,4;4,8;8,5;5,7;7,6];
r = zeros(length(dir),size(idx1,1));
x = R*cos(d_angle);y = R*sin(d_angle);
for i = 1:50
    tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
    r(i,:) = -diff(tao(idx1),1,2);
end
for i = 51:101
    tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
    r(i,:) = -diff(tao(idx2),1,2);
end
figure;
subplot(211);
plot(dir(1:5:length(dir))*180/pi,0.234*ones(1,length(dir(1:5:length(dir)))),'.k',dir*180/pi,r(:,end),'-.k',dir*180/pi,r(:,1:end-1),'k');grid on;
legend('0.234','k=1','k=其他');
title('(a)9元阵投影间隔')
xlabel('入射方向°');ylabel('投影间隔/R');
axis([0,40,-0.01,0.7]);

% ---------------------------------------------------------------------
sensor_num = 10;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;

dir = linspace(0,d_angle(2),101);
%idx=nchoosek(1:sensor_num,2);  % 16取2的组合
idx1 = [1,2;2,10;10,3;3,9;9,4;4,8;8,5;5,7;7,6];
idx2 = [2,1;1,3;3,10;10,4;4,9;9,5;5,8;8,6;6,7];
r = zeros(length(dir),size(idx1,1));
x = R*cos(d_angle);y = R*sin(d_angle);
for i = 1:50
    tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
    r(i,:) = -diff(tao(idx1),1,2);
end
for i = 51:101
    tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
    r(i,:) = -diff(tao(idx2),1,2);
end
subplot(212);
plot(dir(1:5:length(dir))*180/pi,0.3633*ones(1,length(dir(1:5:length(dir)))),'.k',dir*180/pi,r(:,2),'-.k',dir*180/pi,r(:,[1,3:6,7]),'k');grid on;
legend('0.3633','k=1','k=其他');
title('(b)10元阵投影间隔');
xlabel('入射方向/°');ylabel('投影间隔/R');
axis([0,36,-0.01,0.7]);