论文：

贴代码

%2017.6.21
%L型阵列配对
%Pair-matching method for estimating 2d angle of arrival with a cross
%correlation matrix
clc
close all;
clear all

Mx=10;                                                                     %x轴阵元数
Mz=10;                                                                     %z轴阵元数
d_lamda=1/2;                                                               %阵元间距
K=4;                                                                       %信号数
phix=[45,52,120,150];                                                      %x轴方位角
thetaz=[15,25,65,88];                                                      %z轴俯仰角
fc=[200*10^6,300*10^6,500*10^6,400*10^6];                                  %信号频率
fs=1000*10^6;                                                              %采样频率
L=1000;                                                                    %快拍数
snr=20;                                                                    %信噪比

for k=1:L
    S(:,k)=sqrt(10.^(snr/10))*randn(1)*exp(1j*2*pi*fc'/fs*(k-1));
end
for kk=1:K
    A1x(:,kk)=exp(-1j*2*pi*[0:Mx-1]'*d_lamda*cos(phix(kk)/180*pi));
    A1z(:,kk)=exp(-1j*2*pi*[0:Mz-1]'*d_lamda*cos(thetaz(kk)/180*pi));      %1维导向矢量
end
%用ESPRIT算法求Z轴角度
ZZ=A1z*S+wgn(Mz,L,0,'dBm','complex');                                      %z轴接收到的信号
XX=A1x*S+wgn(Mz,L,0,'dBm','complex');                                      %x轴接收到的信号
Z1=ZZ(1:Mz-1,:);                                                           %Z子阵1接受的数据矢量
Z2=ZZ(2:Mz,:);                                                             %Z子阵2接受的数据矢量
X1=XX(1:Mx-1,:);                                                           %x子阵1接受的数据矢量
X2=XX(2:Mx,:);                                                             %x子阵2接受的数据矢量
Z=[Z1;Z2];                                                                 %对两个子阵的模型进行合并
X=[X1;X2];
Rz=Z*Z'/L;
Rx=X*X'/L;
[Uz,Sz,Vz]=svd(Rz);                                                        %对Rz,Rx进行奇异值分解
[Ux,Sx,Vx]=svd(Rx);
Rz=Rz-Sz(2*(Mz-1),2*(Mz-1))*eye(2*(Mz-1));                                 %去噪
Rx=Rx-Sx(2*(Mx-1),2*(Mx-1))*eye(2*(Mx-1));
[Uz1,Sz1,Vz1]=svd(Rz);
Uzs=Uz1(:,1:K);
Uzs1=Uzs(1:Mz-1,:);
Uzs2=Uzs(Mz:2*(Mz-1),:);
Mzz=pinv(Uzs1)*Uzs2;                                                       %按照公式得到旋转不变矩阵M
[Vzm,Dzm]=eig(Mzz);                                                        %对得到的旋转不变矩阵进行特征分解
Dzm=(diag(Dzm)).';
doaz=sort(180-acos(angle(Dzm)/pi)*180/pi,'ascend');
disp('doaz');
disp(doaz);                                                                %得到z轴估计的俯仰角
[Ux1,Sx1,Vx1]=svd(Rx);
Uxs=Ux1(:,1:K);
Uxs1=Uxs(1:Mx-1,:);
Uxs2=Uxs(Mx:2*(Mx-1),:);
Mxx=pinv(Uxs1)*Uxs2;
[Vxm,Dxm]=eig(Mxx);
Dxm=(diag(Dxm)).';
doax=sort(180-acos(angle(Dxm)/pi)*180/pi,'ascend');
disp('doax');
disp(doax);                                                                %得到x轴估计的方位角
rzx=diag((ZZ*XX')/L).';
rzx1=conj(rzx);
Rcc=toeplitz(rzx,rzx1);                                                    %构造Rcc
G=Rcc(:,1:K);
H=Rcc(:,K+1:Mx);
P=pinv(G)*H;                                                               %PM算法构造传播算子
Un=[P;-eye(Mx-K)];                                                         %噪声子空间
Gn=Un*Un';
searching_doa=0:1:180;                                                     %线阵的搜索范围为-90~90度
for i=1:length(searching_doa)
    a_theta=exp(-1j*2*pi*d_lamda*[0:Mx-1]'*sin(pi*searching_doa(i)/180));
    Pmusic(i)=1./abs((a_theta)'*Gn*a_theta);
end
%  plot(searching_doa,20*log10(Pmusic));
[Vw,Iw]=sort(Pmusic,'descend');
Iw=sort(Iw(:,1:K));
for k=1:K
    for kk=1:K
        w(k,kk)=cos(doaz(k)/180*pi)-cos(doax(kk)/180*pi);                  %构造K^2种组合（sin(θ）-sin(φ））
    end
end
w=w.';
n=numel(w);
wnew=reshape(w,1,n);                                                       %把K^2种组合变成一行
for k=1:K
    for kk=1:K^2
        f(k,kk)=abs((exp(1j*wnew(:,kk))-exp(1j*sin(Iw(k)/180*pi))));
    end
end
for k=1:K
    [Vest(k,:),Iest(k,:)]=sort(f(k,:));                                     %配对过程
end
Iest=Iest(1:K,1);
for k=1:K
    a=rem(Iest(k),K);
    if a~=0
        thetaest(k)=doax(a);
        b=(Iest(k)-a)/K;
        phiest(k)=doaz(b+1);
    else  thetaest(k)=doax(K);
        phiest(k)=doaz(K);
    end
end
disp('（俯仰角，方位角）');
for k=1:K
    disp([num2str(phiest(k)) ',' num2str(thetaest(k))]);
end