# 空间谱专题13：联合解算DOA（ML/AP）

该算法基于统计参数估计的思路，不涉及SVD分解或者相关矩阵求逆，因此对于相干信号理论上仍然适用，从理论的结构来看，由于投影矩阵涉及求逆，且有迭代过程，因此耗费资源过大。

function [phi_last,theta_last] = MuCalL_2D(x,srcNum,Array,resolution,lambda_c)
%L阵
sub1 = [1:6];
sub2 = [1,7:11];
J = fliplr(eye(length(sub1)));
x1 = x(sub1,:);
x2 = conj(J*x(sub2,:));
[phi,theta,spec1] = MuCalL_647_1D(x1,srcNum,Array(sub1,:),resolution,lambda_c);
[phi2,alpha,spec2] = MuCalL_647_1D(x2,srcNum,Array(sub1,:),resolution,lambda_c);

%
[val1,phi_pos] = findpeaks(spec1,'minpeakdistance',3);
[val,num_phi] = sort(val1,'descend');
[val2,theta_pos] = findpeaks(spec2,'minpeakdistance',3);
[val,num_theta] = sort(val2,'descend');
phi_est =  phi(phi_pos(num_phi(1:srcNum)));
theta_est =  alpha(theta_pos(num_theta(1:srcNum)));
phi_est = phi_est;
%筛选
snap = size(x,2);
R_all =  x*x'/snap;
para_all = perms([1:srcNum]);
theta_all = kron(theta_est,ones(1,size(para_all,1)));
theta_all([2,4]) = theta_all([4,2]);
phi_all = repmat(phi_est,1,size(para_all,1));
theta_all = asin(sin(theta_all/180*pi)./cos(phi_all/180*pi))/pi*180;

im = sqrt(-1);
Dd = [];
for kkk = 1:size(theta_all,2)/srcNum
nshift = ((kkk-1)*srcNum+1):((kkk)*srcNum);
theta_cache = theta_all(nshift)/180*pi;
phi_cache = phi_all(nshift)/180*pi;

Az = [];
for j = 1:srcNum
r = [sin(phi_cache(j)) cos(phi_cache(j))*sin(theta_cache(j)) cos(phi_cache(j))*cos(theta_cache(j))];
r_rep = repmat(r,size(x,1),1);
dis = sum(r_rep.*Array,2);
am = exp(-im*2*pi*dis/lambda_c);
Az = [Az,am];
end

Pb3 = Az*pinv(Az'*Az)*Az';
Dd(kkk) = abs(trace(Pb3*R_all));
end

[valDd,indexDd]=max(Dd);
n_pos = ((indexDd-1)*srcNum+1):((indexDd)*srcNum);
[phi_last,sort_pos] = sort(phi_all(n_pos),'ascend');
theta_last = theta_all(n_pos);
theta_last = theta_last(sort_pos);


联立解算的思路:

主要代码实现：

Ax = A(sub1,:);
Ay = A(sub2,:);
%利用T矩阵解算
y_sig = x2;
x_sig = x1;
Ryy = y_sig*y_sig'/snapshot;
Rs_hat = pinv(Ay'*Ay)*Ay'*Ryy*pinv(Ay*Ay')*Ay;%eq.5
Rxy = x_sig*y_sig'/snapshot;
Ay_pieH = pinv(Ay*Ay')*Ay;
Rs_pie = pinv(Ax'*Ax)*Ax'*Rxy*Ay_pieH;%eq.9
%构造T矩阵解算
perm = perms([1:srcNum]);
J = zeros(1,size(perm,1));
for i = 1:size(perm,1)
T = zeros(srcNum);
T(perm(i,:)+[0:srcNum-1]*srcNum) = 1;
J(i) = sum(sum(abs(Rs_pie-T*Rs_hat).^2));
end
[minVal,minPos] = min(J);
phi_est = phi_est(perm(minPos,:));
theta_all = theta_est;
%求解
phi_last = phi_est;
theta_last = asin(sin(theta_all/180*pi)./cos(phi_last/180*pi))/pi*180;

当个数不匹配的时候可参考（个人觉得直接拓展效果也可以，就是配对之前添加一个预处理）：

posted @ 2017-10-16 07:53 桂。 阅读(...) 评论(...) 编辑 收藏