DML（matlab代码）

【核心代码】

clear all;
close all;
clc;

source_number=2;%信元数
% sensor_number=13;%阵元数
 sensor_A=[0 4 10 13 17 23]; %52A天线布局
%  sensor_A=[0 1 6 13 14 19]; %52A天线布局
%    sensor_A=[0 2 6 9 11 15]; %博世天线布局
%   sensor_A=[0 1 6 19 20 25]; %均匀天线布局
 % sensor_A=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]; %互质天线布局
 sensor_number=length(sensor_A);
N_x=1; %信号长度
snapshot_number=N_x;%快拍数
w=[pi/4 pi/4].';%信号频率
l=((2*pi*3e8)/w(1) (2*pi*3e8)/w(2))/2;%信号波长  
d=0.5*l;%阵元间距
snr=26;%信噪比
j=sqrt(-1);
for n=1:1:100
source_doa=[0 15];%两个信号的入射角度
% A=[1*exp(-j*(0:sensor_number-1)*d*2*pi*sin(source_doa(1)*pi/180)/l);1*exp(-j*(0:sensor_number-1)*d*2*pi*sin(source_doa(2)*pi/180)/l)].';%阵列流型;以来波方向为source_doa（k）入射
% %M根天线，阵元间距为d的均匀线阵的阵列响应矢量
A=[1*exp(-j*sensor_A*d*2*pi*sin(source_doa(1)*pi/180)/l);0.1*exp(-j*sensor_A*d*2*pi*sin(source_doa(2)*pi/180)/l)].';
s=sqrt(10.^(snr/10))*exp(j*w*(0:N_x-1));%仿真信号
x=awgn(s,snr);
xx=A*s;
xxx=abs(xx);
% plot(xxx);
x=A*s (1/sqrt(2))*(randn(sensor_number,N_x) j*randn(sensor_number,N_x));%加了高斯白噪声后的阵列接收信号
R=x*x'/snapshot_number;
searching_doa=[-20:0.5:20]';
for m=1:1:length(searching_doa)
  for i=1:1:length(searching_doa)
     th=[searching_doa(m) searching_doa(i)];
%     th=[35,40];
     A1=exp(-j*sensor_A'*d*2*pi*sin(th.*pi/180)/l);
     p=A1*pinv(A1'*A1)*A1';
     xx=real(diag(p*R));
     ang(m,i)=sum(xx);
%     ang=sum(xx);
  end
end
[kk1,kk2]=find(ang==max(max(ang)));%寻找最大值所在位置
kk1=searching_doa(kk1);%寻找最大值对应的角度坐标
kk2=searching_doa(kk2);%寻找最大值对应的角度坐标
k1(n)=max(kk1);
k2(n)=min(kk2);
end
figure(2);
mesh(searching_doa,searching_doa,ang);
figure(1);
plot(k1,'y');
hold on;
plot(k2,'r');