DOA代码框架

Posted 2020-11-28 czl-01

-------场景的搭建

M = 15；% the num of antenna
L = 100; % the num of snapshot
theta = [10 30]; %the DOA
K = length(theta); %the num of the signal
SNR = 10;
searchgrid = 0:0.01:90;

------信号模型的搭建

S = (randn(K,L)+j*randn(K,L))/sqrt(2);
loc = zeros(M,3);
loc(:,1) = 0:M-1;
A = zero(M,K)；
for k = 1:K
      rk = [cosd(phi(k))*sind(  theta(k)),sind(phi(k))*sind(theta(k)),cosd(theta(k))].‘;
      % rk = [cos(φk) sin(θk); sin(φk) sin(θk); cos(θk)]T 
       A(:,k) = exp(-j*pi*loc*rk);      
end
X0 = A*S;
Y = awgn(X0,SNR,‘measured‘,‘dB‘);

这里采取二维空间的信号模型，对于A中的元素赋值参考上诉文献。

------信号二阶统计特征的得到

R = Y*Y‘/L;
[V,D] = eig(R); %V：the characterist vector;D: the characterist value :all is matrix
DD = diag(D);
En = V(:,1:M-K);% the noise subspace

------DOA估计（以MUSIC为例）

 Un = En*En‘;%En‘转置操作
 temp = zeros(3,length(searchgrid1));%源代码中运用了三种估计方法
 phi1 = 0;%一维空间
 tic,
 for num1 = 1:length(searchgrid1)
      theta1 = searchgrid1(num1);
      r1 = [cosd(phi1)*sind(theta1),sind(phi1)*sind(theta1),cosd(theta1)].‘;
      aa = exp(-j*pi*loc*r1);% direction matrix
      temp(1,num1)=1/abs(aa‘*Un*aa);%the objective function
end
tem_2(1) = tem_2(1)+toc;% the time spend

------谱峰搜索

maxt1 = max(temp(1,:));
temp2 = 10*log(temp(1,:)/maxt1);%归一化并去对数使谱峰明显
[val2,indx2] =findpeaks(temp2);%val2为峰值大小，indx2为位置
[sval2,pos2] = sort(val2,‘descend‘);%sval2为降序排号的序列，pos2为索引地址
E_theta1 = sort(searchgrid1(indx2(pos2(1:K))));