简介这是一份基于MATLAB的波束成形、DOA估计与RIS可重构智能表面技术仿真资源面向通信、电子信息、数学等专业学生的课程设计、期末大作业与毕业设计场景。代码采用参数化编程参数可灵活调整注释清晰并附带可直接运行的案例数据能帮助使用者快速掌握阵列信号处理与智能超表面辅助通信的核心原理。压缩包共19个文件以18个.m源码文件为主另含1个.md说明文档整体仅39KB轻量精简、便于本地部署与二次开发。目前已有92人学习或下载适合正在学习相关课程或计划开展科研复现的学生参考通过阅读代码与案例可系统理解波束成形、DOA估计及RIS辅助通信等关键步骤也可在此基础上修改参数、替换场景完成自己的仿真实验。1. 一个 ZIP 里同时出现三个词先别急着解压拿到“波束成形-DOA-RIS.zip”这个包多数人的第一反应是解压、找 .m 或 .py 文件直接运行。但真正决定这个包价值的是三个英文缩写之间的依存关系波束成形决定能量怎么打出去DOA 估计决定朝哪个方向打RIS 则是中途用来改变电磁传播路径的一张“可编程反射镜”。三者放在同一个压缩包里通常不是三个独立 demo 的简单合并而是一条“先探测信道 → 估计用户方位 → 设计收发权值 → 协同调 RIS 相位”的链路。这篇文章按我拿到此类仿真包时会走的路径来讲先固定阵列模型再做 DOA再引入 RIS 联合设计最后落到验证方法。无论包里是 MATLAB 还是 Python 工程读完之后你至少能判断这段代码在仿什么、参数改哪里、结果是否可信。2. 阵列模型与波束成形把坐标系和导向矢量先钉死2.1 均匀线阵的导向矢量窄带、远场、平面波是默认前提在波束成形和 DOA 估计的代码里出现最多的结构是均匀线阵ULA。设阵元数为 (M)阵元间距为 (d)信号入射角为 (\theta)采用窄带、远场、平面波假设时第 (m) 个阵元相对参考阵元的传播时延为 (m d \sin\theta / c)。用载波频率归一化后导向矢量写成[ \mathbf{a}(\theta) \left[1,\ e^{j 2\pi (d/\lambda) \sin\theta},\ \ldots,\ e^{j 2\pi (M-1)(d/\lambda) \sin\theta}\right]^T ](d/\lambda) 是最关键的归一化参数。取 0.5半波长时阵列无栅瓣同时角度分辨率达到常用配置的平衡点。(d/\lambda) 小于 0.4 时主瓣变宽大于 0.5 时视场边缘出现栅瓣DOA 估计会出现模糊峰。波束成形与 DOA 这两个模块的精度上限最早就由这个比值决定。坐标系约定同样值得花三十秒确认。多数代码以阵列法线为 0°一侧为正、另一侧为负少数实现用端射方向做参考。若把 (\theta) 直接代入 (\sin\theta)两侧符号换掉波束图会左右翻转MUSIC 谱的峰也会出现在错误象限。拿到包的第一步先找 angle_grid 的定义确认角度范围和步进——这个动作能避免后面几十分钟的无效调试。2.2 从匹配滤波器到 MVDR两种权值设计的取舍常规相控阵直接取 (\mathbf{w} \mathbf{a}(\theta_0))相当于对目标方向的信号做相干累加实现简单但对干扰没有任何抑制能力。MVDR最小方差无失真响应在目标方向保持增益 1同时最小化输出总功率从而把干扰方向“压出凹口”。其闭式解为[ \mathbf{w} \frac{\mathbf{R}^{-1} \mathbf{a}(\theta_0)}{\mathbf{a}(\theta_0)^H \mathbf{R}^{-1} \mathbf{a}(\theta_0)} ]这里 (\mathbf{R}) 是阵列接收数据的协方差矩阵。实际工程中快拍数不足时样本协方差矩阵求逆会放大噪声通常要加对角加载(\mathbf{R} \lambda_0 \mathbf{I})。加载系数 (\lambda_0) 一般取 (\mathbf{R}) 对角线均值乘以 (10^{-2}) 到 (10^{-1})。过小失去稳定作用过大会退化成匹配滤波器。代码包里如果只有 (\mathbf{a}(\theta_0)) 直接加权而没有矩阵求逆那就是常规相控阵谈不上自适应抗干扰。2.3 用 Python 跑通最小波束成形闭环并检查主瓣指向import numpy as np def steering_vector(M, theta_deg, d_over_lambda0.5): 均匀线阵导向矢量theta_deg 是相对阵列法线的角度 theta np.deg2rad(theta_deg) m np.arange(M) return np.exp(1j * 2 * np.pi * d_over_lambda * m * np.sin(theta)) def mvdr_weights(Rxx, a): MVDR 权值带对角加载防止 R 奇异 lam 1e-2 * np.trace(Rxx) / Rxx.shape[0] Rinv np.linalg.inv(Rxx lam * np.eye(Rxx.shape[0])) return Rinv a / (a.conj().T Rinv a) M, N 8, 500 theta0, snr_db 10, 10 s (np.random.randn(N) 1j*np.random.randn(N)) / np.sqrt(2) a0 steering_vector(M, theta0) noise_power 10**(-snr_db/10) X np.outer(a0, s) np.sqrt(noise_power) * ( np.random.randn(M, N) 1j*np.random.randn(M, N)) / np.sqrt(2) Rxx (X X.conj().T) / N w mvdr_weights(Rxx, a0) angles np.linspace(-90, 90, 361) pattern [abs(w.conj() steering_vector(M, ang))**2 for ang in angles] peak angles[np.argmax(pattern)] print(f主瓣峰值: {peak}°, 期望: {theta0}°)逻辑说明先按“单信号加高斯白噪声”生成 8 阵元、500 快拍的数据样本协方差 Rxx 是 (8\times8) 的复数矩阵。MVDR 部分做了一次对角加载加载系数取对角线均值的 1%对 10 dB 信噪比场景足够稳定。归一化由复数除法完成确保期望方向增益为 1。参数可以按场景来改阵元数 M 决定空间自由度M 越大旁瓣越低N 是快拍数从 500 降到 50 时Rxx 估计误差变大主瓣指向偏移会更明显这也是信道估计类项目里最常见的坑。把 snr_db 改成 -5再看主瓣峰值位置能直观理解“低信噪比时波束成形精度不再由孔径决定而由协方差估计误差决定”。3. DOA 估计与 MUSIC 的实用参数子空间维度是第一个检查点3.1 为什么 MUSIC 能做出超分辨率MUSIC 类方法的基本前提是信号子空间与噪声子空间正交。对阵列协方差 (\mathbf{R}) 做特征分解后大特征值对应的特征向量张成信号子空间余下特征向量张成噪声子空间 (\mathbf{E}_n)。由于导向矢量落在信号子空间内它与噪声子空间正交因此在真实入射角附近满足 (\mathbf{a}(\theta)^H \mathbf{E}_n \approx \mathbf{0})。空间谱定义为该内积模长的倒数[ P_{\text{MUSIC}}(\theta) \frac{1}{\mathbf{a}(\theta)^H \mathbf{E}_n \mathbf{E}_n^H \mathbf{a}(\theta)} ]谱峰位置即 DOA 估计结果。超分辨能力的来源在于特征分解能区分角度差小于阵列瑞利波束宽度的两个信号源只要快拍数和信噪比足够。MUSIC 谱不是功率谱它没有物理单位谱峰高度没有绝对意义只有相对意义。注意一个微妙点MUSIC 能分辨多少个源取决于样本协方差矩阵能支持多少个维度而不是真实存在多少个来波。信源数 K 估多了本应属于噪声子空间的向量被划入信号子空间谱会漏峰K 估少了伪峰出现。常见做法是用 MDL 或 AIC 准则估计信源数再传给 MUSIC。代码包里如果直接写死 K1而仿真场景里有两个用户结果必然失真。3.2 SubspaceNet 的思路用网络拟合子空间投影如果做 5G-A 或通感一体化这类场景快拍数经常少于阵元数此时样本协方差不满秩特征分解的性能边界迅速退化。SubspaceNet 这类方法的思路是用深度网络直接学习从阵列数据到信号子空间投影的映射输入可以是复协方差矩阵的实虚部堆叠监督信号则由高快拍数下的特征分解产生。训练完成后推理时避开显式特征分解在低快拍时比传统 MUSIC 更稳。放在 ZIP 包里的位置通常是 DOA 估计模块的一个可选替换传统代码里是 np.linalg.eigh替换后是一个网络推理函数。实际使用中要小心训练场景与测试场景的天线结构不一致——子空间形状和阵元布局强绑定换阵列就要重新训练或微调。3.3 一个不依赖外部库的 MUSIC 最小实现def music_spectrum(X, K, angles): X: M×N 阵列数据; K: 信源个数; angles: 扫描角度列表 M X.shape[0] Rxx (X X.conj().T) / X.shape[1] eigvals, eigvecs np.linalg.eigh(Rxx) En eigvecs[:, :M-K] # 噪声子空间特征值升序排列取前 M-K 列 spec [] for ang in angles: a steering_vector(M, ang) spec.append(1.0 / np.abs(a.conj() En En.conj().T a)) return np.array(spec) angles np.linspace(-90, 90, 361) spec music_spectrum(X, K1, anglesangles) est angles[np.argmax(spec)] print(fDOA 估计结果: {est}°, 真实来波: {theta0}°)逻辑说明np.linalg.eigh 返回的特征值按升序排列因此噪声子空间取前 M-K 列。空间谱每个角度独立计算一次倒数扫描 361 个角度时开销不大工程实现里可以用矩阵一次性完成计算。K 的选取直接影响谱峰质量可以在 -10 dB 信噪比下把 K 改成 2观察伪峰是否出现。把 M、快拍 N、信噪比三个参数串起来看会更直观N 小于 M 时样本协方差秩亏特征分解会把噪声能量摊到后 M-K 个特征值上谱峰抖动明显。工程上我会先保证 N ≥ 4M再谈超分辨率。4. RIS 加入后模型多了一层反射级联信道与联合设计4.1 级联信道模型和反射系数的影响RIS 由 N 个反射单元组成每个单元相当于一个可调相移器反射系数写为 (\beta_n e^{j\varphi_n})(\beta_n) 为幅度理想情况取 1(\varphi_n) 为相移。在远场窄带假设下基站到用户的等效信道由直射径 (\mathbf{h}d)、基站-RIS 信道 (\mathbf{G})(N\times M) 矩阵、RIS-用户信道 (\mathbf{h}{ru})N 维向量三部分组成级联信道表达式为[ \mathbf{h}_{\text{eff}} \mathbf{h}d \mathbf{h}{ru}^H \boldsymbol{\Phi} \mathbf{G} ]其中 (\boldsymbol{\Phi} \mathrm{diag}(\beta_1 e^{j\varphi_1}, \ldots, \beta_N e^{j\varphi_N}))。RIS 本身不放大信号它通过调整相位使来自不同反射单元的径在接收端相干叠加形成类似波束成形的“空间聚焦”。这也是为什么 RIS 相关工作里频繁出现“被动波束成形”这个词——它不消耗发射功率却改变信道增益。DOA 与 RIS 的结合点在于要对用户做定向波束成形先要知道用户方向当直射径被遮挡时波达方向信息只能通过 RIS 反射链路获得。因此这类代码包里会先固定一组相移做探测用 DOA 估计出反射来波方向再据此优化下行权值。这里的 DOA 估计对象是“经过 RIS 反射后的虚拟路径方向”不是直接来波方向。4.2 相移设计和两种常见联合优化策略理想相移是无约束最优解(\varphi_n \arg(\mathbf{h}_{ru,n}) - \arg([\mathbf{G}\mathbf{a}(\theta)]_n))即把每个反射单元的用户信道相位与到达相位对齐。但每个单元相移是量化的例如 2-bit 移相时 (\varphi_n \in {0, \pi/2, \pi, 3\pi/2})量化误差会让波束成形增益下降1-bit 量化损失约 3.9 dB2-bit 损失约 0.9 dB。代码包中看到相位取整操作就是这种离散化。下面给出理想相移计算并对齐的片段N_ris, M 32, 8 theta_user 30 G (np.random.randn(N_ris, M) 1j*np.random.randn(N_ris, M)) / np.sqrt(2) h_ru (np.random.randn(N_ris) 1j*np.random.randn(N_ris)) / np.sqrt(2) g_beamformed G steering_vector(M, theta_user) target_phase np.angle(h_ru * g_beamformed) phi_bits 2 phi_levels 2**phi_bits step 2*np.pi / phi_levels phi_quant np.round(target_phase / step) * step Phi np.diag(np.exp(1j * phi_quant)) cascaded np.abs(h_ru.conj() Phi G steering_vector(M, theta_user)) perfect np.abs(h_ru.conj() np.diag(np.exp(1j * target_phase)) G steering_vector(M, theta_user)) print(f量化后级联增益: {cascaded:.3f}, 理想增益: {perfect:.3f})逻辑说明先计算每个反射单元的到达复增益 g_beamformed再与用户信道 h_ru 逐元素相乘取相位得到理想对齐相位。量化到 2-bit 后构造对角矩阵计算级联信道增益。运行结果中量化增益会略低于理想增益这个差值就是离散化损失。实际系统中更大的损失来自信道估计误差当 G 和 h_ru 存在 10% 幅度误差和 5° 相位误差时算出的相移已偏离最优值此时继续提高相移分辨率对增益提升有限。联合设计方面文献里最常用的两种做法是两阶段法和交替优化。两阶段法先关闭 RIS 或使用正交导频做信道探测完成 DOA/信道估计后一次性设计权值和相移实现简单但存在误差传播交替优化把问题拆成“固定 Φ 优化 w”和“固定 w 优化 Φ”两个子问题循环求解性能更好但收敛性需要验证。zip 包里若存在 while 循环或 for iter 结构多半是交替优化两阶段法通常只有顺序执行的主脚本没有迭代逻辑。4.3 常见仿真参数设置和选择理由以下表格是 RIS 辅助通信仿真中常见的参数范围便于对照 zip 包里的默认配置判断是否合理参数常见取值说明RIS 单元数 N16 ~ 128越大增益越高但 CSI 开销和优化复杂度随之上涨单元间距(\lambda/2)与天线类似过大会出现反射栅瓣相移量化位 b1 ~ 3 bit1-bit 增益损失明显2-bit 是性能和复杂度平衡点基站天线数 M8 ~ 64决定 DOA/波束成形的空间分辨率链路信噪比-5 ~ 15 dB低于 -5 dB 时 DOA 误差对波束增益的影响显著快拍数 N_snap128 ~ 1024匹配 DOA 高分辨需求N_snap ≥ 4M用户移动性低速/准静态RIS 相移更新周期一般按信道相干时间设计表格看完后要核对 zip 包中的三个尺度关系载频决定波长进而决定阵元间距信道相干时间决定相移更新间隔DOA 估计间隔决定波束更新速率。三者若同处一个数量级系统存在严重的闭环风险——估计完 DOA用户已经移动波束指向已经过时。5. 解开 ZIP 包后用三组指标验证链路真正闭环5.1 先做遮挡实验直射径为零时 DOA 是否依然可用遮挡场景是 RIS 参与链路的最大价值所在。手动把直射增益置零只保留反射径跑一遍 DOA观察谱峰是否还出现在真实入射角附近。如果谱峰消失或偏移超过 5°说明代码包中的 DOA 模块默认依赖直射径存在隐藏的链路假设。这一步比任何调试都值得先做它能暴露多数代码包“跑得通但没闭环”的问题。5.2 三组指标估计误差、波束增益、迭代收敛第一组指标是 DOA 估计误差做 200 次蒙特卡洛仿真统计估计角度的均值和标准差。在快拍数 256、信噪比 10 dB 的配置下8 阵元 ULA 对单目标的理论 CRB 大约在 0.3° ~ 1° 量级实际估计的标准差不应超过该数值的 3 ~ 5 倍。第二组指标是波束成形增益与理想值的差距做法是把 MVDR 权值替换成纯匹配滤波权值比较输出信噪比差异若 MVDR 与匹配滤波增益差小于 1 dB说明协方差估计不准需要增加快拍或对角加载系数。第三组指标是 RIS 迭代收敛曲线输出目标函数值检查算法是否能在 50 次迭代内稳定收敛若振荡幅度不下降排查步长设置和相位初始化方式。5.3 用一张图同时验证三个模块把真实入射角、MUSIC 谱峰位置、MVDR 方向图主瓣指向画在同一张图上。正常闭环时这三个角度的偏差应在 1° 以内出现系统性偏移三个角度输出一致但整体偏移说明角度定义有偏差出现 MUSIC 峰在 0° 附近、MVDR 主瓣在另一侧的情况说明两个模块用了不同的角度坐标系假设。这时对比代码中 angle_grid 的生成方式和 steering_vector 的符号约定通常能快速定位问题。对代码包而言靠人工阅读难以发现的坐标系错位用一张图往往能立即看到——这也是我在解开任何波束成形-DOA-RIS 代码包后必做的第一件事。本文还有配套的精品资源点击获取