波束形成基于梯度下降上升法的旁瓣约束稳健自适应波束形成

量大龙神魔

已于 2023-12-01 09:56:33 修改

阅读量1.9k

点赞数 7

文章标签：算法信号处理信息与通信矩阵 matlab

于 2023-06-10 15:37:11 首次发布

本文链接：https://blog.csdn.net/qq_45471796/article/details/131074929

版权

前言

自适应波束形成器中，波束形成器能够在干扰方向自动形成凹槽来抑制干扰。对于某些突发干扰的情况，如干扰持续时间较短，自适应波束形成器来不及对干扰进行抑制。如果波束形成器的旁瓣较高，干扰就会对波束输出产生较大的影响，降低输出 SINR。因此，很多应用中需要对波束旁瓣进行抑制，以获得较低的旁瓣。

早期传统均匀线阵的旁瓣控制是通过窗函数来实现的，其灵感主要来源于均匀线列阵波束形成与时间序列分析的诸多相似之处。后来针对于波束所能获得的最低旁瓣的局限性，一些学者提出了零点展宽的方法，即只需要使波束形成器在干扰方向形成较深凹槽，足以抑制干扰即可，而不是要求整个旁瓣级降低。新世纪以来，凸优化发展迅猛，基于凸规划算法的波束形成器能够同时兼顾束旁瓣级、主瓣宽度、阵增益与稳健性等多个性能指标，在它们之间进行合理折中，以获得最佳的综合性能。

这篇文章便是一篇基于凸规划的旁瓣控制稳健自适应波束形成算法，是本人第一篇接触的旁瓣约束类波束形成算法，因此步骤可能有些错误及繁琐。

问题提出

因为旁瓣约束下的稳健自适应波束形成其核心思想主要是在抑制旁瓣的基础上保持期望信号方向的一个无失真响应，所以该优化问题可以表示为一个min-max问题，如下所示：

$\begin{array}{l}\min_{\boldsymbol{w}}\max_{\theta_m\in\Theta}\boldsymbol{w}^H\boldsymbol{a}(\theta_m)\boldsymbol{a}^H(\theta_m)\boldsymbol{w},s.t.\boldsymbol{w}^H\boldsymbol{\hat{R}_{i+n}}\boldsymbol{w}\leq\eta\eta_c, \boldsymbol{a}^H(\theta_0)\boldsymbol{w}=1\end{array}$

其中 $\Theta$ 为旁瓣所在的角度域，优化项的目的是使得旁瓣的输出功率的最大值最小，从而整体降低波束形成器的波束图旁瓣，而约束项的目的是使得期望信号方向输出无失真的同时抑制干扰信号。该问题很难通过CVX进行求解。

式中， $\boldsymbol{\hat{R}_{i+n}}$ 由期望信号角度邻域的补集进行线性积分而得到，即12年的线积分重构协方差阵的TSP，有：

$\pmb{\hat{R}}_{i+n}=\int_{\bar{\Theta}_0}\frac{\pmb{a}(\theta)\pmb{a}^H(\theta)}{\pmb{a}^H(\theta)\hat{R}_x^{-1}\pmb{a}(\theta)}d\theta$

对于 $\eta_c$ ，可由以下优化问题得到：

$\eta_c=\min_{\boldsymbol{w}}\boldsymbol{w}^H\hat{R}_{i+n}\boldsymbol{w},s.t.\boldsymbol{a}^H(\theta_0)\boldsymbol{w}=1\\~~~\Rightarrow \eta_c=\boldsymbol{w}^H\hat{R}_{i+n}\boldsymbol{w}\big|_{\boldsymbol{w}=\frac{\hat{R}_{i+n}^{-1}\textbf{a}(\theta_0)}{\textbf{a}^H(\theta_0)\hat{R}_{i+n}^{-1}\textbf{a}(\theta_0)}}=\frac{1}{\textbf{a}^H(\theta_0)\hat{R}^{-1}_{i+n}\textbf{a}(\theta_0)}$

对于 $\eta$ ，由于不知道其解集中是否含有最优权向量，一般情况下其值是大于1的。

作者随后将该min-max优化问题进行进一步转化，由重建协方差阵的性质可知，该阵为一正定阵，因此可进行Cholesky分解，有： $L^H L=\hat{R}_{i+n}$ 成立。同时，令 $\boldsymbol w_L=L\boldsymbol w$ ， $\boldsymbol{a}_L(\theta)=(L^H)^{-1}\boldsymbol{a}(\theta)$ ，则上式中的min-max问题可以转化为：

$\underset{\boldsymbol{w}_L}\min\underset{\boldsymbol{y}\in\mathcal{Y}}{\text{max}}\sum\limits_{m=1}^M y_m\boldsymbol{w}_L^H\boldsymbol{a}_L(\theta_m)\boldsymbol{a}_L^H(\theta_m)\boldsymbol{w}_L\\s.t.\boldsymbol{w}_L^H\boldsymbol{w}_L\leq\eta\eta_c,\boldsymbol{a}_L^H(\theta_0)\boldsymbol{w}_L=1,\mathcal{Y}=\{\boldsymbol{y}\mid\boldsymbol{y}\succeq\textbf{0},\sum_{m=1}^M y_m=1\}$

这实际上是对上一个min-max问题进行了进一步松弛，从旁瓣角度域的函数最大值转化为旁瓣角度域中每个元素的非负加权和，同时这个min-max问题是一个只关于 $\boldsymbol{w}_L$ 与 $\boldsymbol{y}$ 的问题，令 $\boldsymbol{x} = \boldsymbol{w}_L$ ，则该复数min-max问题可表示下列实数min-max问题为：

$\underset{\boldsymbol{x}\in\mathcal{X}}{\text{min}}\underset{\boldsymbol{y}\in\mathcal{Y}}{\text{max}}f(\boldsymbol{x},\boldsymbol{y}),\mathcal{X}=\{\boldsymbol{x}\in\mathbb{R}^{2N}|\boldsymbol{x}^T\boldsymbol{x}\leq\eta\eta_c,A_0^T\boldsymbol{x}=\boldsymbol{e}_1,\boldsymbol{e}_1=(1,0)^T\}\\$

$A_m= \left [ \begin{matrix} \mathcal{R}(\boldsymbol{a}_L(\theta_m))&-\mathcal{Z}(\boldsymbol{a}_L(\theta_m)) \\ \mathcal{Z}(\boldsymbol{a}_L(\theta_m))&\mathcal{R}(\boldsymbol{a}_L(\theta_m)) \end{matrix} \right ],\mathcal{Y}=\{\boldsymbol{y}\mid\boldsymbol{y}\succeq\textbf{0},\sum_{m=1}^M y_m=1\}$

其中 $\boldsymbol{x}$ 的前 $N$ 个元素为 $\boldsymbol{w}_L$ 的实部，后 $N$ 个元素为 $\boldsymbol{w}_L$ 的虚部，且若该问题的最优解为 $(\boldsymbol{x}^*,\boldsymbol{y}^*)$ ，则第一个min-max的优化问题的最优解为 $\boldsymbol{w}^*=L^{-1}(\boldsymbol{x}_{1:N}^*+j\boldsymbol{x}_{N+1:2N}^*)$ 。

梯度下降上升算法(Gradient Descent Ascent)及其应用

首先介绍下梯度下降上升算法。针对非凸且凹二元优化问题的梯度下降上升算法于2020年由Tianyi Lin等人在PMLR中提出，该问题及其算法定义如下：

有一个优化问题：

$\min\limits_{\text{x}\in\mathbb{R}^m}\max\limits_{\textbf{y}\in\mathcal{Y}}f(\textbf{x},\textbf{y})$

其中 $f$ 在 $\textbf{x}$ 的定义域内非凸但是在 $\textbf{y}$ 的定义域内为凹，则对于 $\mathbf{x}$ 的求解值 $\mathbf{\hat{x}}$ ，有以下迭代方法求解：

输入： $(\mathbf{x}_0,\mathbf{y}_0)$ ，步长 $(\eta_x,\eta_y)$
迭代： $\textbf{x}_t\leftarrow\textbf{x}_{t-1}-\eta_\textbf{x}\nabla_\textbf{x}f(\textbf{x}_{t-1},\textbf{y}_{t-1})$

$\textbf{y}_t\leftarrow\mathcal{P}_\mathcal{Y}\left(\textbf{y}_{t-1}+\eta_\textbf{y}\nabla_\textbf{y}f(\textbf{x}_{t-1},\textbf{y}_{t-1})\right)$

直至代价函数满足，得到 $\mathbf{\hat{x}}$

其中 $\textbf{x}$ 的迭代步骤为梯度下降法， $\textbf{y}$ 的迭代步骤为投影梯度上升法。介绍下投影梯度法：

针对一个封闭且凸的集合 $\mathcal{C}$ ，及在该集合内为凸的函数 $f$ ，其投影梯度法更新形式如下所示：

$\boldsymbol{x}^{t+1}=\mathcal{P}_{\mathcal{C}}\left(\boldsymbol{x}^t-\eta_t\nabla f\left(\boldsymbol{x}^t\right)\right)$

其中 $\mathcal{P}_{\mathcal{C}}(\cdot)$ 为将括号内的向量在凸集 $\mathcal{C}$ 上进行投影的操作。当 $f$ 在 $\textbf{x}$ 的定义域内为凸时 $\eta_c$ ，梯度下降法亦可换作投影梯度下降法，即：

$\textbf{x}_t\leftarrow\mathcal{P}_\mathcal{X}\left(\textbf{x}_{t-1}+\eta_\textbf{x}\nabla_\textbf{x}f(\textbf{x}_{t-1},\textbf{y}_{t-1})\right)$

之后介绍利普西茨常数(Lipcshitz constant)。对于 $f$ ，如果对其定义域 $\Omega$ 中的任意对元素 $(\mathbf{x},\mathbf{y})$ ，有以下式子满足：

$\|f(\mathbf{x})-f(\mathbf{y})\|\le L\|\mathbf{x}-\mathbf{y}\|$

则称满足该式子的 $L$ 称为利普西茨常数。

在上述定义的条件下，上一个实数min-max问题的梯度下降上升法迭代形式可表示为：

输入：迭代数 $T$ ，初始加权向量 $\boldsymbol{y}^1$ ， $\boldsymbol{y}$ 的步长 $\alpha_y$ ，重建后的干扰加噪声协方差阵 $\boldsymbol{\hat{R}_{i+n}}$ ，期望信号先验方向 $\theta_0$ ，离散旁瓣角度域 $\Theta$ ，松弛系数 $\eta$ ，约束参数 $\gamma$
计算：重建干扰加噪声协方差阵的Cholesky分解阵 $L$ ，松弛系数 $\eta_c$ ， $\boldsymbol{x}$ 的步长 $\alpha_{x}=\frac{1}{L_{x}}$ ，初始等价阵列权向量 $\boldsymbol{x}^1 = (\mathcal{R}(\boldsymbol{w}_c)^T,\mathcal{Z}(\boldsymbol{w}_c)^T)^T$
迭代： $\boldsymbol{x}^{t+1}\leftarrow\mathcal{P}_{\mathcal{X}}(\boldsymbol{x}^t-\alpha_{\chi}\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^t,\boldsymbol{y}^t))$

$\textbf{y}^{t+1}\leftarrow\mathcal{P}_\mathcal{Y}(\textbf{y}^t+\alpha_y\frac{\boldsymbol{\eta}^t}{\|\boldsymbol{\eta}^t\|_1})$

$p^t\gets\text{max}_m||A_m^T\pmb{x}^{t+1}||^2$

输出： $\boldsymbol{\hat{w}}=L^{-1}(\textbf{x}_{1:N}^{T_{m}+1}+j\textbf{x}_{N+1:2N}^{T_{m}+1})$ ，其中 $T_m=\text{arg}\text{min}_tp^t$

针对该流程，首先考虑迭代步长，关于迭代步长的选取及其与梯度下降法、投影下降法建议参考以下三个专栏：

梯度下降的收敛性分析 (the Convergence Analysis of Gradient Descent)

【非线性优化】梯度的Lipschitz连续性

优化算法（四）——投影梯度法解约束问题

本人在此仅针对当前问题所需要涉及到的进行简单说明。

考虑迭代步长 $\alpha_x$ ，针对该优化问题的优化项函数 $f(\boldsymbol{x},\boldsymbol{y})$ ，根据梯度下降引理，当 $0 \leq \alpha_x \leq 2/L$ 时，目标函数值的大小始终保持单调递减状态，当优化步长选取 $1/L$ 时，目标函数值的缩减速度最快。而目标函数值的二阶偏导与 $L$ 存在着以下关系： $\|\nabla^2_{\boldsymbol{x}}f(\boldsymbol{x},\boldsymbol{y})\|\leq L$ ，因此可通过二阶偏导大致确定迭代步长，具体如下：

$\|\nabla_{\boldsymbol{x}}^2f(\boldsymbol{x},\boldsymbol{y}))\|=2\left \| \sum\limits_{m=1}^M y_mA_mA_m^T \right \|=2\lambda_{max}\left(\sum\limits_{m=1}^M y_mA_mA_m^T\right)\\=2\underset{\|\boldsymbol{v}\|=1}{\text{max}}\sum\limits_{m=1}^M y_m\boldsymbol{v}^T\boldsymbol{a}_L(\theta_m)\boldsymbol{a}_L^H(\theta_m)\boldsymbol{v}\leq2\text{max}_m\|\pmb a_L(\theta_m)\|^2\\\Rightarrow L_x=2\text{max}_m\|\pmb a_L(\theta_m)\|^2\Rightarrow \alpha_X=\frac{1}{L_X}$

考虑迭代步长 $\alpha_y$ ，因其二阶偏导数为0，难以快速确定其迭代步长，因此作者此处设定其为0.005。

考虑 $\boldsymbol{x}$ 的迭代操作 $\boldsymbol{x}^{t+1}\leftarrow\mathcal{P}_{\mathcal{X}}(\boldsymbol{x}^t-\alpha_{\chi}\nabla_{\boldsymbol{x}}f(\boldsymbol{x}^t,\boldsymbol{y}^t))$ ，其中 $\boldsymbol{x}$ 的投影操作可表示为：

$\mathcal{P}_{\mathcal{X}}(\pmb{x}_0)=\min_{\pmb{x}\in\mathcal{X}}\|\pmb{x}-\pmb{x}_0\|$

即对于集合 $\mathcal{X}$ ，该集合外的一向量 $\boldsymbol{x}_0$ 在该集合上的投影即为其欧式距离最小的 $\boldsymbol{x}\in\mathcal{X}$ 。在考虑约束后，该问题为一凸问题，可用CVX进行求解，作者此处根据KKT条件对其进行直接求解。简单介绍下KKT条件：

有一有约束优化问题，为：

$\min f_0(x)~~~~s.t.~~f_i(x) \leq 0,i=1,...,m~~~~h_i(x) = 0,i=1,...,p$

易知该问题的对偶函数为： $g(\lambda,v) =\inf _x \{f_0(x)+\sum_{i=1}^{m}\lambda_if_i(x)+\sum_{j=1}^{p}v_ih_i(x)\}$

该问题的拉格朗日函数为： $L(x,\lambda,v) ={f_0(x)+\sum_{i=1}^{m}\lambda_if_i(x)+\sum_{j=1}^{p}v_ih_i(x)}$

则其五个KKT条件如下所示：

原问题可行性1： $f_i(x) \leq 0,i=1,...,m$
原问题可行性2： $h_i(x) = 0,i=1,...,p$
对偶问题可行性： $\lambda_i \geq 0,i=1,...,m$
互补松弛条件： $\lambda_i^*f_i(x^*) = 0,i=1,...,m$
稳定性条件： $\frac{\partial L(x,\lambda^*,v^*)}{\partial x}\big|_{x=x^*} = 0$

其中 $()^*$ 为其最优解，且当该问题为一般问题时，解不一定满足KKT，且满足KKT的解不一定最优；当该问题为凸问题，且slater条件成立时，最优解一定满足KKT。

则考虑约束后的该问题可表示为：

$\min_{\boldsymbol{x}\in\mathcal{X}}\|\boldsymbol{x}-\boldsymbol{x}_0\|^2~~~~s.t.\boldsymbol{x}^T\boldsymbol{x}\leq \eta\eta_c~~~~A^T_0 \boldsymbol{x} = \boldsymbol{e}_1$

首先证Slater条件成立：

将 $\boldsymbol{x}_c = (\mathcal{R}(\boldsymbol{w}_c)^T,\mathcal{Z}(\boldsymbol{w}_c)^T)^T$ 带入约束函数，即得： $\boldsymbol{x}_c^T\boldsymbol{x}_c = \eta_c\leq \eta\eta_c$ ， $A^T_0 \boldsymbol{x}_c = \boldsymbol{e}_1$ ，即证Slater条件成立，满足KKT条件的解为最优解的充要条件，因此该优化问题的拉格朗日函数为：

$L(\boldsymbol{x},\mu,\boldsymbol{v}) =\|\boldsymbol{x}-\boldsymbol{x}_0\|^2+\mu(\boldsymbol{x}^T\boldsymbol{x}- \eta\eta_c)+\boldsymbol{v}^T(A^T_0 \boldsymbol{x} - \boldsymbol{e}_1)$

KKT条件的解为：

$\boldsymbol{x}^{*T}\boldsymbol{x}^*\leq \eta\eta_c$
$A^T_0 \boldsymbol{x}^* = \boldsymbol{e}_1$
$\mu^* \geq 0$
$\mu^*(\boldsymbol{x}^{*T}\boldsymbol{x}^*- \eta\eta_c)=0$
$2(\boldsymbol{x}^*-\boldsymbol{x}_0)+A_0\boldsymbol{v}^*+2\mu^*\boldsymbol{x}^* = 0$

考虑5.，得到 $\boldsymbol{x}^*$ 的初步解，有：

$\boldsymbol{x}^* = (2+2\mu^*)^{-1}(2\boldsymbol{x}_0-A_0\boldsymbol{v}^*)$

带入2.，得：

$\boldsymbol{v}^* = (A_0^TA_0)^{-1}(2A_0^T\boldsymbol{x}_0-(2+2\mu^*)\boldsymbol{e}_1)$

继续考虑5.，有：

$\boldsymbol{x}^* = (2+2\mu^*)^{-1}(2\boldsymbol{x}_0-A_0( (A_0^TA_0)^{-1}(2A_0^T\boldsymbol{x}_0-(2+2\mu^*)\boldsymbol{e}_1)))\\~~~~=(1+\mu^*)^{-1}\boldsymbol{x}_0-(1+\mu^*)^{-1}A_0 (A_0^TA_0)^{-1}A_0^T\boldsymbol{x}_0+A_0 (A_0^TA_0)^{-1}\boldsymbol{e}_1$

作者令 $P_{\mathcal{N}(A_0^T)}=I-A_0(A_0^T A_0)^{-1}A_0^T$ ， $\boldsymbol{b}_p=A_0(A_0^T A_0)^{-1}\boldsymbol{e}_1$ ，则 $\boldsymbol{x}^*$ 的进一步解为：

$\boldsymbol{{x}}^*=(1+{\mu}^*)^{-1}P_{\mathcal{N}(A_0^T)}\boldsymbol{x}_0+\boldsymbol{b}_p$

注意到 $\boldsymbol{b}_p^T\boldsymbol{b}_p=\boldsymbol{e}_1^T((A_0^T A_0)^{-1})^TA_0^TA_0(A_0^T A_0)^{-1}\boldsymbol{e}_1=\eta_c$ ， $\boldsymbol{b}_p^T P_{\mathcal{N}(A_0^T)}\boldsymbol{x}_0=\boldsymbol{e}_1^T((A_0^T A_0)^{-1})^TA_0^T(I-A_0(A_0^T A_0)^{-1}A_0^T)=0$ ，有：

${\boldsymbol{x}}^{*T}{\boldsymbol{x}^*}=(1+{\mu}^*)^{-2}\boldsymbol{x}_0^{T} P_{\mathcal{N}(A_0^T)}^T P_{\mathcal{N}(A_0^T)}\boldsymbol{x}_0+\eta_c=\eta\eta_c$

即得：

$\mu^* = \sqrt{\frac{\boldsymbol{x}_0^T P_{N(A_0^T)}\boldsymbol{x}_0}{(\eta-1)\eta_c}}-1$

带入 $\boldsymbol{x}^*$ 的进一步解，即得：

$\boldsymbol{{x}}^*=(\sqrt{\frac{\boldsymbol{x}_0^T P_{N(A_0^T)}\boldsymbol{x}_0}{(\eta-1)\eta_c}})^{-1}P_{\mathcal{N}(A_0^T)}\boldsymbol{x}_0+\boldsymbol{b}_p=~\sqrt{\frac{(\eta-1)\eta_c}{\boldsymbol{x}_0^T P_{N(A_0^T)}\boldsymbol{x}_0}}(I-~A_0(A_0^T A_0)^{-1}A_0^T)+~A_0(A_0^T A_0)^{-1}\boldsymbol{e}_1$

其中：

$\mu^* \geq 0 \Rightarrow \sqrt{\frac{(\eta-1)\eta_c}{\boldsymbol{x}_0^T P_{N(A_0^T)}\boldsymbol{x}_0}} \leq 1$

考虑 $\boldsymbol{y}$ 的迭代操作 $\boldsymbol{y}^{t+1}\leftarrow\mathcal{P}_\mathcal{Y}(\boldsymbol{y}^t+\alpha_y\frac{\boldsymbol{\eta}^t}{\|\boldsymbol{\eta}^t\|_1})$ ，因 $f(\cdot,\textbf{y})$ 为一凹函数，则其极值点为最大处，因此梯度方向为其梯度上升方向。首先对min-max问题的优化函数求偏导，有：

$\nabla_{\boldsymbol{y}}f(\boldsymbol{x}^t,\boldsymbol{y})=(\|A_1^T\pmb{x}^t\|^2,\|A_2^T\pmb{x}^t\|^2,\dots,\|A_M^T\pmb{x}^t\|^2)^T$

该偏导数实际上反应的是旁瓣区域各离散角度的空域功率，为了进一步抑制旁瓣区域功率较大的角度，作者将该梯度向量改写为：

$\eta_m^t=\begin{cases}\|A_m^T\boldsymbol{x}^t\|^2&\|A_m^T\boldsymbol{x}^t\|^2\geq\gamma\max\limits_k\|A_k^T\boldsymbol{x}^t\|^2\\ 0&{otherwise}\end{cases}$

即设置一门限函数，使得该优化问题得到的权向量对超出门限比例 $\gamma$ 的部分进行优先抑制。同时，投影操作 $\mathcal{P}_\mathcal{Y}(\cdot)$ 可表示为：

$\mathcal{P}_\mathcal{Y}(\boldsymbol{y}_0){=}\frac{\boldsymbol{y}_0}{\sum_{m=1}^M y_{0m}}$

考虑 $p^t\Leftarrow\text{max}_m||A_m^T\pmb{x}^{t+1}||^2$ ，该函数实际上对每次迭代的旁瓣区域的最大功率进行记录。

因此 $t$ 次迭代过程中最小的 $p^t$ 后的迭代对应的权向量即为该算法所得的最优权向量。

此外，作者还考虑了附加旁瓣约束的前提，即人为构造波束图，设其旁瓣附加约束函数为 $DSL(\theta)$ ，则原二元优化问题的优化项可表示为：

$f(\boldsymbol{x},\boldsymbol{y})=\sum_{m=1}^{M}y_m10^{-\frac{DSL(\theta_m)}{10}}\boldsymbol{x}^TA_mA_m\boldsymbol{x}$

仿真

部分仿真条件如表所示：

阵元数	32
信噪比	10
快拍数	256
扰噪比	30
阵元间距	半波长
信源数	1
干扰数	2
信源到达角	5°
干扰到达角	-30°，45°
信号域宽度	10°
角度采样步长	0.2°
附加旁瓣约束	-40dB，[-35°,-25°]，[40°,50°]
迭代次数	100
旁瓣功率谱峰门限 $\gamma$	0.95
权重向量 $\boldsymbol{y}$ 步长	0.005
松弛系数 $\eta$	1.1

代码如下所示：

%% 初始化及设定参数
clc;
clear;
jk = 0;
snr = 10;
array_num = 32;%%阵元数
deg_dev_predf = 4;
sample_acc = 0.2;
ang_mismatch1 = 0;
ang_mismatch2 = 0;
ang_mismatch3 = 0; %%三个不同方向的DOA误差
snapshot_num = 256;%%快拍数
source_aoa = [5,-30,45];%%信源到达角
tgt_ang = 1;
source_dev = source_aoa + [ang_mismatch1,ang_mismatch2,ang_mismatch3];%%DOA得到的带有误差的
c = 340;%%波速
f = 1000;%%频率
lambda = c/f;%%波长
d = 0.5*lambda;
inr1 = 30;
inr2 = 30;
source_num = length(source_aoa);%%信源数
sig_nr = [snr,inr1,inr2];%%信噪比、扰噪比
sig_range = 5;
desire_side_lobe_lvl = -40;
other_side_lobe_lvl = -30;
desire_side_lobe_wide = 5;
cons_gamma = 0.95;
alpha_y = 0.005;
eta = 1.1;
iter_times = 100;
%% 导向矢量
A_bar = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(source_dev));%%阵列响应矩阵--每个导向矢量都带误差
A = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(source_aoa));%%阵列响应矩阵
tar_sig = wgn(1,snapshot_num,sig_nr(tgt_ang));
inf1 = wgn(1,snapshot_num,sig_nr(2));
inf2 = wgn(1,snapshot_num,sig_nr(3));
n = (randn(array_num,snapshot_num)+randn(array_num,snapshot_num)*1i)/sqrt(2);
rec_sig = A(:,1)*tar_sig+A(:,2)*inf1+A(:,3)*inf2+n;%%接收信号
R = rec_sig*rec_sig'/snapshot_num;%%协方差矩阵
Rs = A(:,1)*tar_sig*(A(:,1)*tar_sig)'/snapshot_num;
Ri = (A(:,2)*inf1+A(:,3)*inf2)*(A(:,2)*inf1+A(:,3)*inf2)'/snapshot_num;
Rn = n*n'/snapshot_num;
%% 算法
sample_set = [-90:sample_acc:source_aoa(1)-sig_range,source_aoa(1)+sig_range:sample_acc:90];
%%DSL计算
DSL = zeros(1,length(sample_set));
for ik = 1:length(sample_set)
    for jk = 2:source_num
        if sample_set(ik) <= source_aoa(jk)+desire_side_lobe_wide && sample_set(ik) >= source_aoa(jk)-desire_side_lobe_wide
            DSL(ik) = 10^(-1*desire_side_lobe_lvl/10);
        end
    end
end
for ik = 1:length(sample_set)
    if DSL(ik) == 0
        DSL(ik) = 10^(-1*other_side_lobe_lvl/10);
    end
end
%%INCM重建
R_recons = zeros(array_num,array_num);
for ik = 1:length(sample_set)
    a_sample_set = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(sample_set(ik)));
    R_recons = R_recons+a_sample_set*a_sample_set'*abs(1/(a_sample_set'*inv(R)*a_sample_set));
end
a_0 = exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(source_aoa(1)));
w_c = inv(R_recons)*a_0/(a_0'*inv(R_recons)*a_0);
eta_c = 1/(a_0'*inv(R_recons)*a_0);
%%Cholsky分解
[L_chol] = chol(R_recons);
A_L_theta = inv(L_chol')*exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(sample_set));
A_L_theta_norm = zeros(1,length(sample_set));
for ik = 1:length(sample_set)
    A_L_theta_norm(ik) = norm(A_L_theta(:,ik),2);
end
%%x步长计算
L_x = 2*max(A_L_theta_norm);
alpha_x = 1/L_x;
x_iter = zeros(2*array_num,iter_times+1);
y_iter = zeros(length(sample_set),iter_times+1);
y_iter(1,1) = 1;%%初始值设定为(1,0,0,...,0)^T
x_iter(:,1) = [(real(w_c))',(imag(w_c))']';
p_t = zeros(1,iter_times);%%旁瓣区域谱峰记录
A_L_theta_0 = inv(L_chol')*exp(-1i*(0:array_num-1)'*2*pi*(d/lambda)*sind(source_aoa(1)));
A_0_theta = [real(A_L_theta_0),-1*imag(A_L_theta_0);imag(A_L_theta_0),real(A_L_theta_0)];%%A_0
e_1 = [1;0];
for iter = 1:iter_times
    disp(iter);
    %%求x偏导
    diff_x_f_x_y = zeros(2*array_num,1);
    for ik = 1:length(sample_set)
        A_m_theta = [real(A_L_theta(:,ik)),-1*imag(A_L_theta(:,ik));imag(A_L_theta(:,ik)),real(A_L_theta(:,ik))];
        diff_x_f_x_y = diff_x_f_x_y+2*DSL(ik)*y_iter(ik,iter)*(A_m_theta*A_m_theta')*x_iter(:,iter);
    end
    x_iter_mid = x_iter(:,iter)-alpha_x*diff_x_f_x_y;
    %%x投影
    cvx_begin quiet
    variables proj_x_iter(array_num*2,1) 
    minimize((proj_x_iter-x_iter_mid)'*(proj_x_iter-x_iter_mid));
    subject to 
        proj_x_iter'*proj_x_iter <= abs(eta*eta_c);
        A_0_theta'*proj_x_iter == e_1;
    cvx_end
    x_iter(:,iter+1) = proj_x_iter;
    %%求y偏导
    diff_y_f_x_y = zeros(1,length(sample_set));
    for ik = 1:length(sample_set)
        A_m_theta = [real(A_L_theta(:,ik)),-1*imag(A_L_theta(:,ik));imag(A_L_theta(:,ik)),real(A_L_theta(:,ik))];
        diff_y_f_x_y(ik) = DSL(ik)*norm(A_m_theta'*x_iter(:,iter),2)^2;
    end
    eta_t = zeros(length(sample_set),1);
    for ik = 1:length(sample_set)
        if diff_y_f_x_y(ik) >= cons_gamma*max(diff_y_f_x_y)
            eta_t(ik) = diff_y_f_x_y(ik);
        else
            eta_t(ik) = 0;
        end
    end
    y_iter_mid = y_iter(:,iter)+alpha_y*eta_t/sum(abs(eta_t));
    %%y投影
    y_iter(:,iter+1) = y_iter_mid/sum(abs(y_iter_mid));
    %%求p_t
    p_t_mid = zeros(1,length(sample_set));
    for ik = 1:length(sample_set)
        A_m_theta = [real(A_L_theta(:,ik)),-1*imag(A_L_theta(:,ik));imag(A_L_theta(:,ik)),real(A_L_theta(:,ik))];
        p_t_mid(:,ik) = (A_m_theta'*x_iter(:,iter+1))'*(A_m_theta'*x_iter(:,iter+1));
    end
    p_t(iter) = max(p_t_mid);
end
[~,index_m] = min(p_t);
x_target = x_iter(:,index_m(1)+1);
w0 = inv(L_chol)*(x_target(1:array_num,:)+1i*x_target(array_num+1:end,:));
beam_plot(w0);

function [] = beam_plot(input_weight_vector)
    array_num = length(input_weight_vector);
    theta = -90:0.01:90;
    p = exp(-1j*2*pi*(0:array_num-1)'*sind(theta)/2);
    y = input_weight_vector'*p;
    outputval = 20*log10(abs(y)/max(abs(y)));
    figure()
    plot(theta,outputval,'LineWidth',2);
end

仿真结果如下所示：