【ISAC】paper_NOMA Empowered Integrated Sensing and Communication

NOMA: Enhanced Integrated Sensing and Communication System, Non-Public Conference Proceedings

[1]. 关于MIMO雷达信号探测的设计研究
[2]. 对信号的频谱分析研究
[3]. D. R. Fuhrmann和G. San Antonio合著，《使用信号交叉相关实现MIMO雷达系统传输波束形成的》发表于IEEE transactions on aerospace and electronic systems期刊第44卷第1期（pp 171-186），时间是2008年一月
[4] 线性代数及其在工程中的应用研究

文章目录

• 1-Model

• • Communication Model
  • Sensing Model
  • Problem Formulation

• 2-Solution

• • 计算信噪比

第三代编码技术（3-code）研究与优化设计

基于信道状态信息的信道质量评估模型

基于信道状态信息的信道质量评估模型

非正交多址（NOMA）速率优化算法研究

非正交多址（NOMA）速率优化算法研究

基于信道状态信息的信道质量评估模型

非正交多址（NOMA）速率优化算法研究

• 4-附录
• • 4.1- The quadratic form of the covariance matrix is non-convex.
• • 4.2- Estimating the cross-correlation between the target signal.

1-Model

双功能基站（BS）配备了一组由N个天线组成的均匀线性阵列（ULA）。

Communication Model

• \mathcal{K}=\{1,\ldots,K\}: K single-antenna users indexed.

• \mathcal{M}=\{1,\ldots,M\}: M radar targets indexed.

• \mathbf{w}_is_i\text{ for }\forall i\in\mathcal{K}: where \mathbf{w}_i\in\mathbb{C}^{N\times1} are beamformers for delivering the information symbol s_i to user i. BS transmits the superimposed signals to all users in downlink. \mathbf{w}_i是波束形成器，s_i是发给第i个用户的信息。

• 第k个用户的接收信号，y_k=\mathrm{h}_k^H\sum_{i\in\mathcal{K}}\mathrm{w}_is_i+n_k=\sum_{i\in\mathcal{K}}\mathrm{h}_k^H\mathrm{w}_is_i+n_k

• 信道增益系数，\mathrm{h}_k=\Lambda_k^{-1/2}\widetilde{\mathsf{h}}_k,\forall k\in\mathcal{K} denotes the BS-user channel, \Lambda_k^{-1/2} and \widetilde{\mathbf{h}}_k\in\mathbb{C}^{N\times1} denote the large and small scale fading, respectively, and n_k denotes the circularly symmetric complex Gaussian noise with variance \sigma_n^2.

• 不同用户的大尺度衰弱系数不同，即，\Lambda _1^{- 1}\leq \Lambda _2^{- 1}\leq \cdots \leq \Lambda _K^{- 1}. We assume that the users’ indexes are in an increasing order with respect to their large-scale channel strength. 用户1最弱.

• R_{k\to k}=\log_2\left(1+\frac{|\mathrm{h}_k^H\mathrm{w}_k|^2}{\sum_{i\in\mathcal{K},i>k}|\mathrm{h}_k^H\mathrm{w}_i|^2+\sigma_n^2}\right), the achievable rate of s_k after SIC at user k for \forall k\in\mathcal{K},k\neq K. 第K个用户的可达码率.

• Note: 要解码第k个用户，先利用将前k-1移除，然后将k之后的用户当做噪声. NOMA解码是逐个的过程，e.g., 解第k个，那么k后面的都是噪声，可以得到R_{k\to k}，再解第k+1个，利用SIC，为什么可以利用SIC，因为前面k个已经知道了可以直接减去，这也是要做"不同用户的大尺度衰弱系数不同"假设的原因. 有个特殊的用户，就是最后一个解码的用户，那么前面都已经解出来了，噪声只剩下AWGN，没有干扰.
  在 SIC 过程中，先解码强用户的信号，再依次解码弱用户的信号，并将已解码的信号视为已知干扰，从接收到的信号中减去。

• R_k = \min \left \{ R_{k \to k}, \ldots, R_{k \to K} \right \}. 整体可达速率对于每个k \in \mathcal{K}且k \neq K来说是s_k. 为何取最小值？若任一用户无法正确解码(s_k), 则该信号将对后续用户的正常工作造成干扰并破坏整个系统运行机制. 换句话说, 系统的整体传输性能受限于最差情况下的传输效果.

• 因此, 系统整体速率取决于该信号在最不利情况下（即最小速率）的信道状态.

• R_K = \log_2 \left( 1 + \frac{|{\bf h}_K^H {\bf w}_K|^2}{\sigma_n^2} } \right) 当处理最后一个接入者时, 前面所有接入者均已成功完成数据传输. 因此可达码率不受其他接入者的干扰.

• 其吞吐量计算公式为R = \sum_{k \in {\cal K}} R_k

Sensing Model

• \mathrm{R}_\mathbf{w}=\sum_{i\in\mathcal{K}}\mathrm{w}_i\mathrm{w}_i^H：通过优化发射信号来实现感知与生成过程之间的等同性。
  • 其中涉及的技术包括解释协方差以及其相关的性质。
    • 详细解释了协方差的概念
    • 详细阐述了协方差矩阵的概念及其作用
    • 共轭转置是一种重要的线性代数操作，在信息处理中被广泛应用
  • P(\theta_m)=\mathrm{a}^H(\theta_m)\mathrm{R}_\mathbf{w}\mathrm{a}(\theta_m)：旨在最大化目标探测信号的能量。
    • 探讨这一形式的原因及背后的数学原理
    • 见[1]中的详细推导过程

\begin{aligned}\sum_{m=1}^Me^{-j2\pi f_0\tau_m(\theta)}x_m(n)&\triangleq\mathbf{a}^*(\theta)\mathbf{x}(n),\quad n=1,\ldots,N\end{aligned}
在M根天线的MIMO雷达系统中，基带信号\mathbf{x}(n)打到被感知的目标（那么信号就带有目标参数\theta），这个信息包含在e^{-j2\pi f_0\tau_m(\theta)}中，f_0是发射雷达信号的载频，\tau_m(\theta)是信号从发射到目标上所需要的时间，注意这个信号不是回波. 左边是逐元素点积，等价右边向量乘积。上述定义首次定义在[2, pp270]。
此时信号的功率是，P(\theta)=\mathbf{a}^*(\theta)\mathbf{Ra}(\theta), \mathbf{R}是\mathbf{x}(n)的协方差矩阵，i.e.,\mathbf{R}=E\{\mathbf{x}(n)\mathbf{x}^*(n)\}.，为什么感知信号的功率是这样二次型的形式，见[3]，什么是二次型，见线性代数及其应用.

• \theta_m,\forall m\in\mathcal{M} are target directions
• a(\theta_m)= [1,e^{j\frac{2\pi}\lambda d\sin(\theta_m)},\ldots,e^{j\frac{2\pi}\lambda d(N-1)\sin(\theta_m)}]^T is the steering （转向）vector（用于描述天线阵列对来自特定方向）, where \lambda and d denote the carrier wavelength and antenna spacing（波长和天线间距）, respectively. 有N更天线同时发射信号感知，所以有N项.

其中：

• \lambda 是信号的波长。

• d 是天线阵列中相邻天线之间的距离。

• N 是天线的数量。

• \theta_m 是信号入射的角度。

解释：

• 1 ：第一个元素 1 对应于第一个天线元的相位响应，假设参考相位为零。

• e^{j\frac{2\pi}{\lambda}d\sin(\theta_m)} ：第二个元素表示第二个天线元接收到的信号相对于第一个天线元的相位偏移。相位偏移取决于信号的波长 \lambda 、天线元之间的距离 d 以及信号入射角 \theta_m 。

• e^{j\frac{2\pi}{\lambda}d(N-1)\sin(\theta_m)} ：最后一个元素表示第 N 个天线元的相位响应。
  转向矢量的作用：
  转向矢量在阵列处理中用于构建波束形成器、执行信号方向估计和其他空间滤波任务。通过调整转向矢量的相位，可以指向或聚焦到特定的入射角 \theta_m 上，从而增强特定方向的信号，或抑制来自其他方向的干扰。

此式表明，在给定方向θ时所接收信号的能量度量即为该式的结果。
该能量值是由转向矢量\bm a(\theta)与协方差矩阵\bm R进行点积运算所得的结果。
从直观上讲，
其中相关函数矩阵\bm R通常被定义为信号样本经过均值归一化后的自相关矩阵，
其物理意义即为该方向上的信号强度。
将不同方向上的能量值进行加权求和即可获得整体感知能力。

信号的相关矩阵与信号功率之间的关系
相关矩阵的迹与信号功率 ：信号的相关矩阵 \mathbf{R}_s 包含了信号的全部统计信息，包括其分量间的相互作用。而信号的总功率可以通过相关矩阵的迹来计算。因此，相关矩阵的迹等于信号的平均功率。
相关矩阵的元素 ：相关矩阵的对角线元素 R_{s,ii} 表示信号第 i 个分量的自相关，即这个分量的平均功率。而非对角线元素 R_{s,ij} 则表示信号第 i 和第 j 个分量之间的互相关（交叉功率）。
信号功率的分布 ：相关矩阵的特征值可以用来 描述信号功率在不同方向上的分布 。特征值越大，意味着信号在对应方向上的功率越强。

• C( \theta _{k}, \theta _{p}) = | \mathrm{a} ^{H}( \theta _{k}) \mathbf{R} _{\mathbf{w} }a(\theta_{p})|,\forall k\neq p\in\mathcal{M}.，两个感知目标参数之间的交叉相关， cross-correlation，这是一个强假设，两个感知目标信号的互相关很小。
• \bar{C}=\frac2{M^2-M}\sum_{k=1}^{M-1}\sum_{n=k+1}^MC(\theta_k,\theta_p)^2，均方交叉相关，mean-squared cross-correlation

Problem Formulation

目标是在通信中实现加权求和的最大化 \max _{\mathbf{w}} \rho_c \rho_{k \in \mathcal{K}} R_k+\rho_r \sum_{m \in \mathcal{M}} P\left(\theta_m\right)
受限于：每个子信道k满足最低性能R_{\min , k}, \forall k \in \mathcal{K}
同时保证：相邻子信道间的性能波动不超过P_{\text {diff }}, \forall k \neq p \in \mathcal{M}
其中各子信道权重向量模长平方后的对角矩阵等于均匀分配功率谱\frac{P_t \mathbf{1}^{N \times 1}}{N}
此外还满足平均载波效率\bar{C} \leqslant\xi

如何理解这个优化问题？优化就是在资源有限的情况下进行最有效配置,不然我当然希望通信速率和感知性能最大化。
在此系统中,考虑下行 NOMA 系统,基站能够发射 ISAC 双功能波形,服务于 K 个用户.很自然地,在这个系统中我的优化目标就应为所有用户的总传输速率与感知能力之和（否则你说在通信方面你该优化啥）,
因为这个 beamformers 还具备感知功能,我们当然希望在存在目标方向上的感知功率达到最大值,两者加权构成目标函数.第一个约束很自然地就必须满足每个用户都能接收到信号吧,第二个约束:不在目标方向上的感知功率需实现公平分配,第三个约束:每个天线的功率约束,第四个约束:实现自适应感知.

对角线元素 \left[\mathbf{P}_{\text{total}}\right]_{kk} 是总功率矩阵的第 k 个对角线元素.这个元素表示第 k 个天线接收到的总功率.具体来说:\left[\mathbf{P}_{\text{total}}\right]_{kk} = \left(\sum_{i \in \mathcal{K}} \mathbf{w}_i \mathbf{w}_i^H \right)_{kk}

该优化问题的主要挑战体现在两个方面：其一，在目标函数中所涉及的目标码率既不具有凸性也不具备凹性；其二，在第2、3条约束条件中均为非-convex形式。

连续凸近似体现在两方面：其一为目标函数逼近；其二则通过定义\mathbf{W}_k\triangleq\mathbf{w}_k\mathbf{w}_k^H使约束变得凸特性的同时引入了秩1约束非凸性。尽管可以通过半正定松弛方法解决这一问题但该方法存在性能代价以及可能失去解的有效性；因此在优化过程中选择将秩1约束纳入惩罚项与目标函数联合优化以平衡这两方面的问题

• \rho_c与\rho_r通过调节通感性能实现平衡，并被定义为正则系数。

  • \left|P(\theta_k)-P(\theta_p)\right|\leq P_{\text {diff }}, 对所有k≠p∈𝓜, 确保了不同感知目标之间的感知功率相近。
  • 矩阵\operatorname{diag}\left(\sum_{i∈𝓚}\bm w_i{\bm w_i}^\mathrm H\right)的对角线元素定义了每个天线的约束常数, 并将其平均功率分配给每根天线。

\mathbf{w}_i : 表示第 i 个用户的波束形成向量，其维度通常为 N× ，其中 N 表示天线的数量。

\mathbf{w}_i \mathbf{w}_i^H 代表的是波束形成向量的外积运算结果，并形成了一个 N \times N 的矩阵来表示该波束形成向量所对应的协方差矩阵。

\sum_{i \in \mathcal{K}} \mathbf{w}_i \mathbf{w}_i^H 是用于计算所有用户协方差矩阵的总和，并生成一个大小为 N \times N 的矩阵；该矩阵代表了在所有用户波束形成向量作用下天线阵列的整体空间信道状态信息。

\operatorname{diag}(\cdot)

\frac{P_t \mathbf{1}^{N \times 1}}{N} 表示一个常量向量，在其各个元素位置上都赋值为 \frac{P_t}{N} 的值，并用于表示各天线均分的平均功率分配。

\bar{C}被设定为均方互相关的一个理想上界，在上述假设下得以确定。
由于目标函数仅依赖于可达码率这一单一指标，在这种情况下其凸性特征变得不明显。
此外，在这种情况下该约束条件组成为非凸形式。
协方差矩阵的二次形式导致相邻方向之间的功率差异受限，
即\left|P\left(\theta_k\right)-P\left(\theta_p\right)\right| \leq P_{\text {diff }}对于所有k≠p∈M成立，
同时该矩阵对角线元素满足特定条件：
\operatorname{diag}\left(\sum_{i∈K} w_i w_i^H\right)=\frac{P_t 1^{N×1}}{N}，
从而使得整体结构呈现非凸特性。

2-Solution

称\bm W_k为由向量\bm w_k构成的外积矩阵；该矩阵满足半正定性、对称性以及秩为1等性质；为何会引入这样的定义？

为了确定\mathbf{w}_k\mathbf{w}_k^H是否为半正定矩阵（positive semidefinite），我们注意到该矩阵是symmetric的，并基于二次型理论将其代入半正定性的定义中进行验证。具体证明参考文献中的第4卷第7-2章习题25。

为了最大化函数f的最大值\max _{\gamma_k, W_k} f（\rho_c乘以\sum_{k∈𝒦} \gamma_k）加上\rho_r乘以\sum_{m∈𝒴} P(θ_m)），受限于以下条件：对于所有的k∈𝒦有R_k≥γ_k；对于所有的k∈𝒦有γ_k≥R_{min,k}；对于所有的k∈𝒦有\mathrm{W}_k⪰0并且\mathrm{W}_k=\mathrm{W}_k^H；以及对于所有的k∈𝒦有\mathrm rank(\mathrm W_k)=1

为什么|\mathrm{h}_j^H \mathbf{w}_k|^2 =\operatorname{Tr}(\mathrm{H}_j \mathbf{W}_k)

基于迹运算的基本属性\operatorname Tr({\bm A}{\bm B})=\operatorname Tr({\bm B}{\bm A})以及\operatorname Tr({\bm a}{\bm b}\llap{\rm H}})={\bm b}\llap{\rm H}}{\bm a}

计算结果表明：
经过推导得出：
\operatorname{Tr}(\mathrm{h}_{j}\mathrm{h}_{j}^{H}\mathbf{w}_{k}\mathbf{w}_{k}^{H})=\operatorname{linspace}(\mathbf{\omega}_{k}^{H}\mathrm{\theta}_{j}\mathrm{\theta}_{j}^{H}\mathbf{\omega}_{k})=\left(\\right)
利用求迹的线性性质进行交换，
将其分解为两个向量的内积运算，
最终得出其模平方的结果

\text{SINR}_{k\to j}等于分子部分\operatorname{Tr}(\mathrm{H}_j\mathbf{W}_k)与分母中的累加项以及噪声方差\sigma_n^2之和。

经由通分操作，则上式可进一步推导得出：
\begin{aligned} & R_{k→j}=log_2(σ² + ∑_{i∈K, i≥k} Tr(H_j W_i)) \\ & underbrace{-log_2(σ² + ∑_{i∈K, i>k} Tr(H_j W_i))}_{F_{j,k}} ≥ γ_k. \\ & \\ & 结合迹运算特性，则有 R_{k→j}=log_2(...) - F_{j,k} ≥ γ_k$. \ & 结束推导过程。

上式中第一项为凸函数，则\sum_{i \in \mathcal{K}, i \geq k} \operatorname{Tr}(\mathbf{H}_j\mathbf{W}_i)是关于\mathbf{W}_i的一个线性组合具有凸性的性质；而F_{j,k}是非凸函数的原因在于其前部带有负号的作用

收敛速度与目标函数值之间存在权衡关系，在优化过程中选择不同的惩罚项权重系数

• |\mathrm{h}_k^H\mathrm{w}_i|^2=\mathrm{Tr}(\mathbf{H}_k\mathbf{W}_i) 下面解释这个式子，

定义，\mathrm{Tr}(\mathbf{H}_k\mathbf{W}_i)=\mathrm{Tr}(\mathrm{h}_k(\mathrm{h}_k^H\mathrm{w}_i)\mathrm{w}_i^H)
标量提出来和迹的性质\mathrm{Tr}(\mathbf{AB})=\mathrm{Tr}(\mathbf{BA})，\mathrm{Tr}(\mathbf{H}_k\mathbf{W}_i)=\mathrm{h}_k^H\mathrm{w}_i\cdot\mathrm{Tr}(\mathrm{w}_i^H\mathrm{h}_k)
标量的迹就是本身，\mathrm{Tr}(\mathbf{H}_{k}\mathbf{W}_{i})=\mathrm{h}_{k}^{H}\mathrm{w}_{i}\cdot\mathrm{w}_{i}^{H}\mathrm{h}_{k}=(\mathrm{h}_{k}^{H}\mathrm{w}_{i})^{2}

秩为1的约束可以用the semidefinite relaxation (SDR)省略

PRML Eq10.125

PRML Eq10.125

\frac{d}{dx} \log_2(x) = \frac{1}{x \ln(2)}

对 -\log_2(x) 在 x = x_n 处的泰勒展开为：-\log_2(x) \approx -\log_2(x_n) - \frac{1}{x_n \ln 2} (x - x_n)

\begin{aligned} &\text{-}\log_2 (\sigma_n^2 + &\sum_{i ∈ 𝒦, i > k} &\operatorname*{Tr}(H_j W_i)) ≈ \\ &-\log_2 (\sigma_n^2 + &\sum_{i ∈ 𝒦, i > k} &\operatorname*{Tr}(H_j W_i^n)) \\ &- [∑_{i ∈ 𝒦, i > k} Tr(H_j (W_i - W_i^n))] / [(σ_n² + ∑_{i ∈ 𝒦, i > k} Tr(H_j W_iⁿ)) · ln ₂] \end{aligned}

计算信噪比

10.^((-para.noise - path_loss)/10)：将以dB表示的电平转换为线性比例。10^(x/10) 用于将以dB表示的电平转换成线性的数值

在电功率领域中，我们通常使用分贝（dB）这一对数单位来量化能量或信号的强弱程度（想了解dBW的含义吗？）。具体而言：

• dBm 是基于对数尺度衡量功率相对于1毫瓦(milliwatt, mW)大小的一种单位。
• 转换公式如下：
  P\:(\mathrm{dBm})=10\log_{10}\left(\frac{P\:(\mathrm{mW})}{1\:\mathrm{mW}}\right)
• 典型换算关系：
  1 \:\mathrm{W} = 30 \:\mathrm{dBm}
  例如：
• 30 \:\mathrm{dBm} \leftrightarrow 3.684 \:\mathrm{V/m}（此处补充了具体数值）
• 2 \:\mathrm{dB} 的增益意味着信号强度增加了约4倍
• -6 \:\mathrm{dB} 的衰减意味着信号强度减少了约25%
• 典型换算关系：
  0 \:\mathrm{dBm} = 1 \:\mathrm{mW}
  例如：
  -30 \:\mathrm{dB}= 4.999999997\times {1e^{-8}}

dB是一种对数尺度指标，在表现功率相对变化方面具有重要作用。增益为+3dB时表明输出功率是输入的两倍；当输出功率降至-3dB时，则相当于变为原来的一半；而0dB则意味着输出信号与输入信号的强度完全相同。此外，在声学领域中常用分贝（decibel, dB）来量化声音强弱的程度。分贝 (dB) 被定义为衡量功率比值的对数值。具体而言，在工程应用中通常采用以下公式计算两个功率值P₁和P₂之间的差异：\mathrm{dB}=10\log_{10}\left(\frac{P_2}{P_1}\right)。其中P_{1}代表初始功率水平而P_{2}代表经过处理后的测量值。

d 表示分贝（decibel）
m 表示毫瓦（milliwatt）

将功率转换为瓦特单位：
发送功率P_t=20 dBm；
噪声功率值\sigma_n^2=-120 dBm；
其计算公式如下：

P\left(\text{瓦特}\right)= 10^{\frac{P\left(\mathrm{dBm}\right)}{10}} \times 1e-3

2. 计算发射功率$P_{t}$在瓦特中：  

P_t=10^{\frac{20}{10}}\times10^{-3}\text{瓦特}=10^2\times10^{-3}

3. 计算噪声功率$\sigma_n^2$在瓦特中：  

\sigma_n^2=10^{\frac{-120}{10}}\times10^{-3}\text{瓦特}=10^{-12}\times10^{-3}=10^{-15}

4. 计算信噪比(SNR):  

\text{SNR}=\dfrac{P_t}{\sigma_n^2}=\dfrac{0.1}{10^{-15}}=10^{14}

在信号处理中将信噪比（SNR）从线性单位转换为对数单位dB时：
\mathrm{SNR~(dB)}=10\log_{10}(\mathrm{SNR})=10\log_{10}(10^{14})=10\times14=140\:\mathrm{dB}
当系统在给定发射功率P_t=20 dBm以及噪声功率\sigma_n^2=-120 dBm的情况下进行工作时：
因此，在这种情况下，
信噪比\mathrm{SNR}被计算得出为140 dB。

3-code

此两类优化问题必须提前告知CVX平台

二次规划：其所有特征值必须是非负的。
二阶锥SOCP(second-order cone program)：对于每个二阶锥约束,其前部分元素的平方和小于等于最后一个元素。
SDP(semidefinite program): 半正定规划问题中的决策变量为对称矩阵。
如果没有cvx begin sdp, X >= 0表示所有矩阵元素非负。

Solver precision设置求解精度

SCA_algorithm

beampattern

    function [pattern] = beampattern(R,M,theta_degree)
    %The beampattern of ULA in a DFRC BS with seperate depolyment
    %  [pattern] = beampattern(P,Rx,Mc,Mr,theta_degree)
    %Inputs:
    %   R: 发射信号的协方差矩阵。这个矩阵表示信号在不同天线间的相关性和功率分布。
    %   M: 通信天线的数量。这个参数决定了天线阵列的大小。
    %   theta_degree: 角度范围，以度为单位。这个范围指定了你希望计算波束图的方向。
    %Outputs:
    %   pattern: beampattern
    
    theta = theta_degree*pi/180;
    pattern = zeros(length(theta),1);
    % C = [Rx, zeros(Ntr, Ntc); zeros(Ntr,Ntc), P*P'];
    
    for i = 1:length(theta)
    t = theta(i);
    a = ULA_func(t,M); % 通过 ULA_func 函数计算方向向量 a
    pattern(i) = real(a'*R*a); % 计算的是信号在该方向上的功率响应，并取实部来表示实际的功率。
    end
    
    % pattern = 10*log10(pattern);
    end
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    代码解释

rate_NOMA

    function [R] = rate_NOMA(para, W, H)
    % 计算NOMA-DFRC系统的可达速率
    %  [R] = rate_NOMA(para, W, H)
    % 输入:
    %   para: 系统参数
    %   W: 波束形成矩阵
    %   H: 信道增益矩阵
    % 输出:
    %   R: 可达速率
    % 初始化速率矩阵为正无穷
    R = inf * ones(para.K, para.K);
    
    % 遍历每个用户组中的用户
    for k = 1:para.K-1
       W_k = W(:,:,k);  % 获取第k个用户的波束形成矩阵
       
       % 对第k个用户的信号在第j个用户处进行解码
       for j = k:para.K
       H_j = H(:,:,j);  % 获取第j个用户的信道矩阵
       
       % 计算干扰
       I = 0;
       for i = k+1:para.K
           W_i = W(:,:,i);  % 获取第i个用户的波束形成矩阵
           I = I + real(trace(H_j * W_i));  % 计算干扰信号的总功率
       end     
       % 计算第k个用户在第j个用户处的可达速率
       R(k,j) = log2( 1 + real(trace(H_j * W_k)) / (I + para.n) );
       end
    end
    
    % 处理最后一个用户的速率
    W_K = W(:,:,end); 
    H_K = H(:,:,end);
    R(end,end) = log2( 1 + real(trace(H_K * W_K)) );
    
    end
    
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    代码解释

initial_rate函数

因为缺乏明确的目标函数定义,CVX默认处理为一个可行问题.这意味着该系统将寻求一组能够满足所有约束条件的变量值,从而无需针对任何特定目标函数进行优化.

    function [ineqn] = initial_rate(para, H, W)
    % 为 SCA 算法中的速率初始化设置不等式约束
    % 该函数用于根据系统参数和信道增益，初始化 SCA 算法中的速率不等式约束。
    % 输入:
    %   para: 系统参数，包括天线数量、用户数量、噪声功率等
    %   H: 信道增益矩阵，大小为 (N, N, K)，其中 N 是天线数，K 是用户数
    %   W: 波束形成矩阵，大小为 (N, N, K)
    % 输出:
    %   ineqn: 包含速率约束的不等式集合，是 CVX 变量的数组
    % 日期: 03/10/2021
    % 作者: Zhaolin Wang
    
    ineqn = [];  % 初始化不等式约束数组
    
    % 对于每个用户 k（除了最后一个用户 K）计算速率不等式约束
    for k = 1:para.K-1
    W_k = W(:,:,k);  % 当前用户 k 的波束形成矩阵 W_k
    
    % 对于用户 k 至最后一个用户 j 计算信道 H_j 和不等式
    for j = k:para.K
        H_j = H(:,:,j);  % 当前信道矩阵 H_j
        
        % 干扰项初始化
        I = 0;
        
        % 计算来自用户 k+1 至 K 的干扰项
        for i = k+1:para.K
            W_i = W(:,:,i);  % 干扰用户 i 的波束形成矩阵 W_i
            I = I + real(trace(H_j * W_i));  % 计算干扰的真实部分并累加
        end
        
        % 将当前用户 k 对用户 j 的速率不等式添加到 ineqn 数组中
        % 不等式形式：-real(trace(H_j*W_k)) + (2^para.R_min - 1) * (I + para.n)
        % 解释: 速率约束确保用户 j 在给定干扰条件下达到最小速率 R_min
        ineqn = [ineqn, -real(trace(H_j * W_k)) + (2^para.R_min - 1) * (I + para.n)];
    end   
    end
    
    % 最后一个用户 K 的特殊处理
    H_K = H(:,:,end);  % 获取最后一个用户 K 的信道增益矩阵
    W_K = W(:,:,end);  % 获取最后一个用户 K 的波束形成矩阵
    ineqn = [ineqn, -real(trace(H_K * W_K)) + (2^para.R_min - 1) * para.n];  % 添加速率不等式
    
    end
    
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    代码解释

initial_point函数

    %% initialize - 初始化波束形成矩阵的初始值
    % 该函数用于在优化过程中提供初始波束形成矩阵 W。
    % 输入:
    %   para - 系统参数，包括天线数量、用户数量、发射功率等
    %   H    - 信道矩阵，大小为 (N, N, K)，其中 N 是天线数，K 是用户数
    % 输出:
    %   W         - 初始化后的波束形成矩阵，大小为 (N, N, K)
    %   cvx_status - CVX 优化状态，用于判断优化是否成功
    
    function [W, cvx_status] = initial_point(para, H)
    
    cvx_begin quiet  % 启动 CVX 优化，quiet 选项用于抑制输出信息
        
        % 定义复数波束形成矩阵变量 W，大小为 (N, N, K)
        variable W(para.N, para.N, para.K) complex
    
        % 计算总波束形成矩阵的协方差矩阵 Rw
        Rw = sum(W,3);
    
        % 计算探测功率 P
        [P] = probing_power(para, Rw);
    
        % 计算初始速率不等式 ineqn
        [ineqn] = initial_rate(para, H, W);
    
        % 计算交叉相关性约束 corr
        [corr] = cross_correlation(para, Rw);
    
        % 约束1：每个用户的波束形成矩阵 W(:,:,k) 必须是 Hermitian 半正定矩阵
        for k = 1:para.K
            W(:,:,k) == hermitian_semidefinite(para.N);
        end
    
        % 约束2：确保所有速率不等式 ineqn(i) <= 0
        for i = 1:length(ineqn)
            ineqn(i) <= 0;
        end
    
        % 约束3：相邻角度探测功率之差在 [-10, 10] 之间
        M = length(para.theta);  % 参数 theta 表示角度的集合
        for k = 1:M
            for p = k+1:M
                P(k) - P(p) <= 10;
                P(k) - P(p) >= -10;
            end
        end
    
        % 约束4：波束形成矩阵的对角线元素等于发射功率 Pt 平均分配到 N 个天线上的功率
        diag(Rw) == para.Pt / para.N;
    
        % 约束5：交叉相关性约束小于等于给定的最小值
        corr <= para.correlation_min;
    
    cvx_end  % 结束 CVX 优化，并返回优化结果
    
    end
    
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    代码解释

SCA_SINR

    function [sca_all] = SCA_SINR(para, H, W, W_n, r)
    % SCA_SINR - 针对信干噪比 (SINR) 约束的逐步凸逼近 (SCA) 方法实现。
    %  [sca_all] = SCA_SINR(para, H, W, W_n, r)
    %
    % 输入:
    %   para - 系统参数，包括天线数量、用户数量等
    %   H    - 信道增益矩阵，大小为 (N, N, K)，其中 N 是天线数，K 是用户数
    %   W    - 当前步骤中的波束形成矩阵 (cvx 变量)，大小为 (N, N, K)
    %   W_n  - 上一步中的波束形成矩阵
    %   r    - 当前步骤中的辅助变量 (cvx 变量)，表示 SINR 的最小值
    %
    % 输出:
    %   sca_all - 为每个用户生成的 SCA 凸 cvx 变量，用于构造 SINR 约束
    %
    % 日期: 2021年6月10日
    % 作者: Zhaolin Wang
    
    sca_all = [];  % 初始化用于存储所有 SCA 变量的数组
    
    % 遍历前 K-1 个用户，逐个用户计算 SINR 的凸近似
    for k = 1:para.K-1
    W_k = W(:,:,k);  % 提取第 k 个用户的波束形成矩阵
    r_k = r(k);      % 提取第 k 个用户的 SINR 辅助变量
    
    % 解码过程
    for j = k:para.K
        H_j = H(:,:,j);  % 提取第 j 个用户的信道矩阵
        
        % 计算干扰项 (interference)
        I = 0; I_n = 0; I_diff = 0;  % 初始化干扰及其相关的变量
        for i = k+1:para.K
            W_i = W(:,:,i);        % 第 i 个用户的波束形成矩阵
            W_n_i = W_n(:,:,i);    % 上一步骤中第 i 个用户的波束形成矩阵
            I = I + real(trace(H_j*W_i));  % 当前干扰项的累计
            I_n = I_n + real(trace(H_j*W_n_i));  % 上一步骤干扰项的累计
            I_diff = I_diff + real(trace(H_j * (W_i - W_n_i)));  % 干扰项的差异累计
        end
        
        % 计算 j 用户对 k 用户的 SCA 凸近似项
        % F_j_k 是通过 Taylor 展开得到的第一项近似
        F_j_k = - log(para.n + I_n)/log(2) - I_diff / (para.n + I_n) / log(2);
        
        % 计算当前 SCA 约束
        % sca 是对 log(1 + SINR) 函数的凸近似表示
        sca = log(para.n + real(trace(H_j * W_k)) + I)/log(2) + F_j_k - r_k;
        
        % 将当前 SCA 约束添加到 sca_all 中
        sca_all = [sca_all, sca];
    end   
    end
    
    % 对最后一个用户 (第 K 个用户) 计算 SINR 的凸近似
    H_K = H(:,:,end);  % 提取最后一个用户的信道矩阵
    W_K = W(:,:,end);  % 提取最后一个用户的波束形成矩阵
    r_K = r(end);      % 提取最后一个用户的 SINR 辅助变量
    
    % 计算最后一个用户的 SCA 约束，并添加到 sca_all 中
    % 这里使用了 log(1 + SINR) 的直接表示，而不需要凸近似，因为最后一个用户没有干扰项
    sca_all = [sca_all, log( 1 + real(trace(H_K*W_K))/para.n )/log(2)-r_K];
    
    end
    
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    代码解释

SCA_step

    % 开始CVX优化环境，'quiet'选项表示抑制CVX输出信息
    cvx_begin quiet
    
    % 定义一个复数变量 W，大小为(N x N x K)，表示为每个用户分配的波束成形矩阵
    variable W(para.N,para.N,para.K) complex
    
    % 定义一个变量 r，大小为(K x 1)，表示每个用户的信噪比
    variable r(para.K, 1)
    
    % 将W沿第三维度相加，得到总的波束成形矩阵 Rw
    Rw = sum(W, 3);
    
    % 计算探测功率 P，输入参数为系统参数 para 和波束成形矩阵 Rw
    [P] = probing_power(para, Rw);
    
    % 计算信噪比的近似值 sca，输入参数为系统参数 para、信道矩阵 H、波束成形矩阵 W、噪声矩阵 W_n 和信噪比变量 r
    [sca] = SCA_SINR(para, H, W, W_n, r);
    
    % 计算交叉相关性 corr，输入参数为系统参数 para 和波束成形矩阵 Rw
    [corr] = cross_correlation(para, Rw);
    
    % 约束条件
    
    % 1. 对于每个用户 k，波束成形矩阵 W(:,:,k) 必须是 N 阶的 Hermitian 半正定矩阵
    for k = 1:para.K
        W(:,:,k) == hermitian_semidefinite(para.N);
        % 2. 每个用户的信噪比 r(k) 必须大于等于最小需求 R_min
        r(k) >= para.R_min;
    end
    
    % 3. 确保每个用户的信噪比近似值 sca(i) 非负
    for i = 1:length(sca)
        sca(i) >= 0;
    end
    
    % 4. 控制不同方向的探测功率 P(k) 之间的差异
    M = length(para.theta); % 方向角 theta 的数量
    for k = 1:M
        for p = k+1:M
            % 控制探测功率的差异在 [-10, 10] 范围内
            P(k) - P(p) <= 10;
            P(k) - P(p) >= -10;
        end
    end
    
    % 5. 控制总波束成形矩阵 Rw 的对角元素，使其与每个天线的平均功率一致
    diag(Rw) == para.Pt / para.N;
    
    % 6. 控制交叉相关性小于等于设定的最小相关性门限 correlation_min
    corr <= para.correlation_min;
    
    % 目标函数
    
    % 计算目标函数 f，输入参数为波束成形矩阵 W、噪声矩阵 W_n、信噪比 r、探测功率 P、用户数 K、
    % 通信与探测的权重系数 rho_c 和 rho_r，以及系统效率参数 eta
    [f] = obj_func(W, W_n, r, P, para.K, rho_c, rho_r, eta);
    
    % 最小化目标函数 f
    minimize(f);
    
    % 结束CVX优化环境
    cvx_end
    
    
    
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
    
    代码解释

total throughput（总吞吐量）– 25, total sensing power（总感知功率）= 77, target function value（目标函数值）= -1, penalty coefficient（惩罚系数）= -8 × 10^-5；
在第③次测试中，
total throughput（总吞吐量）– = -,
total sensing power（总感知功率）= -,
target function value（目标函数值）= -, penalty coefficient（惩罚系数）= -；
测试次数④中，
total throughput（总吞吐量）– = -,
total sensing power（总感知功率）= -,
target function value（目标函数值）= -, penalty coefficient（惩罚系数）= -

sum rate – 16.3811 / sum probing power-- 759.7147 / objective value-923.5164
penalty term – 6.4237e-05

sum rate – 16.3332 / sum probing power-- 772.9037 / objective value-936.2359
penalty term – 2.5672e-05
sum rate – 16.0283 / sum probing power-- 773.2016 / objective value-933.4848
penalty term – 6.1369e-05
sum rate – 15.9629 / sum probing power-- 772.6916 / objective value-932.3208
penalty term – 8.2638e-05
sum rate – 15.2548 / sum probing power-- 773.0922 / objective value-925.6404
penalty term – 8.4755e-05

sum rate – 15 / sum probing power-- 697.2146 / objective value-847.1164
penalty term – 1.0831e-05

sum rate – 13.6034 / sum probing power-- 770.2418 / objective value-906.276
penalty term – 5.393e-05

该参数达到了...的表现：
数值为8.7554的总速率，
数值为772.5604的总探bing能力，
数值为860.1141的目标函数值，
惩罚项值为1.615e-06。
第二段：
该参数达到了...的表现：
数值为
的总速率,
数值
的总探bing能力,
数值
的目标函数值,
惩罚项
。

该值对应的线性SINR值为10^{\frac{e^{k-1}}{10}} dB。
该值对应的信道容量C为\log_2(1 + 2^k - 1)。

4-附录

4.1- the quadratic form of the covariance matrix是非凸的

关于协方差矩阵的二次形式通常表现为非凸性质。这种现象可以从以下几个角度进行分析和解释：

本节将介绍二次型及其相关性质。考虑 \bm{\mathrm{\ W}} 作为一个向量或矩阵，则该二次型可表示为 f(\bm{\mathrm{\ W}}) = \bm{\mathrm{\ W^{H}}} \cdot \bm{\mathrm{\ R}} \cdot \bm{\mathrm{\ W}}。其中符号\bm{\mathrm{\ W^{H}}}表示\bm{\mathrm{\ W}}的共轭转置，并且\bm{\mathrm{\ R}}被定义为协方差矩阵的一个实例；特别地，在实际应用中我们通常假设其为半正定或正定矩阵以确保该函数具有良好的优化特性。

非凸性分析
在优化问题中，在某个范围内Hessian矩阵（即二阶导数矩阵）不是正定的情况下，则称该函数为非凸函数。协方差矩阵对应的二次型可能出现非凸性的情况，请问这是什么情况？原因如下：

Hessian矩阵存在非正定性质 :
针对二次型表达式 \mathbf{W}^H \mathbf{R} \mathbf{W}, 其对应的Hessian矩阵计算结果为 2\mathbf{R}. 然而协方差矩阵 \mathbf{R} 原本就是半正定或全正定的, 因此优化问题出现非凸性主要源于约束条件或其他变量特性的影响, 而与二次型自身并无直接关联.

非凸约束的影响 :
在实际应用场景中, 二次型优化问题往往伴随着多种非凸约束条件（例如, 秩约束、非线性约束等）。这些约束条件可能使整个优化问题成为一个典型的非凸优化问题。

3. 非凸优化问题的难度
   相较于凸优化问题而言，在复杂性上来说是非凸优化问题更为突出的问题。这类问题是由于其可能存在多个局部最优点而难以求得真正的全局最优点。在实际应用中，在面对这类复杂情况时, 传统的方法(例如梯度下降法)往往只能得到一个局部最优点, 而无法保证得到的是真正的全局最优点。

总结
二次型的非凸性主要表现为涉及协方差矩阵时的情况。当优化问题中存在由协方差矩阵带来的复杂性或附加了非凸约束以及其他类型的复杂结构时，则会导致整个优化过程呈现出高度的不确定性与挑战性，并最终形成一个全局性的NP难（NP-Hard）组合优化问题模型。这类复杂的问题普遍存在于实际工程领域，并且往往需要利用特定类型的非线性规划方法来解决。

4.2-感知目标信号的互相关

在 MIMO 雷达系统中，在降低发射信号之间的相互关联度（cross-correlation）方面存在几个主要原因。这些因素包括与目标识别能力相关的指标、图像清晰度（即分辨率）、抗干扰性能以及整体系统效率提升密切相关。因此，在实际设计和应用中应综合考虑这些影响因素。

1. 提高目标检测能力

• 降低多径干扰影响：当发射信号之间的互相关值较小，在不同目标或路径返回的回波信号更容易被分辨并识别。
  *MIMO雷达系统增强分辨率能力：通过使用低互相关特性的信号，在MIMO雷达中能够有效提升目标的角度分辨率和距离分辨率。即使多个目标处于相似的方向或距离范围内，在雷达系统中也能实现更加有效的区分与分离。

  2. 减少自干扰和旁瓣
  • 自干扰的最小化 ：互相关较低的信号减少了自干扰现象（即同一系统的信号干扰自身接收过程），从而提高了雷达接收的信号质量。自干扰会降低信噪比（SNR），影响系统的检测和测量精度。
  • 旁瓣抑制 ：低互相关的信号有助于降低雷达天线方向图的旁瓣效应，从而减少旁瓣中的误检测。旁瓣的存在可能导致虚假目标的出现，影响目标检测的可靠性。
  3. 提升多目标检测的能力

• 降低各目标间的相互干扰 ：在存在多个探测目标的环境中，在不同发射信号间的互相关性较低的情况下，在接收机所观测到的回波信号间的相互干扰会得到降低。此时雷达接收机所观测到的回波信号间的相对相位差异会更加显著，在多目标检测任务中能够更好地分辨出各个独立的目标，并且提高整体探测性能水平。

  4. 支持更复杂的信号处理技术

• 优化信号解码 ：在接收端的信号处理系统中，具有低互相关特性的信号在解码过程中表现出较高的可识别性。特别是在采用先进的波束形成技术和基于距离-多普勒效应的处理方法时，在复杂噪声环境下依然能够有效识别出目标信号参数并实现精准定位。

  5. 增强抗干扰能力

• 抗干扰技术：在雷达对抗干扰场景中采用具有低互相关特性的信号能够有效抵御有意或无意产生的电磁干扰；这是因为这些信号经过解相关处理后能够有效去除不想要的干扰。

综上所述，在MIMO雷达系统中减小发射信号间的自相关度能够将该系统的检测性能、分辨能力等关键指标均得到显著提升，并进而进一步提升了该系统的整体效能。