简介面向IRS辅助NOMA上行链路研究者的论文复现代码包围绕Sum Rate Maximization for IRS-Assisted Uplink NOMA图2与图3结果进行仿真适合通信工程、信息与信号处理方向的硕博生、科研人员用于算法验证、曲线对比和课堂研讨。压缩包内共12个文件其中8个Matlab脚本涵盖信道生成、IRS相位设计、NOMA功率分配、SDR优化求解及OMA基准对比等关键环节4个fig图形文件可直接查看迭代收敛与性能变化趋势。整包大小约129KB轻量易用下载后配合主脚本运行即可快速产出论文对应结果。目前已有584人学习下载属于典型的小型复现类学习资源。通过运行两个主脚本可分别得到图2与图3帮助读者从代码层面理解IRS辅助NOMA和速率最大化的建模思路与实现细节节省自行推导与编程调试的时间。1. 这套IRS辅助上行NOMA链路代码到底在仿真什么仿真跑到第三天我一度以为是 SIC 解码的索引写错了——同一套 IRS 辅助上行 NOMA 链路代码换个相移初始化方式两条用户速率曲线就交叉得完全没法解释。最后定位到的问题不在解码顺序而在信道生成时把级联信道的方向搞反了。这类仿真表面上是把 IRS 相移和 NOMA 功率分配叠起来算用户速率实际上几乎所有的坑都藏在信道维度、线性/对数单位换算和解码顺序这三处。这套 IRS 辅助上行 NOMA 链路代码就是用来量化智能反射面对上行非正交多址接入系统速率增益的蒙特卡洛仿真框架适合做物理层算法预研、课程设计和毕业设计的开发者。它不碰硬件原型但能在你投入真实验证平台之前先把 IRS 的增益边界、NOMA 功率分配策略和相移量化的影响摸清楚。2. 先理解仿真对象IRS辅助上行NOMA的链路模型与三个核心模块2.1 链路里到底有几条信道直接链路与级联信道的建模差异上行 NOMA 场景里K 个单天线用户同时向基站发射信号基站附近部署一个带有 N 个反射单元的 IRS。每个用户到基站有两条路径一条是直接信道用户→基站另一条是反射信道用户→IRS→基站。反射信道是两跳信道的乘积不能像直接信道那样一次性生成完事代码里必须把「用户到 IRS」和「IRS 到基站」分开建模再合成等效级联信道。维度约定是第一个容易翻车的地方。我一般统一用这套约定直接信道 h_d 存成 1×K 行向量用户到 IRS 的信道 h_ru 存成 K×N 矩阵IRS 到基站的信道 h_br 存成 N×1 列向量相移向量 phi 也是 N×1 列向量。这样等效反射信道可以写成(h_ru .* phi.) * h_br得到 K×1 再转置成 1×K和 h_d 直接相加就是用户 k 的完整等效信道。符号维度含义h_d1×K用户直接到基站的复信道h_ruK×N各用户到 IRS 各反射单元的复信道h_brN×1IRS 到基站的复信道phiN×1IRS 各单元相移单位模约束衰落模型的选择也要区分两条链路。直接链路通常用瑞利衰落就够路径损耗指数按城市宏蜂窝取 3.5 左右级联链路因为 IRS 一般部署在建筑高处视距概率高建议用莱斯衰落莱斯因子取 5~10 dB。代码里如果统一用瑞利会把 IRS 的增益低估不少因为少了视距分量的稳定贡献。路径损耗方面级联链路要算两段损耗的和用户到 IRS 一段IRS 到基站一段。初版代码可以合并成一个等效路径损耗 PL_r这样参数少、跑得快等后面要对比 IRS 部署位置时再把两段拆开。小尺度衰落生成时注意每个复数元素的实部和虚部要独立同分布、方差各为 1/2这样信道模平方的期望才是 1路径损耗才能单独叠加。2.2 无源波束成形与功率分配代码里两个可调旋钮IRS 的本质是一个无源反射阵列它不能放大信号只能对每个反射单元的信号做相位旋转这就是「无源波束成形」。代码里每个单元的相移是一个单位模复数exp(1j*theta)整个 IRS 的相移向量 phi 必须满足逐元素模为 1 的约束。如果代码里出现幅度不为 1 的相移仿真结果就会凭空多出增益或损耗这是物理上不允许的。最常见的启发式相移设计是相位对齐取某个用户到基站的级联信道对每个反射单元的相位取负让所有反射信号在基站处同相叠加。这个做法不是最优的但它是所有对比的基线。真正的最优相移需要联合优化功率分配属于非凸问题仿真初期先跑启发式后面再上优化算法对比否则你很难判断优化带来的增益是真实的还是 SIC 写错了。NOMA 部分的旋钮是用户发射功率。上行 NOMA 和下行 NOMA 有本质区别下行是基站总功率约束基站给不同用户分配功率上行是每个用户各自有功率约束用户之间不能转移功率。代码里把所有用户的发射功率存成一个 1×K 数组初版直接全部设成等功率先把链路跑通再改成功率分配策略。2.3 代码组织方式一个主脚本加三个函数文件这套仿真不要写成一个 800 行的脚本我习惯按功能拆成四个文件放在同一个目录下main_irs_noma_ul.m主脚本负责参数设置、蒙特卡洛循环、结果累计和画图gen_channels.m生成直接信道和级联信道返回 h_d、h_ru、h_brdesign_phase.m设计 IRS 相移向量返回 phicompute_rate_ul_noma.m给定等效信道、功率和噪声按 SIC 顺序计算各用户速率函数化的理由很简单信道生成、相移设计、速率计算这三个模块在调试阶段会被反复调用你改参数、换相移策略、加信道估计误差时只需要动对应函数的内部实现主脚本不用改。每个函数头部我都习惯写三行注释输入维度、输出维度、单位。IRS 辅助 NOMA 这种多维复数运算维度一致性和单位一致性决定了仿真代码的命。3. 从零搭好IRS辅助上行NOMA链路代码参数表、信道生成与相移初始化3.1 仿真参数先定死一张可抄的参数表参数不统一的代码最容易翻车。我建议先把下面这张表当作初始配置写进主脚本开头跑通之后再按需调整。K 先取 2双用户场景的 SIC 索引关系最直观排查问题容易N 取 32 起步既能看到 IRS 增益又不会因为复杂度卡住。参数符号建议取值说明用户数K2先跑双用户SIC 索引不容易错IRS 单元数N32单元数翻倍时理想对齐增益约 3 dB用户发射功率P10 dBm每用户独立代码内用线性功率存带宽B1 MHz决定噪声底噪声系数NF5 dB接收机底噪抬升直接链路损耗PL_d100 dB用户到基站的平均路径损耗级联链路损耗合并PL_r80 dB用户-IRS-基站两段合计莱斯因子Rice10 dB级联信道视距分量占比蒙特卡洛次数MC10000少于 5000 曲线抖动会很严重相移量化q_bitinf连续先连续相移后加 2 bit 量化对比噪声功率的换算这里最容易出错。标准做法是用noise_power 10^((-174 10*log10(B) NF)/10)算出瓦特数再和用户功率一起参与 SINR 计算。更省事的做法是把噪声功率归一化为 1用信噪比直接反推用户线性功率但无论哪种都要在参数区注释里写清楚用的哪种约定不然过两周回来看代码一定会懵。3.2 生成信道瑞利衰落与莱斯信道的代码写法信道生成函数返回三个量直接信道、用户到 IRS 的信道、IRS 到基站的信道。路径损耗在线性域计算小尺度衰落部分要区分瑞利和莱斯两种写法。function [h_d, h_ru, h_br] gen_channels(params) K params.K; N params.N; % 路径损耗 dB - 线性 PL_d 10^(-params.PL_d_dB/10); PL_r 10^(-params.PL_r_dB/10); rice_lin 10^(params.Rice_dB/10); % 莱斯因子转线性 % 直接链路瑞利衰落实部虚部分别是方差1/2的高斯 h_d sqrt(PL_d/2) * (randn(1,K) 1j*randn(1,K)); % 用户-IRS莱斯衰落含视距分量 los_ru exp(1j*2*pi*rand(K,N)); % 简化视距相位 nlos_ru sqrt(PL_r/(2*(rice_lin1))) * (randn(K,N) 1j*randn(K,N)); los_ru sqrt(PL_r*rice_lin/(rice_lin1)) * los_ru; h_ru nlos_ru los_ru; % IRS-基站同样的莱斯结构 los_br exp(1j*2*pi*rand(N,1)); nlos_br sqrt(PL_r/(2*(rice_lin1))) * (randn(N,1) 1j*randn(N,1)); los_br sqrt(PL_r*rice_lin/(rice_lin1)) * los_br; h_br nlos_br los_br; end这段代码的关键在莱斯衰落的功率分配nlos 分量除以(rice_lin1)再开方los 分量按rice_lin/(rice_lin1)缩放这样合成后信道平均功率恰好是PL_r路径损耗和小尺度衰落互不干扰。如果这里写反IRS 的平均增益会被低估一半以上。参数说明randn生成的实部虚部各有方差 1所以除以sqrt(2)让每个分量方差变成 1/2合成信道模平方期望为 1。莱斯因子建议先传 dB 值再在函数内转线性参数区看着更直观。视距分量的相位这里用均匀随机量简化严格仿真应该用阵列响应向量但初版这个精度足够看趋势。3.3 相移初始化先对齐相位再看NOMA增益相移设计函数只做一件事把某个用户的级联信道相位取负让各反射单元的反射信号在基站处同相叠加。function phi design_phase(h_ru, h_br, target_user) % 输入: h_ru (KxN), h_br (Nx1), target_user 对齐目标用户 % 输出: phi (Nx1) 单位模相移向量 cascade h_ru(target_user,:). .* h_br; % Nx1 各单元级联信道 phi exp(-1j * angle(cascade)); % 相位取负实现对齐 % 单位模约束abs(phi) 全为 1 endcascade是逐元素乘积不是矩阵乘法含义是每个反射单元从目标用户到基站各自形成的级联信道。angle取出相位exp(-1j * angle)把它旋转到零度这样每个单元的反射信号到达基站时相位一致、直接相加。参数说明target_user选谁直接决定多用户公平性——对齐强用户弱用户反射增益会很小对齐弱用户弱用户速率提升明显但和速率不一定最高。初版先固定对齐用户 1跟随机相移做对比。随机相移的基准也很简单phi exp(1j*2*pi*rand(N,1))就行做一个函数或直接在循环里生成都可以。4. 把NOMA解码和速率统计写进链路叠加编码、SIC顺序与蒙特卡洛循环4.1 上行NOMA的叠加信号与SIC解码顺序基站接收信号是所有用户发射信号经过各自等效信道后的叠加表达式为 y Σ sqrt(P_k) h_k s_k n。每个用户的等效信道是直接信道加反射信道h_eff(k) h_d(k) sum(h_ru(k,:) .* phi. .* h_br.)代码里用矩阵乘法一次算完。上行 NOMA 的 SIC 解码顺序和下行相反。下行基站给弱用户分配更多功率所以基站先解强用户信号上行每个用户功率独立常见做法是按信道增益升序解码也就是信道最弱的用户最先解。原因是弱用户不先解的话会被强用户信号完全压住第一步判决就出错后面误差传播更严重。解码次序解码对象此时干扰来源第 1 个信道增益最小的用户其余所有用户第 2 个信道增益次小的用户除去已解码用户的其余用户第 K 个信道增益最大的用户无只剩噪声这个顺序不是唯一选择有些方案会按速率公平性动态调整但初版代码用「弱先解」最容易验证。SINR 计算时要注意对当前解码用户干扰只来自排在它后面的未解码用户已经解掉的用户信号要从干扰项里剔除。4.2 速率计算的完整函数一个蒙特卡洛循环的骨架速率计算函数接收等效信道、线性功率、噪声功率和解码顺序输出各用户速率。这个函数把 SIC 的干扰逻辑封装好主脚本只需要调用它。function rate compute_rate_ul_noma(h_eff, Ptx_lin, noise_power, order) K length(h_eff); rate zeros(1,K); for idx 1:K m order(idx); % 当前要解码的用户 inter 0; % 干扰只来自还没解码的用户 for j idx1:K l order(j); inter inter Ptx_lin(l) * abs(h_eff(l))^2; end sinr Ptx_lin(m) * abs(h_eff(m))^2 / (inter noise_power); rate(m) log2(1 sinr); % bps/Hz end end这里用「累计剩余干扰」的写法而不是先算完整信号再逐次相减避免了完美 SIC 下信号重建的数值误差。如果有人把 SIC 写成了每次减去Ptx_lin(m)*h_eff(m)再继续碰到信道条件差、SINR 为负的情况速率会变成 NaN。主脚本的蒙特卡洛骨架如下%% main_irs_noma_ul.m clear; clc; rng(20240618); % 固定种子结果可复现 params.K 2; params.N 32; params.Ptx_dBm 10; params.noise_dBm -109; params.PL_d_dB 100; params.PL_r_dB 80; params.Rice_dB 10; Ptx_lin 10^(params.Ptx_dBm/10) / 1000 * ones(1, params.K); noise_power 10^(params.noise_dBm/10) / 1000; MC 10000; rate_acc zeros(1, params.K); for mc 1:MC [h_d, h_ru, h_br] gen_channels(params); phi design_phase(h_ru, h_br, 1); % 对齐用户1 h_eff h_d ((h_ru .* phi.) * h_br).; % 1xK [~, order] sort(abs(h_eff).^2, ascend); rate compute_rate_ul_noma(h_eff, Ptx_lin, noise_power, order); rate_acc rate_acc rate; end rate_avg rate_acc / MC;(h_ru .* phi.) * h_br这一步做了三件事phi.把 N×1 转成 1×N点乘让每个用户的信道逐单元乘上相移再乘h_br完成反射链路叠加最后转置成 1×K 和h_d对齐。MATLAB 的隐式扩展需要 R2016b 之后的版本如果用的是老版本要手动repmat(phi., K, 1)。参数说明Ptx_lin的单位是瓦特noise_power也是瓦特两者量纲一致才能算 SINR。sort(..., ascend)对应弱用户先解。速率累计除以 MC 就是遍历速率横轴扫发射功率时这个循环结构完全不用改只动参数区。4.3 对比基准怎么选OMA、随机相移、全用户等功率没有基准的仿真曲线没有说服力。最少要跑四个场景放在一张图里对齐相移加 NOMA、随机相移加 NOMA、对齐相移加 OMA、无 IRS 直接链路。每个场景换一下design_phase的返回值和接入方式即可。对比场景设置看什么对齐相移 等功率 NOMAphi 对齐用户 1Ptx 等功率本文主推场景随机相移 NOMAphi 随机生成IRS 无源波束成形的增益对齐相移 OMA每个用户独占总资源NOMA 相对 OMA 的频谱效率增益无 IRS 直接链路去掉反射项h_eff h_dIRS 部署本身的价值OMA 基准有个新手必踩的坑每个用户只能分到 1/K 的资源速率要乘 1/K即rate_oma (1/K) * log2(1 Ptx*abs(h_d)^2 / noise_power)。很多人的 OMA 曲线比 NOMA 还高就是忘了除资源份额。横轴一般扫用户发射功率从 0 dBm 到 20 dBm纵轴画和速率或最弱用户速率。这时候能清楚看到两个问题IRS 对齐带来的增益随功率怎么变化NOMA 相对 OMA 的增益在哪个功率区间被放大或吞掉。5. IRS辅助上行NOMA链路代码避坑5 个高频翻车点与修正方法5.1 翻车点1相移向量维度写错级联信道拼不起来现象运行时报矩阵维度不匹配或者h_eff变成 N×K 的矩阵画出来的速率曲线乱成一团。这事我遇到过不止一次基本都是phi用了行向量、h_br是列向量点乘时隐式广播把维度搞爆炸。原因h_ru是 K×Nh_br是 N×1phi如果是 1×N那么h_ru .* phi得到 K×N 没问题但(h_ru .* phi) * h_br维度是 K×1如果你忘了转置后面和h_d的 1×K 对不上如果phi用了 N×1直接h_ru .* phi就会报错。解决统一约定phi永远存 N×1 列向量h_br也是 N×1。在design_phase函数头部加断言让问题在源头暴露而不是等到主脚本里报错assert(size(phi,1) N size(phi,2) 1, phi必须为Nx1列向量); assert(size(h_br,1) N size(h_br,2) 1, h_br必须为Nx1列向量);5.2 翻车点2功率用dBm直接参与线性叠加SINR变成负数现象计算出来的 SINR 是负值速率全是 NaN曲线根本没法看。第一次跑 NOMA 仿真的人十个有八个栽在这。原因dBm 是对数量纲两个 dBm 值相加相当于线性值相乘不是叠加。把 Ptx_dBm 和 noise_dBm 直接带进Ptx*abs(h)^2 / noise分母分子量纲混乱结果完全失真。解决所有功率在进入公式前统一转成瓦特。dBm 转瓦特是P_watt 10^(P_dBm/10) / 1000只在最后画图、标纵轴刻度时才转回 dBm。主脚本里换算完加一行注释避免下次手滑。5.3 翻车点3SIC把想要的信号一起消掉了现象强用户的速率算出来比弱用户还低甚至低于 OMA 基准曲线完全违背直觉。原因解码顺序写反了从强用户开始解且重建信号时把当前用户也减掉了。更隐蔽的原因是干扰项里包含了已解码用户导致强用户 SINR 被恶意压低。解决算 SINR 用「剩余干扰」写法当前用户的信号功率永不出现在干扰项里。解码顺序用排序索引取不要手写order [weak_user, strong_user]这种硬编码因为信道排序每次蒙特卡洛实现都不一样。% 正确当前用户m不在干扰项里 inter 0; for j idx1:K l order(j); inter inter Ptx_lin(l) * abs(h_eff(l))^2; end5.4 翻车点4蒙特卡洛次数不够相移增益看不出来现象对齐相移的曲线和随机相移交叉甚至低于随机相移看起来像是 IRS 没有用。这种问题最容易让人怀疑「无源波束成形是不是玄学」。原因相位对齐本身就是针对随机信道实现的单次信道里对齐带来的增益方差很大MC 次数少时随机信道没有充分平均曲线抖动超过实际增益差。解决MC 至少 10000 起步曲线不平滑就加到 50000。固定rng种子确保每次跑同参数得到同一条曲线否则你无法判断曲线抖动是算法问题还是随机性。调参阶段先用 N16、MC2000 快速验证趋势确认无误再放大参数。5.5 翻车点5IRS单元数N一大就跑不动复杂度爆炸现象N64 还能忍N256 时一个蒙特卡洛循环跑几分钟内存也不断攀升。原因有人在相移设计里写了for n 1:N逐单元循环或者把每次信道实现都存进一个大矩阵MC 还没跑完内存先爆了。解决信道生成和相移设计全部矢量化不要对 IRS 单元写循环。主循环里只累加速率不存中间信道矩阵内存占用和 MC 次数无关。N 很大时可以把 MC 分块每块 1000 次块结束清空临时变量防止 MATLAB 内存碎片化。6. 让仿真结果更可信验证曲线、基准对比与进阶调试习惯拿到第一版曲线先别急着加复杂算法先做退化验证。把 K 设成 1NOMA 退化成单用户上行链路此时没有 SIC信道只剩h_eff h_d sum(|cascade|)。固定一个信道实现手算Ptx * abs(h_eff)^2 / noise_power对应的速率和代码打印值对比。对不上就逐个检查信道生成和相移函数对得上再去跑多用户。第二个习惯是固定信道快照调试。把蒙特卡洛循环关掉只跑一次gen_channels→design_phase→compute_rate_ul_noma打印每个中间变量的disp。重点核对三处h_eff 的实部虚部量级是否和路径损耗匹配排序后的 order 是不是符合信道强弱直觉SINR 公式里干扰项拉出来手算一遍。这套操作能过滤掉八成以上的逻辑错误。验证通过后进阶方向有三个相移量化、信道估计误差、用户配对。相移量化是把连续phi改成 2 bit 或 1 bit 量化观察 IRS 增益损失这个改动很小但非常值得做信道估计误差是在 h_eff 上叠加一个百分比误差用来评估 IRS 辅助系统的鲁棒性K 大于 2 时还要加用户配对策略不同用户组合对 NOMA 增益影响很大。我自己的习惯是每个函数文件头部固定写三行注释输入维度、输出维度、单位。IRS 辅助 NOMA 这种到处是复数矩阵、功率换算的仿真维度一致性和单位一致性就是代码的生命线。这个习惯帮我少走了很多弯路也希望帮到你。本文还有配套的精品资源点击获取