基于原子范数最小化的互质阵列相干源DOA估计

From : DOA Estimation of Coherent Sources Using Coprime Array via Atomic Norm Minimization

IEEE SIGNAL PROCESSING LETTERS, VOL. 29, 2022

目录

主要内容

具体模型实现

MATLAB算法如下：

---

主要内容

相干源互质阵列DOA估计。首先，通过互质阵列插值生成一个虚拟的均匀线阵（ULA）。随后，推导出一个增强的无噪声协方差矩阵，通过求解MMV ANM来恢复Hermitian Toeplitz矩阵。最后， MUSIC算法进行DOA估计。

优点：与其他算法相比，该算法对源之间的相位差不敏感。 

具体模型实现

互质阵结构：

观测信号结构：

注解：S为互质阵的实际阵元结构。不妨设M=3，N=5，则有S={0，3，6，9，12，5，10}。Xs则是由S得到的接收信号。

虚拟阵元结构：

 虚拟阵元观测信号：

 注解：V为虚拟阵列结构。同样设M=3，N=5，则有V={0，1，2，3，4，5，6，7，8，9，10，11，12}。

协方差矩阵为：

 选择Rv中的非零列：

 注解：设M=3，N=5，则Rv是一个13x13维的矩阵，R则是由Rv的第{1，4，7，10，13，6，11}=S+1列组成的一个新的13x7维的矩阵。

 Xv可初始化如下：

则协方差矩阵为：

 R则为：

R可以表示为：

 经过一系列推导证明（证明过程略，可参考原文）：DOA的估计问题等效于求R0的原子0范数。

 上式是一个NP难题，对其凸松弛：

 由于插值，R中具有零项，因此需要恢复R中缺少的元素。

可由以下ANM问题恢复期望的R0：

 G为选择矩阵，R中元素为0时，对应位置的G元素为0，非0时为1。

该ANM问题可转换为SDP问题：

 通过CVX可对其求解。

最后，对T应用music算法即可实现对相关信号的DOA估计。

MATLAB算法如下：

%% DOA Estimation of Coherent Sources Using Coprime Array via Atomic Norm Minimization
%% IEEE SIGNAL PROCESSING LETTERS, VOL. 29, 2022


clc;
clear all;
close all;


M = 7;
N = 5;
L = M+N-1;
array1 = 0:M:(N-1)*M;
array2 = 0:N:(M-1)*N;
array = [array1 array2];
array = unique(array);
arraymax = max(array);

d = 0.5;
theta = [ -1 1 ];
K = length(theta);
A = exp(-1i*2*pi*d*array'*sind(theta));
snap = 500;
% s = randn(length(theta),snap);
s = complex(rand(length(theta),1),rand(length(theta),1));
phi = rand(1,snap) + 0.0;
s = s*exp(-1i*2*pi*phi);

% t=1:snap;
% f0 = 0.5; 
% s = 2.*(ones(K,1)*exp(j*2*pi*(f0*t)));%独立信源

x= A*s;
snr = 20;
x = awgn (x,snr);

%% x重排为xv
xv = zeros(arraymax+1,snap);
for ii = 0:arraymax
    count = find(ii == array);
    if count
        xv(ii+1,:) = x(count,:);
    end
end

Rv = xv*xv'/snap;
R = Rv(:,array+1);
% R = Rv;

G = ones(size(R));
[MM,NN]=size(R);

pos = find(R==0);
G(pos) = 0;

%cvx
%%%%%利用CVX工具箱求解凸优化问题%%%%%
mu = 1;
cvx_begin sdp quiet
% cvx_precision high
cvx_solver sdpt3
    variable T(MM,MM) hermitian toeplitz semidefinite
    variable W(NN, NN) hermitian 
    variable Z(MM, NN) complex
    % minimize(  square_pos(norm(Z .* G - R,'fro')) + mu/(2*sqrt(arraymax)) * (trace(T)+trace(W)) )	%目标函数
    minimize(  0.5*sum_square_abs(vec(Z .* G - R )) + mu/(2*sqrt(arraymax+1)) * (trace(T)+trace(W)) )	%目标函数
    [W Z'; Z T] >= 0;
cvx_end

 
derad = pi/180;
[EV,Dv] = eig(T);%特征值分解
DD = diag(Dv);%将特征值变为向量形式
[DD,I] = sort(DD);%从小到大
DD = fliplr(DD');%翻转函数，从大到小
EV = fliplr(EV(:,I));
En = EV(:,K+1:end);%噪声子空间
dm_ss = 0:arraymax;
dm_ss = dm_ss*d;
for ii = 1:2001
    angle(ii) = (ii-1001)*90/1000;
    phim = derad*angle(ii);
    a = exp(-1j*2*pi*dm_ss*sin(phim) ).';
    Pmusic(ii) = 1/(a'*En*En'*a);
end
Pmusic = abs(Pmusic);
Pmax = max(Pmusic);
Pmusic_db = 10*log10(Pmusic/Pmax);

 
plot(angle,Pmusic_db);
hold on;
plot([theta(1),theta(1)],ylim,'m-.');
plot([theta(2),theta(2)],ylim,'m-.');

 

 仿真结果：-1和30

 而且文章作者能够实现-1度和1度的分离：

-1和1时实验结果：

 

 

请各位大佬们提点提点。