Direction-of-Arrival Estimation of Coherent Signals via Coprime Array Interpolation

基于互质阵列的相干信源DOA估计（NNM）

相干信号模型：
$$\begin{gathered} {\mathbf{x}}_{\mathbf{S}}(t)=s_1(t)\sum _{k=1}^K{\alpha }_k{\mathbf{a}}_{\mathbf{S}}({\theta }_k)+{\mathbf{n}}_{\mathbf{S}}(t) \\ ={\mathbf{A}}_{\mathbf{S}}\alpha s_1(t)+{\mathbf{n}}_{\mathbf{S}}(t) \end{gathered}$$
对互质阵列进行内插：
U = { l ∣ 0 ≤ l ≤ max ⁡ ( S ) } . \mathbb{U}=\{l|0\leq l\leq\max(\mathbb{S})\}. U={l∣0≤l≤max(S)}.
$$\begin{array}{c}{\mathbf{x}}_{\mathrm{U}}(t)=s_1(t)\sum _{k=1}^K{\alpha }_k{\mathbf{a}}_{\mathbf{U}}({\theta }_k)+{\mathbf{n}}_{\mathbf{U}}(t)\\ ={\mathbf{A}}_{\mathbf{U}}\mathbf{\alpha }{s}_1(t)+{\mathbf{n}}_{\mathbf{U}}(t)\end{array}$$
R x U = E { x U ( t ) x U H ( t ) } = p 1 A U α α H A U H + σ n 2 I ∣ U ∣ \mathbf{R}_{\mathbf{x}_{\mathbf{U}}}=E\{\mathbf{x}_{\mathbf{U}}(t)\mathbf{x}_{\mathbf{U}}^{H}(t)\}=p_{1}\mathbf{A}_{\mathbf{U}}\mathbf{\alpha}\mathbf{\alpha}^{H}\mathbf{A}_{\mathbf{U}}^{H}+\sigma_{n}^{2}\mathbf{I}_{|\mathbf{U}|} RxU​​=E{xU​(t)xUH​(t)}=p1​AU​ααHAUH​+σn2​I∣U∣​
从任意一行$R_{xU}$，我们得到如下的Toeplitz矩阵:
$$\mathbf{F}(m)=\left[\begin{array}{cccc}{\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,L-1}& {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,L}& \cdots & {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,|\mathbf{U}|-1}\\ {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,L-2}& {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,L-1}& \cdots & {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,|\mathbf{U}|-2}\\ ⋮& ⋮& \ddots & ⋮\\ {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,0}& {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,1}& \cdots & {\left[{\mathbf{R}}_{{\mathbf{x}}_{\mathbf{U}}}\right]}_{m,L-1}\end{array}\right]$$
通过求解核范数最小化问题恢复托普利兹矩阵：
$$\begin{array}{rl}& \underset{\tilde{\mathbf{F}}(m)\in {\mathbb{C}}^{L\times L}}{\min }\Vert \tilde{\mathbf{F}}(m){\Vert }_{\ast }\\ & \mathrm{subject}\mathrm{to}\tilde{\mathbf{F}}(m)\circ \mathbf{C}=\hat{\mathbf{F}}(m)\end{array}$$
仿真代码：

%% 《Direction-of-Arrival Estimation of Coherent Signalsvia Coprime Array Interpolation》仿真
%% IEEE SIGNAL PROCESSING LETTERS, VOL. 27, 2020
%% 相干信源DOA估计
clc;
clear all;
close all;
M = 3;
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 = [10 30 ];
K = length(theta);
snap = 500;
f0 = 1e4;%载频
fs = 100*f0;%采样频率
d = 0.5;
snr =20;
t = [1:snap]/fs;
s = exp(1j*2*pi*(f0.*t));
A = exp(-1i*pi*array'*sin(theta/180*pi));
alpha = [1;exp(1j*pi/4)];
x= A*alpha*s;
x = awgn (x,snr,'measured');
%% 阵列内插
x_v = zeros(arraymax+1,snap);
for i = 0:arraymax
    count = find(i == array);
    if count
        x_v(i+1,:) = x(count,:);
    end
end
l =  (length(0:arraymax)+1)/2;
Rv = (x_v*x_v')/snap;
zD = Rv(1,:);
%% 重构Teoplitz矩阵
% F = zeros(l,l);%重构的协方差矩阵
F = toeplitz(fliplr(zD(1:l)),zD(l:end));
G = ones(l,l);
posi0V = find(F==0);
G(posi0V) = 0;
%% 求解凸优化问题
mu = 0.25;
cvx_begin quiet
variable T(l, l)  complex toeplitz 
minimize(norm_nuc(T) )
subject to
T.* G==F;
cvx_end
%% musicDOA估计
array = 0:l-1;
[Pmusic_db,angle] = MUSIC(T,K,array);
plot(angle,Pmusic_db,'Linewidth',1);
hold on
for n = 1:K
    plot([theta(n),theta(n)],ylim,'r-.','Linewidth',1);
end
xlabel('角度(°)');
ylabel('归一化幅度');
legend('NNM','波达方向');
[pks,locs] = findpeaks(Pmusic_db);
[a,b] = sort(pks,'descend');
angle_mer=sort(angle(locs(b(1:K))));
disp(angle_mer);
% % [U,D] = eig(T);                                       % 特征值分解
% % [D,I] = sort(diag(D));                                % 将特征值排序从小到大
% % U = fliplr(U(:, I));                                  % 对应特征矢量排序，fliplr 之后，较大特征值对应的特征矢量在前面
% % Un = U(:, K+1:end);                                     % 噪声子空间
% % Gn = Un*Un';
% % % 提取多项式系数并对多项式求根
% % len = arraymax+1;
% % coe = zeros(1, 2*len-1);
% % for i = -(len-1):(len-1)
% %     coe(-i+len) = sum(diag(Gn,i));
% % end
% % r = roots(coe);                                       % 利用roots函数求多项式的根
% % r = r(abs(r)<1);                                      % 找出在单位圆里的根
% % [lamda,I] = sort(abs(abs(r)-1));                      % 挑选出最接近单位圆的K个根
% % Theta = r(I(1:K));
% % Theta = asin(-angle(Theta)/pi)*(180/pi);                 % 计算信号到达方向角
% % Theta = sort(Theta).';
% % disp('估计结果：');
% % disp(Theta);

函数MUSIC代码如下：

function [Pmusic_db,angle] = MUSIC(R,K,Array)

[EV,Dv] = eig(R);%特征值分解
% [EV,Dv] = eig(T);%特征值分解
DD = diag(Dv);%将特征值变为向量形式
[DD,I] = sort(DD);%从小到大
DD = fliplr(DD');%翻转函数，从大到小
EV = fliplr(EV(:,I));
En = EV(:,K+1:end);%噪声子空间
angle=-90:0.01:90;
for m=1:length(angle)
    a=exp(-1i*pi*Array'*sin(angle(m)*pi/180));
    Pmusic(m)=1/(a'*En*En'*a);
end
Pmusic = abs(Pmusic);
Pmax = max(Pmusic);
Pmusic_db =pow2db(Pmusic/Pmax);
end

仿真结果如下图：