【阵列信号处理】DOA估计算法

DOA估计中的ESPRIT算法

ESPRIT算法时一种利用子空间旋转法估计DOA参数的方法，其算法的基本思想是将阵列在结构上分成两个完全一致的子列，两个子列相应阵元偏移的距离相等，也就是说阵列的阵元被分成一对对的形式，而且每一对之间具有相同的平移距离，这样入射角在两个子阵列上仅相差一个旋转不变因子，该因子包含了各个入射信号的波达方向信息，因此通过求解一个广义的特征值，就可以得到入射信号的DOA。该算法与MUSIC算法相比具有以下优点：

*无需精确知道阵列流形向量，仅需要各子阵列之间保持一致，因此降低了对阵列较准的严格性
*不需要在整个空间上进行DOA谱峰搜索从而极大地降低了计算量和存储量

信号建模

从上式我们可以看出两个子列阵列流形之间的关系，将这一关系带入子阵2的接收信号中可以得到下式

ESPRIT算法的关键就在于旋转因子的求解。

主要步骤分为

旋转因子、建立方程
使用特征值分解工具：1)求非零广义特征值; 2)求伪逆; 3)LS-Esprit; 4)TLS-Esprit

LS-Esprit方法步骤

 由阵列数据估计相关矩阵
 对相关矩阵进行特征值分解
 用P个大特征值对应的特征矢量构成子空间，并将其分成E_X和E_Y两部分
 计算矩阵inv(E_X)*E_Y的特征值ξ
 计算到达角

仿真结果

仿真程序

%LS_ESPRIT ALGORITHM
clear all;
%close all;
clc;
source_number=2;%信元数
sensor_number=8;%原阵元数
sub_sensor_number=7;%子阵元数
N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=[pi/4 pi/6].';%信号频率
l=sum(2*pi*3e8./w)/length(source_number);%信号波长  
d=0.5*l;%阵元间距

snr=0;%信噪比/dB
source_doa=[-10 60];%两个信号的入射角度
A=[exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(1)*pi/180)/l);exp(-1j*(0:sensor_number-1)*d*2*pi*sin(source_doa(2)*pi/180)/l)].';%阵列流型

s=sqrt(10.^(snr/10))*exp(1j*w*[0:N_x-1]);%仿真信号
%x=A*s+(1/sqrt(2))*(randn(sensor_number,N_x)+1j*randn(sensor_number,N_x));%加了高斯白噪声后的阵列接收信号
x=awgn(A*s,10);
x1=x(1:sub_sensor_number,:);%子阵1接受的数据矢量
x2=x(2:(sub_sensor_number+1),:);%子阵2接受的数据矢量

%对两个子阵的模型进行合并
X=[x1;x2];
R=X*X'/snapshot_number;
%对R进行奇异值分解
[U,S,~]=svd(R);
%选择大的特征值构成空间并分成两部分
Us=U(:,1:source_number);
disp(Us);
Us1=Us(1:sub_sensor_number,:);
Us2=Us((sub_sensor_number+1):2*sub_sensor_number,:);
%按照公式得到旋转不变矩阵M
E=pinv(Us1)*Us2;
disp('E');
disp(E);
%对得到的旋转不变矩阵进行特征分解
[V,D]=eig(E);
disp('D');
disp(D);
D=(diag(D)).';
doa=-asin(angle(D)/pi)*180/pi;
doa=sort(doa);
polarplot(doa(1)*pi/180,1,'*',doa(2)*pi/180,1,'square');
grid on;
title('DOA Estimation by LS ESPRIT Algorithm with Monte Carlo');
hold on;
drawnow;

比较简单，从代码中我们可以看出算法中并未涉及阵列流形的计算，所以此方法将适用于阵列流形未知的情况。

TLS-Esprit方法步骤

由阵列数据估计相关矩阵
对相关矩阵进行特征值分解
用P个大特征值对应的特征矢量构成子空间，并将其分成E_X和E_Y两部分
计算矩阵(E_X)’*E_Y的特征向量，并将其分解为PXP的子阵[V_11 V_12; V_21 V_22]
计算V_12*inv(V_22)的特征值
计算到达角

与上述LS-Esprit方法多了几步而已，下面是仿真结果

下面仅给出仿真算法部分

%对两个子阵的模型进行合并
X=[x1;x2];
R=X*X'/snapshot_number;

%对R进行奇异值分解
[U,S,~]=svd(R);
Us=U(:,1:source_number);
disp(Us);
Us1=Us(1:sub_sensor_number,:);
Us2=Us((sub_sensor_number+1):2*sub_sensor_number,:);

%形成矩阵Us12
Us12=[Us1,Us2];
%对“Us12'*Us12”进行特征分解，得到矩阵F
[F,Sa,Va]=svd(Us12'*Us12);
disp('F');
disp(F);
disp(Sa);
%将F分解为四个小矩阵
F11=F(1:2,1:2);
F12=F(1:2,3:4);
F21=F(3:4,1:2);
F22=F(3:4,3:4);
%按照公式得到旋转不变矩阵E
E=-(F12*(inv(F22)));
disp('E');
disp(E);
%对得到的旋转不变矩阵进行特征分解
[V,D]=eig(E);
disp(D);
D=(diag(D)).';
doa=-asin(angle(D)/pi)*180/pi;
doa=sort(doa);

LS与TLS的Esprit算法对比
针对单一信号源的DOA估计，根据快拍数以及信噪比的不同分别进行1000次蒙特卡洛仿真，后求取估计的RMSE的平均值，实验代码较为简单这里不再给出了，下面直接给出结果。

没太看出来有多大的差别。。。

总结

本文主要实现了基础的EESPRIT算法仿真，包括LS跟TLS，并对两种方法进行了蒙特卡洛仿真，从结果上来看，ESPRIT算法由于不需要角度搜索从而大大提升了算法的运行效率；LS与TLS两者效果差不多，如果后续有些新的理解再来进行修改。
本文未给出一些详细的公式推导，如果后续有时间，将对这一部分进行补充。