写在最前

本文章的代码已开源，基于凸优化法（CVX）实现DOA时，依赖CVX工具箱，如果你的MATLAB没有安装，请前往这里下载，解压后在MATLAB命令行，cd到解压目录并执行其中的“cvx_setup.m”文件进行安装。如果不做CVX部分的仿真，可忽略这一步。

前言

波达方向估计(Direction Of Arrival, DOA)也称为测向、空间谱估计，为利用电磁波来获取目标或信源对天线阵列的角度信息，主要应用于雷达、通信、电子侦察与对抗等领域。
本文利用MIMO天线阵列实现DOA相关算法的总结，主要仿真实现了常规波束形成(CBF)、Capon和最大似然估计(ML)三种常规方法，多重信号分类法(MUSIC)、LS-ESPRIT和TLS-ESPRIT三种子空间方法，欠定系统聚焦法(FOCUSS)、正交匹配追踪法(OMP)、凸优化法(CVX)、伪逆法(PINV)和期望最大化-稀疏贝叶斯学习法(EM-SBL)等稀疏恢复方法。对比了上述方法在常规、低信噪比、低快拍以及信源相干情况下的性能，并研究了空间平滑算法在处理相干信源问题上的表现。

空间谱估计的历史发展

本文不涉及深度学习法。

仿真原理

离散时间阵列信号模型

天线阵元示意图如下图所示，在仿真中，默认最左边的阵元为相位中心，仿真的前提为信源发射的信号为窄带信号且为远场信号，阵元为全向天线其阵元之间没有互耦和误差。

为了方便数字化处理，仿真采用了离散时间阵列信号模型：
$$\begin{gathered} \left[\begin{array}{c}{x}_1(n)\\ x_2(n)\\ ⋮\\ x_M(n)\end{array}\right]=\left[\begin{array}{ccccc}& 1& 1& \cdots & 1\\ & e^{-j2\pi d\sin {\theta }_1/\lambda }& e^{-j2\pi d\sin {\theta }_2/\lambda }& \cdots & e^{-j2\pi d\sin {\theta }_p/\lambda }\\ & ⋮& ⋮& \ddots & ⋮\\ & e^{-j2\pi (M-1)d\sin {\theta }_1/\lambda }& e^{-j2\pi (M-1)d\sin {\theta }_2/\lambda }& \cdots & e^{-j2\pi (M-1)d\sin {\theta }_p/\lambda }\end{array}\right]\left[\begin{array}{c}{s}_1(n)\\ s_2(n)\\ ⋮\\ s_p(n)\end{array}\right]+\left[\begin{array}{c}{n}_1(n)\\ n_2(n)\\ ⋮\\ n_M(n)\end{array}\right] \\ \mathbf{x}=\mathbf{A}(\theta )\mathbf{s}+\mathbf{n} \end{gathered}$$
其中，$\mathbf{x}\in {\mathbb{C}}^{M\times N}$为天线阵列接收到的信号；$\mathbf{A}\in {\mathbb{C}}^{M\times P}$为空域导向矢量；$\mathbf{s}\in {\mathbb{C}}^{p\times N}$为目标发出的信号；$\mathbf{n}\in {\mathbb{C}}^{M\times N}$为噪声；$d$为阵元间隔；$\lambda$为载频信号的波长；$N$为快拍数；$M$为阵元数量；$P$为目标数量。

波束形成矩阵(完备字典)

该矩阵也称为DOA栅格矩阵，用于估算在每个角度上的功率，该矩阵的表达式为：
$$a({\theta }_0,{\theta }_1,\dots ,{\theta }_k)=\left[\begin{array}{ccccc}& 1& 1& \cdots & 1\\ & e^{j2\pi d\sin {\theta }_0/\lambda }& e^{j2\pi d\sin {\theta }_1/\lambda }& \cdots & e^{j2\pi d\sin {\theta }_k/\lambda }\\ & e^{j2\pi (2d)\sin {\theta }_0/\lambda }& e^{j2\pi (2d)\sin {\theta }_1/\lambda }& \cdots & e^{j2\pi (2d)\sin {\theta }_k/\lambda }\\ & ⋮& ⋮& \ddots & ⋮\\ & e^{j2\pi (M-1)d\sin {\theta }_0/\lambda }& e^{j2\pi (M-1)d\sin {\theta }_1/\lambda }& \cdots & e^{j2\pi (M-1)d\sin {\theta }_k/\lambda }\end{array}\right]$$
其中，${\theta }_0,{\theta }_1,\dots ,{\theta }_k$为需要做功率估计的等间距的来向角度；$a\in {\mathbb{C}}^{M\times (k+1)}$为波束形成矩阵。

回波生成

为了保证信号的正交性，省去上变频步骤，本仿真使用随机数+带通滤波方法实现指定带宽与中心频率的信源生成。针对每个信源在$[f_0-bw,f_0+bw]$之间等间隔地划分k个频段，用高斯分布的随机数初始化矩阵$s$后，对矩阵的每一行执行对应频段的带通滤波以获得频域正交的信号，是不是有种OFDM的感觉了。

空间平滑滤波

在实际使用中，由于存在多径传输等影响因素。阵列接收到的信源信号往往是相干而非完全独立的，为了解决相干源的DOA估计问题，需要使用空间平滑滤波破坏掉接收信号的相干性。该方法将等数量的相邻的实际阵列构成子阵列，并利用所有子阵列的自相关矩阵的平均数作为阵列接收信号的自相关矩阵估计，从而解决实际自相关矩阵的相干问题。当子阵大小L越大时，去相关能力越强，但等效阵列数会越小，从而降低空间分辨率。该方法可表示为：
R x x ^ = E { R i } \hat {R_{xx}}=E\{ R_i \} Rxx​^​=E{Ri​}
其中$R_i$为第$i$个子阵列的自相关矩阵。

传统方法

CBF~常规波束成型

该方法大约在二战期间被提出，其本质是傅里叶变换在空域的直接应用，其分辨率受限于瑞丽限。该方法可表示为：
$$J_{CBF}(\theta )=\frac{a^H(\theta )R_{xx}a(\theta )}{|a^H(\theta )a(\theta ){|}^2}$$

Capon~最小方差无失真响应法

这是一种自适应方法，提出于60年代，它将维纳滤波的思想应用于空域处理，相比于CBF法，分辨率得到了一定的提高。Capon法可表示为：
$$J_{Capon}(\theta )=\frac{1}{a^H(\theta )R_{xx}^{-1}a(\theta )}$$

ML~最大似然估计法

最大似然估计在九十年代提出，实质为给定观测数据，评估模型参数。在DOA估计中，ML法可表示为：
$$J_{ML}(\theta )=\mathrm{trac}(\frac{a(\theta )a^H(\theta )}{(a^H(\theta )a(\theta ){)}^{-1}}\ast R_{xx})/M$$
在只有一个目标信源时，该方法等效于CBF法，在多个信源时，ML法为多目标优化问题，运算量较大。

子空间方法

MUSIC~多重信号分类法

时间回到上世纪八十年代，多重信号分类法在这个时候被提出。相比于CBF和Capon法仅使用接收数据的自相关矩阵直接处理，该方法考虑了信号和噪声的分布特性，通过对自相关矩阵做特征分解得到信号子空间和噪声子空间，利用两个子空间的正交特性进行DOA估计。
由于接收信号的自相关矩阵为赫米特(Hermite)矩阵，满足$R^H=R$，根据相关引理，$R$酉相似于对角矩阵，可进行酉对角化分解，即存在 酉矩阵$U$和对角阵$V$，使$U^{H}RU=V$，且V主对角线上的元素为$R$的特征值。
MUSIC方法可表示为：
$$\begin{gathered} J_{MUSIC}(\theta )=\frac{a^H(\theta )a(\theta )}{a^H(\theta )(I-U_s{U}_s^H)a(\theta )}或 \\ {J}_{MUSIC}(\theta )=\frac{a^H(\theta )a(\theta )}{a^H(\theta )(U_n{U}_n^H)a(\theta )} \end{gathered}$$
其中，$U_s\in {\mathbb{C}}^{M\times k}$为信号子空间，为$U$中$k$个最大特征值所对应的特征列向量组成的矩阵；$U_n\in {\mathbb{C}}^{M\times (M-k)}$为噪声子空间；$k$表示信源个数。在应用中可根据信源数量选择信号子空间解法或噪声子空间解法。

ESPRIT~旋转不变子空间法

与MUSIC法在同一时期被提出的子空间方法还有旋转不变子空间法，区别于MUSIC使用酉对角化后信号和噪声子空间的正交性，ESPRIT利用信号子空间的旋转不变特性来求解DOA。
假设存在两个完全相同的子阵列，且子阵列的间距$\Delta$已知，假设两个子阵列接收到数据分别为$x_1$和$x_2$，则：
$$\begin{gathered} x_1=As+n_1 \\ {x}_2=A\Phi s+n_2 \\ \Phi =\mathrm{diag}(e^{j{\varphi }_1}{e}^{j{\varphi }_2}\cdots e^{j{\varphi }_p}) \end{gathered}$$
只要得到两个子阵列的旋转不变关系$\Phi$，就能得到信号的到达角，完成DOA估计。
仿真时，令$x_1$和$x_2$分别为接收到的信号$\mathbf{x}$的前$M-1$行和后$M-1$行所构成的两个子矩阵，令回波子阵列合并为$x_{exp}=[x_1;x_2]$，计算$x_{exp}$的自相关矩阵 R e s p = E { x e x p ∗ x e x p H } R_{esp}=E\{x_{exp}*x_{exp}^H\} Resp​=E{xexp​∗xexpH​}，对$R_{esp}$做酉对角化分解，得到酉矩阵$U_{esp}$和特征值对角矩阵$V_{exp}$。根据目标数量从$U_{esp}$中提取信号子空间并进行拆分得到$U_s$的前M-1行$U_{s1}$和后M-1行$U_{s2}$，接下来需要求解ESPRIT矩阵$\Xi$，使：
$$U_{s1}\Xi =U_{s2}$$
求解出ESPRIT矩阵$\Xi$后，通过求取该矩阵的特征值并求其复角度即可得到来向角的正弦值，即：
$$DO{A}_{exp}=\arcsin (\frac{-\mathrm{angle}(V_{esp})\ast \lambda }{2\pi d})$$
其中，$V_{esp}\in {\mathbb{C}}^{P\times 1}$是由$\Xi$的特征值构成的向量。
求解ESPRIT矩阵$\Xi$可以使用最小二乘准则或总体最小二乘准则。

最小二乘准则

不难证明，在最小二乘准则下，使用下式对ESPRIT矩阵${\Xi }_{LS}$进行求解：
$${\Xi }_{LS}=\mathrm{pinv}(U_{s1})\ast U_{s2}$$

总体最小二乘准则

在总体最小二乘准则下，令$V_{TLS-esp}=[U_{s1},U_{s2}]$并对V做奇异值分解得到右奇异矩阵$V_{TLS-ESP}\in {\mathbb{C}}^{(P\ast 2)\times (P\ast 2)}$，将$V_{TLS-ESP}$等分为4个方阵并取右上角方阵和右下角方阵分别记为$E_{12}$和$E_{22}$，即可解出ESPRIT矩阵${\Xi }_{TLS}$：
$${\Xi }_{TLS}=-E_{12}\ast \mathrm{inv}(E_{22})$$

稀疏恢复方法

_{终于写到我最感兴趣的部分了。}
2004年，陶哲轩、Emmanuel Candes和David Donoho等人证明了“如果信号是稀疏的，那么它可以由远低于采样定理要求的采样点恢复信号”，紧接着提出了“压缩感知”概念。总得来讲，稀疏恢复的前提条件之一为信号的稀疏性，即信号在某个域（常用作信号处理的领域包括时域、频域、多普勒域和空域等）上是稀疏的，在该域上零点的数量远远大于非零点的数量；前提条件之二为观测矩阵和稀疏表示基之间的不相关。如果一个信号在某个变换域是稀疏的，那么可以用一个与变换基不相关观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题可从这些少量投影中以高概率重构出信号。
压缩感知的步骤包括压缩采样和信号重构：
压缩采样：
$$y=\Phi \Psi s=\Theta s$$
$\Phi$和$\Psi$分别为观测矩阵和稀疏表示基，它们是不相关的；$\Theta$为传感矩阵
信号重构：
$$\underset{s}{\min }\Vert s{\Vert }_{0}s.t.y=\Theta s$$
其中$\Vert \cdot {\Vert }_0$为$p_0$范数，表示$x$中非零元素的个数，在实际求解中，一般将该约束松弛为$p_1$或$p_2$等范数。
对于方程：$y=\Theta s$，已知$y$和$\Theta$求解未知数$x$时，如果超定方程，问题实质为拟合问题；如果$\Theta$为满秩的方阵，问题则为线性方程组求解；当$\Theta$为欠定时，才是稀疏恢复问题。
在仿真中，对接收信号的自相关矩阵$R_{xx}$做酉对角化后，取酉矩阵中最大特征值对应的特征向量作为稀疏恢复的观测值矩阵y，将波束形成矩阵$a$作为传感器矩阵，稀疏恢复得到的$s$即为DOA估计结果，即求解方程：
$$y=as$$

FOCUSS~欠定系统聚焦法

在该方法中，首先确定$s_0$的迭代初始值为：
$$s_0=a^H\mathrm{inv}(a{a}^H)y$$
迭代过程为：
$$\begin{gathered} W=\mathrm{diag}(s_k^{1-{\lambda }_{spe}/2}) \\ {s}_{k+1}=W{W}^H{a}^H\ast \mathrm{inv}(aW{W}^H{a}^H+{\lambda }_{reg}I)\ast y \end{gathered}$$
其中，${\lambda }_{reg}$为正则化因子，过大会趋于0解，过小结果发散；${\lambda }_{spe}$为稀疏因子，效果等效于稀疏解的范数约束。
迭代次数超过最大给定次数，或误差小于误差限制${\lambda }_{err}$（即$\frac{\Vert s_{k+1}-s_k{|}_2}{\Vert s_k\Vert }<{\lambda }_{err}$）时，结束迭代。

OMP~正交匹配追踪法

该方法的迭代步骤如：

1. 初始化，令标签集${\Omega }_{mask}$为空集并初始化残差向量$r_0=y$。
2. 辨识，求取传感器矩阵$a$中与残差向量$r_k$最强相关的原子，即$j_k\in \arg max|<{\mathrm{r}}_{\mathrm{k}},{\varphi }_{\mathrm{j}}>|$，往${\Omega }_{mask}$中添加该原子的索引。
3. 估计，最小化问题$mi{n}_x\Vert y-{\Phi }_{{\Omega }_k}x{\Vert }_2^2$的最小二乘解为$x_k=\mathrm{inv}({\Phi }_{{\Omega }_k}^H{\Phi }_{{\Omega }_k})\ast {\Phi }_{{\Omega }_k}^{H}y$。
4. 更新残差，$r_{k+1}=[I-{\Phi }_{{\Omega }_k}({\Phi }_{{\Omega }_k}^T{\Phi }_{{\Omega }_k}{)}^{-1}{\Phi }_{{\Omega }_k}^T]y$
5. 重复步骤2~4直到满足迭代停止条件，退出条件包括迭代次数达到字典的栅格数，或者$a(:,{\Omega }_{mask}{)}^{\prime }\ast a(:,{\Omega }_{mask})$为奇异矩阵。
6. 输出结果，
   $$x(i)=\{\begin{array}{cccc}& x_k(i)& ,& i\in {\Omega }_{mask}\\ & 0& ,& 其它\end{array}$$

伪逆法

对于方程：$y=as$，令解s为：
$$s=y^H\ast {\mathrm{pinv}(a)}^H$$
该方法因为运算量较大且产生的不是稀疏解，在实际工程中很少使用，在仿真中仅用于效果对比。

EM-SBL~最大期望-稀疏贝叶斯学习法

SBL算法是稀疏信号重构的重要算法，涉及高斯分布、最大似然估计、贝叶斯公式等数学理论，由于推导过程比较复杂，_{主要是我懒得敲了，后续如果读者有强烈需要我再更新一下吧^_^}，这里直接介绍推导结果。该方法迭代的步骤为：

1. 初始化，令${\Gamma }_0=0.1\ast I$。
2. E-step，令${\sigma }_x=\mathrm{pinv}(\sigma (a^{H}a+\mathrm{pinv}({\Gamma }_k))$，令${\mu }_x={\sigma }_x{\sigma }^{-1}{a}^H{R}_{xx}$，其中$\sigma$表示阵列接收信号的估计噪声方差，$R_{xx}$为接收信号的自相关矩阵。
3. M-step，令${\Gamma }_{k+1}={\mu }_x(\theta )\ast {\mu }_x(\theta {)}^H/N+{\sigma }_x(\theta )$
4. 重复步骤2和3，直到达到最大迭代次数或$\Vert {\Gamma }_{k+1}-{\Gamma }_k{\Vert }_1<误差限$

CVX~凸优化法

由于方程“$y=as$”的求解可以转化为凸优化问题：
$$argmi{n}_s\Vert s{\Vert }_p,s.t.\Vert y-as{\Vert }_2<{\lambda }_{lim}$$
因此可以使用凸优化工具箱进行求解，并且还可以根据需要调整用于约束稀疏解$s$的范数，当使用$p_1$范数时，上述优化将产生稀疏解，当范数$p$增大时，相当于增加了松弛量，解密度会越来越大，从而变得越来越不稀疏。
${\lambda }_{lim}$为收敛控制量，过小可能产生发散的结果，过大则可能产生非稀疏解。

仿真实现

代码下载

代码已上传到GitHub：https://github.com/highskyno1/MIMO_DOA

可修改的参数

1. 载波频率$f_0$：信源发射信号的中心频率。
2. 系统全局采样率$f_g$：由于仿真是数字系统，该采样率为仿真所使用的全局采样率。
3. 信号带宽$BW$：信源发射信号的总带宽，当存在多个信源时，如果使用非相干生成模式，会根据总带宽为每个信源划分不重合的等长信号频段；否则所有信源使用整个带宽的频段。
4. 目标来向角度$sr{c}_{angle}$：目标真实的相对于阵元平面法线的入射角度，内部元素数量被计算为目标数。考虑到实际的MIMIO阵列的扫描角度一般只有120°，取值范围为[-60°，60°]，如果要在这个范围外取值，请连同$sca{n}_{angle}$一起修改。
5. 阵元数量$M$：MIMO阵列阵元数量。
6. 快拍数$L$：每个阵元的采样点个数。
7. 信噪比snr：信源发出的信号被收到时的信噪比。
8. 扫描角度$sca{n}_{angle}$：控制波束形成矩阵的角度，也称为完备字典或DOA栅格，可以改变间隔以提高空域采样率，注意空域分辨率不因空域采样率的改变而改变。
9. 是否使用空间平滑滤波is_space_smooth：逻辑变量，置为ture时使能空间平滑滤波以解决信源的相干性问题。
10. 是否产生相干的信源is_coherent：逻辑变量，置为true时将产生频域重叠的相干信源。

仿真结果

常规情况

在信源不相干，信噪比达到20dB，阵元数和快拍数分别为32和128的情况下，DOA结果如下图所示，可见CBF法主瓣最宽、旁瓣最高，其次为EM-SBL方法。CBF、EM-SBL、Capon、MUSIC的主瓣宽度依次减少。在稀疏恢复类型的算法中，FOCUSS算法效果最佳，OMP算法产生了过多的假峰，PINV法和ML法则出现了与传统方法相当的主瓣宽度。CVX方法则没有得到正确的来向角估计。两种准则下的ESPRIT算法均估计出了正确的来向角。
尽管CVX算法的结果不太理想，但仍可以发现在p1范数约束下将产生较为稀疏的解，随着范数的增加，稀疏恢复结果变得越来越不稀疏。

低信噪比情况

修改信噪比为0dB，仿真结果如下图所示，可见此时常规方法和MUSIC法仍然可以正常工作；ESPRIT算法估计的DOA角度出现偏差；稀疏恢复部分，FOCUSS方法仍然表现出良好的性能，OMP法的假峰消失，CVX方法的零陷刚好落在真实来波方向上。

修改信噪比为-10dB，仿真结果如下图所示，此时大部分算法都只能找到其中的一个来向角。

相干信源情况

信源发射信号相干且不使用平滑滤波时，仿真结果如下图所示，可见子空间类型算法和Capon法只能发现一个来向的角度，但稀疏恢复算法和CBF

信源发射信号相干且使用平滑滤波时，设置等效阵元数量为14，仿真结果如下图所示，此时Capon法、MUSIC法变得正常，可以分辨两个来向的信号，但ESPRIT仍然无法使用。但由于等效阵元数的降低，DOA结果的主瓣变换，降低了空域分辨率，这与理论是符合的。

少快拍情况

基于正常情况设置快拍数为32个（刚好等于阵元数）时，仿真结果如下图所示，可见此时算法仍然可以正常工作。

观察更极限的情况，设置快拍数为2，仿真结果如下图所示，可见此时CBF法仍然可以正常工作，Capon法出现了相当多的假峰，ESPRIT算法已无法正常工作。稀疏恢复类型的算法在此种情况下表现出了极佳的性能，FOCUSS算法可以准确地估计出来向角，CVX法可以正常工作但有一定误差，OMP算法在真实来向附近产生了多个假峰。

仿真小结

1. CBF算法是万精油，可以用在低信噪比、低快拍和信源相干情况下，低时间复杂度，少空间占用使得它时至今天仍然是DOA的主流算法之一，在你不知道使用什么算法时，不妨先试试CBF吧。
2. 稀疏恢复算法是概率类型的算法，每次执行都可能产生不一样的结果，这种算法优点在于可以获得很高的分辨率（超分辨率识别），并且可以在低快拍和信源相干情况下使用，是当前除了机器学习算法外的研究主流，稀疏算法唯一的缺点在于它的鲁棒性欠佳，且只能在高信噪比下使用，高信噪比也对应了空域“稀疏”。

后语

心里话：终于写完了！！！这篇长文从代码重构，到撰写与敲公式前前后后画了差不多一个星期，文中的所有公式、符号和矩阵等都是用$\mathcal{L}atex$手敲的，通过撰写这篇文章我的Latex技术又得到了精进。
这是三年后复出的第二弹，后面还有，敬请期待。
看到这里了，_点个赞再走吧，让我知道有人把文章看完也是我写下去的动力！

附录代码

仅附录DOA主要部分的代码，完整代码的下载请前往GitHub。

DOA_CBF.m

function P_cbf = DOA_CBF(scan_a, R)
% DOA_CBF 基于常规波束形成法实现DOA估计
%   scan_a  需要估计的来波方向
%   R       快拍的协方差矩阵，L*L维，L为快拍数
%   P_cbf   各个scan_a栅格上的来波归一化强度

scan_len = size(scan_a,2);
P_cbf = zeros(1,scan_len);
for i = 1:scan_len
    foo = scan_a(:,i);
    P_cbf(i) = foo' * R * foo / (foo' * foo)^2;
end
P_cbf = abs(P_cbf);
P_cbf = P_cbf ./ max(P_cbf);
end

DOA_Capon.m

function P_capon = DOA_Capon(scan_a, R)
% DOA_Capon 基于Capon法实现DOA估计
%   scan_a  需要估计的来波方向
%   R       快拍的协方差矩阵，L*L维，L为快拍数
%   P_capon 各个scan_a栅格上的来波归一化强度

scan_len = size(scan_a,2);
P_capon = zeros(1,scan_len);
for i = 1:scan_len
    foo = scan_a(:,i);
    P_capon(i) = 1 / (foo' * pinv(R) * foo);
end
P_capon = abs(P_capon);
P_capon = P_capon ./ max(P_capon);
end

DOA_MUSIC

function P_music = DOA_MUSIC(U, target_len, scan_a)
% DOA_MUSIC 基于多重信号分类法实现DOA估计
%   U           对接收信号的自相关矩阵做酉对角化分解后的左酉矩阵
%   target_len  目标数量
%   rec_len     阵元数量
%   scan_a      DOA估计栅格

% 计算阵元数和栅格数
[rec_len,scan_len] = size(scan_a);
W = fliplr(U);
% 已知只有两个信源，用信号子空间法
U_s = W(:,1:target_len);
U_music = eye(rec_len) - U_s * U_s';
P_music = zeros(1,scan_len);
for i = 1:scan_len
    foo = scan_a(:,i);
    P_music(i) = 1 / (foo' * U_music * foo);
end
P_music = abs(P_music);
P_music = P_music ./ max(P_music);
end

DOA_ESPRIT.m

function [DOA_esp_ml, DOA_esp_tls] = DOA_ESPRIT(x_sig, target_len, lamda, d)
% DOA_ESPRIT 基于旋转不变子空间法实现DOA
%   x_sig       每个阵元接收到的信号矩阵，阵元数*快拍数
%   target_len  目标数量
%   lamda       载波波长
%   d           阵元间隔
%   DOA_esp_ml  基于最大似然估计准则得到的估计结果
%   DOA_esp_tls 基于最小二乘准则得到的估计结果

% 计算阵元数
rec_len = size(x_sig,1);
% 回波子阵列合并
x_esp = [x_sig(1:rec_len-1,:);x_sig(2:rec_len,:)];
% 计算协方差
R_esp = cov(x_esp');
% 特征分解
[~,~,W] = eig(R_esp);
% 获取信号子空间
W = fliplr(W);
U_s = W(:,1:target_len);
% 拆分
U_s1 = U_s(1:rec_len-1,:);
U_s2 = U_s(rec_len:end,:);

%% LS-ESPRIT法
mat_esp_ml = pinv(U_s1) * U_s2;
% 获取对角线元素并解算来向角
DOA_esp_ml = -angle(eig(mat_esp_ml));
DOA_esp_ml = asin(DOA_esp_ml .* lamda ./ 2 ./ pi ./ d);
DOA_esp_ml = rad2deg(DOA_esp_ml);

%% TLS-ESPRIT
V = [U_s1,U_s2];
[~,~,V] = svd(V);
% 提取E12和E22
E12 = V(1:target_len,target_len+1:end);
E22 = V(target_len+1:end,target_len+1:end);
mat_esp_tls = - E12 / E22;
% 获取对角线元素并解算来向角
DOA_esp_tls = -angle(eig(mat_esp_tls));
DOA_esp_tls = asin(DOA_esp_tls .* lamda ./ 2 ./ pi ./ d);
DOA_esp_tls = rad2deg(DOA_esp_tls);
end

DOA_ML.m

function P_ml = DOA_ML(scan_a, R_sig)
% DOA_ML 基于最大似然估计法得到DOA估计
%   scan_a  DOA估计栅格矩阵
%   R_sig   接收信号的自相关矩阵，阵元数*阵元数
%   P_ml    通过ML法得到的归一化估计结果

% 计算阵元数和栅格数
[rec_len,scan_len] = size(scan_a);
P_ml = zeros(1,scan_len);
for i = 1:scan_len
    scan = scan_a(:,i);
    Pa = scan / (scan'*scan) * scan';
    P_ml(i) = trace(Pa*R_sig) / rec_len;
end
P_ml = abs(P_ml);
P_ml = P_ml ./ max(P_ml);
end

DOA_FOCUSS.m

function P_focuss = DOA_FOCUSS(scan_a, u, lamda_spe, lamda_reg, lamda_err)
% DOA_FOCUSS 基于欠定系统局灶解法(Focal Under determinedSystem Solver)
% 实现稀疏恢复获得DOA估计结果
%   scan_a      DOA估计的栅格，在稀疏恢复理论中也称为"超完备字典"
%   u           对回波自相关矩阵做酉对角化后，最大特征值对应的酉向量
%   lamda_spe   稀疏因子，效果类似于结果的范数约束
%   lamda_reg   正则化因子，过大会趋于0解，过小结果发散
%   lamda_err   迭代结束误差
%   P_focuss    通过FOCUSS法得到的归一化来波方向功率估计

% 计算阵元数
rec_len = size(u,1);
% 确定s的初始值
Dg = scan_a;
s0 = Dg' / (Dg * Dg') * u;
for i = 1:1000
    W = diag(s0.^(1-lamda_spe/2));
    s = W * W' * Dg' / (Dg * (W * W') * Dg' + lamda_reg .* eye(rec_len)) * u;
    if norm(s - s0,2) / norm(s0,2) < lamda_err
        break;
    end
    s0 = s;
end
P_focuss = abs(s);
% 数据饱和钳制
P_focuss = P_focuss ./ max(P_focuss);
P_focuss(P_focuss < 1e-4) = 1e-4;
end

DOA_OMP.m

function P_omp = DOA_OMP(u, scan_a, omp_toler)
% 基于正交匹配追踪法实现来波的DOA估计
%   u           对回波自相关矩阵做酉对角化后，最大特征值对应的酉向量
%   scan_a      DOA估计的栅格，在稀疏恢复理论中也称为"超完备字典"
%   omp_toler   容许的最小误差，低于此误差时结束迭代
%   P_omp       基于OMP法得到的不同方向来波功率估计

% 计算阵元数量和栅格数量
[rec_len,scan_len] = size(scan_a);
% 初始化标签集合
omp_omiga_mask = zeros(1,scan_len,'logical');
% 初始化初始向量
omp_r = u;
for i = 1:scan_len
    % 求字典矩阵中与残差向量rk-1最强相关的原子
    foo = abs(omp_r' * scan_a);
    [~,idx] = max(foo);
    % 往omiga中添加该列索引
    omp_omiga_mask(idx) = true;
    % 依据索引从字典中选出需要的列向量
    omp_omiga_foo = scan_a(:,omp_omiga_mask);
    % 奇异判断
    rms = det(omp_omiga_foo' * omp_omiga_foo);
    % 退出条件
    if rms < omp_toler
        break;
    end
    % 最小化min(x) ||y-fa*x|| L2范数的平方最小
    omp_xk = (omp_omiga_foo' * omp_omiga_foo) \ omp_omiga_foo' * u;
    % 更新残差
    omp_r = (eye(rec_len) - omp_omiga_foo / (omp_omiga_foo.'*omp_omiga_foo) * omp_omiga_foo.') * u;
end
P_omp = zeros(1,scan_len);
% 提取结果
foo = find(omp_omiga_mask==true);
P_omp(foo(1:length(omp_xk))) = abs(omp_xk)';
% 归一化
P_omp = P_omp ./ max(P_omp);
P_omp(P_omp < 1e-4) = 1e-4;
end

DOA_PINV.m

function s_pinv = DOA_PINV(u, scan_a)
% DOA_PINV 基于伪逆法实现稀疏恢复，效果等同于最大似然估计法
%   u       对回波自相关矩阵做酉对角化后，最大特征值对应的酉向量
%   scan_a  DOA估计的栅格，在稀疏恢复理论中也称为"超完备字典"
%   s_pinv  基于PINV法得到的不同来波方向的归一化功率

s_pinv = u' * pinv(scan_a)';
s_pinv = abs(s_pinv);
s_pinv = s_pinv ./ max(s_pinv);
end

DOA_EM_SBL.m

function s_sbl = DOA_EM_SBL(sigma, scan_a, R_sig, shot_len, err_lim, times_lim)
% 基于期望最大化-稀疏贝叶斯学习方法实现DOA估计
%   sigma       估计的噪声方差
%   scan_a      DOA估计的栅格，在稀疏恢复理论中也称为"超完备字典"
%   R_sig       回波自相关矩阵
%   shot_len    快拍数
%   err_lim     迭代误差限，迭代退出的条件之一
%   times_lim   迭代次数限制，迭代退出的条件之二

% 计算DOA栅格数量
scan_len = size(scan_a,2);
times_cnt = 0;
Gamma = eye(scan_len)*0.1;  % 初始Gamma矩阵
while true 
    times_cnt = times_cnt + 1;
    % E-step
    Sigma_x = pinv(sigma * (scan_a'*scan_a) + pinv(Gamma));
    Mu_x = Sigma_x / sigma * scan_a' * R_sig;
    % M-step
    Gamma_new = Gamma;
    for i = 1:scan_len
        mu_xn = Mu_x(i,:); 
        Gamma_new(i,i) = mu_xn*mu_xn'/shot_len + Sigma_x(i,i);
    end
    if sum(abs(diag(Gamma_new - Gamma))) < err_lim || times_cnt > times_lim
        break;
    end
    Gamma = Gamma_new;
end
Gamma_new = abs(diag(Gamma_new));
s_sbl = Gamma_new ./ max(Gamma_new);
end

DOA_CVX.m

function s_cvx = DOA_CVX(u, scan_a, p_norm, tor_lim)
% DOA_CVX 尝试利用凸优化方法实现稀疏恢复，获得DOA估计
%   !!使用前必须先安装CVX凸优化工具箱，下载地址为：http://cvxr.com/cvx/download/
%   u       对回波自相关矩阵做酉对角化后，最大特征值对应的酉向量
%   scan_a  DOA估计的栅格，在稀疏恢复理论中也称为"超完备字典"
%   s_cvx   基于CVX法得到的不同方向来波功率估计
%   p_norm  约束结果的P范数
%   tor_lim 对稀疏恢复整体(u - scan_a * s_cvx)的容许限度
%           如果结果为+inf，说明发散了，请增大该值
% 注意在CVX工具箱的scope中
% 尽量调用工具箱自己的函数（定义在CVX解压目录/functions 文件夹下）
% 使用MATLAB定义的函数可能会有问题

% 计算DOA栅格数量
scan_len = size(scan_a,2);
%% CVX工具箱调用语法~开始
cvx_begin
variable s_cvx(scan_len,1)
minimize(  sum(pow_abs(s_cvx, p_norm)) );
subject to
    sum(pow_abs(u - scan_a * s_cvx,2)) <= tor_lim;
cvx_end
% CVX工具箱调用语法~结束
s_cvx = abs(s_cvx);
s_cvx = s_cvx ./ max(s_cvx);
end

space_smooth.m

function R_sig = space_smooth(rec_len, equ_l, x_sig)
% 基于空间平滑法解决相干信源问题
%   rec_len 真实阵元数量
%   equ_l   等效阵元数量
%   x_sig   回波矩阵，真实阵元数量*快拍数

% 计算平滑阵元数量
rec_len_equ = rec_len - equ_l + 1;
% 计算子阵的自相关
R_sig = zeros(equ_l,equ_l);
for i = 1:rec_len_equ
    % 第i个子阵
    foo = x_sig(i:equ_l+i-1,:);
    % 计算协方差
    R_sig = R_sig + (foo * foo')./rec_len;
end
R_sig = R_sig ./ rec_len_equ;
end