简介基于Matlab的多种MPA多用户检测算法仿真资源面向通信工程、电子信息等相关专业本科生、研究生及科研人员重点解决稀疏码多址SCMA等场景下多用户检测算法的实现、对比与性能评估问题。资源共38个文件以36个Matlab脚本.m为主辅以1个说明文档.txt和1个结果示意图.png压缩包仅56KB短小精悍便于直接下载和运行。代码覆盖传统消息传递算法MPA及其改进版本包括Max-log MPA、MSMPA、Threshold MPA、SUS MPA等并按不同用户数、星座点配置、排序与非排序场景预置了多组仿真入口同时提供编码器、译码器、软信息计算、收敛性分析及复杂度计算等辅助脚本模块化程度高便于理解算法原理和修改参数复现结果。已有105人学习适合作为通信方向本科毕业设计、研究生课程项目及科研入门对比调参的参考资料。1. MPA 多用户检测的复杂度瓶颈与 Matlab 实现切入点在 SCMA 非正交多址系统的上行链路里多个用户共享同一组时频资源接收端拿到的每个符号都是多个用户星座点经过信道后的叠加结果。如果直接做最大后验MAP检测需要遍历所有用户星座点的笛卡尔积复杂度随用户数和调制阶数指数上涨J6、M4 时组合数就已经到 4096J 拉到 16 后普通 PC 上仿真基本停摆。MPAMessage Passing Algorithm利用资源块与用户之间的稀疏连接把联合检测拆成变量节点与功能节点之间反复交换外信息复杂度从指数级降到与迭代轮数和局部连接数相关的多项式级。这份 Matlab 资源把标准 MPA、Max-Log MPA、Threshold MPA、SUS 串行 MPA 及 sorted/unsorted 版本都拆成了独立脚本配合 sinfoalt、trellis、decoder、jacobian 等辅助函数在 MATLAB 2014/2019a 下可以直接跑 BER、收敛性和复杂度三类验证。适合通信系统方向的高年级本科生、刚接触非正交多址的研究生也适合想把 SCMA 检测器作为算法模块评估的工程师快速做横向对比。2. 从因子图到 MPA4消息传递的建模与 Matlab 实现先把 File 清单里的 MPA4.m 单独拿出来看因为它是最接近教材的标准实现。MPA 的根基是因子图每个用户对应一个变量节点每个资源块对应一个功能节点用户 j 在资源块 k 上有非零扩频系数就在两个节点之间连一条边。由于 SCMA 码本的稀疏性每个资源块上实际参与叠加的用户只有两到三个所以整张图可以被拆成若干很小的局部更新子图这决定了算法不需要做全排列遍历。2.1 变量节点、功能节点与消息数组设计在 MPA4.m 的多用户检测实现里消息被组织成两个三维数组变量节点消息 V 的维度是(用户数 J, 星座点数 M, 资源块数 K)第三个下标表示这条消息是发给哪个资源块的功能节点消息 F 的维度是(资源块数 K, 星座点数 M, 用户数 J)。这样设计最直接的好处是更新某条边时只需要按 H 矩阵中的非零位置索引邻居不需要维护显式的图对象。初始化时所有星座点的先验概率相等所以 V 和 F 都可以填 1/M。真正有信息量的东西从功能节点更新开始对资源块 k遍历它连接的用户集合固定其中某个用户 u 的星座点枚举其余相关用户所有可能的星座组合计算接收信号 y_k 与“假设组合重构出来的叠加星座点”之间的高斯似然作为功能节点发给用户 u 的消息。随后变量节点将来自多个资源块的消息逐点相乘再归一化完成一轮迭代。MPA4.m 文件本身并没有把整个算法包成单个函数而是把组合枚举、因子更新、硬判决拆到多个局部函数里尤其是 trellis.m 和 decoder.m 承担了星座状态转移和最终软比特解码。我一般会把核心逻辑压缩成下面这个教学版本便于在命令行里单步调function [L_e] simp_mpa4(y, H, constellation, sigma2, Niter) % y : Kx1 接收信号向量 % H : KxJ 等效信道矩阵, 非零元素表示用户 j 占用资源块 k % constellation : Mx1 归一化星座点 % sigma2 : 高斯白噪声方差 % Niter : 消息传递迭代次数 [K, J] size(H); M length(constellation); V ones(J, M, K) / M; % 变量节点消息, 初始等概率 F ones(K, M, J) / M; % 功能节点消息 for it 1:Niter % 功能节点更新: 枚举资源块 k 上的活跃用户 for k 1:K users find(H(k, :) ~ 0); for u users others users(users ~ u); if isempty(others) F(k, :, u) V(u, :, k); continue; end for mi 1:M acc 0; % 构造 others 用户的星座索引组合 grids cell(1, length(others)); [grids{:}] ndgrid(1:M); combos cellfun((x) x(:), grids, UniformOutput, false); combos cell2mat(combos); % 每行是一种干扰组合 for c 1:size(combos, 1) sym constellation(mi) * H(k, u); for t 1:length(others) sym sym constellation(combos(c, t)) * H(k, others(t)); end prob exp(-abs(y(k) - sym)^2 / sigma2); for t 1:length(others) prob prob * V(others(t), combos(c, t), k); end acc acc prob; end F(k, mi, u) acc; end F(k, :, u) F(k, :, u) / sum(F(k, :, u)); end end % 变量节点更新: 相乘再归一化 for j 1:J res find(H(:, j) ~ 0); for k res tmp ones(1, M); for q res(res ~ k) tmp tmp .* F(q, :, j); end V(j, :, k) tmp / sum(tmp); end end end % 计算每个用户的软比特对数似然比 L_e zeros(J, log2(M)); for j 1:J post prod(F(find(H(:, j) ~ 0), :, j), 1); post post / sum(post); for b 1:log2(M) bits de2bi(0:M-1, log2(M)); p1 sum(post(bits(:, b) 1)); L_e(j, b) log(p1 / (1 - p1 eps)); end end end这段代码的每个参数都有明确的物理含义H是等效信道矩阵包含了扩频码和信道增益而不是原始 SCMA 码本sigma2必须和y的量纲一致如果发射端做了归一化噪声功率要按同一基准折算Niter在串行 SUS 算法里通常可以取小一些并行 MPA 则建议 5 到 8 轮。ndgrid枚举只适合用户数较少的场景真正跑 16 用户时组合数会迅速膨胀所以代码包里出现了专门针对 16 用户优化过的MSMPA16_6s1u_sorted.m它不再做全枚举而是用 Max-log 近似把组合搜索截断到高概率子集上。2.2 MPA 迭代中的归一化与对数域切换标准 MPA 每轮迭代结束后都需要把消息归一化否则连乘之后数值会无限趋近 0。变量节点更新时我习惯在prod之后马上除以sum(tmp)这一步不是可选项。另一个容易被忽略的问题是exp(-abs(y(k) - sym)^2 / sigma2)在 SNR 较高时会出现 underflow两个概率都很小直接相除会得到 NaN。文件清单里的 jacobian.m 就是为对数域准备的修正项后面章节会单独拆开讲。2.3 文件命名与参数规模的对应关系拿到代码包时先不用急着跑 main我建议先把文件名里的数字和字母的含义按下面这张表对齐否则很容易把 MPA4 的参数套到 MPA16 上出现维度不匹配的报错。文件片段含义典型对应配置我实际使用时注意的点MPA4 / MPA16因子图规模或星座配置4 点/16 点星座或 4 用户/16 用户先看 H 和 constellation 长度不要凭名字猜3s2u / 6s1u资源块数 s 和用户数 u 的一种缩写3 资源 2 用户6 资源 1 用户类子图这种命名常出现在内部迭代子函数里sorted / unsorted更新顺序是否按信道增益排序串行 MPA 调度策略sorted 版本通常收敛更快但不保证 BER 更优Threshold阈值提前终止判断相邻迭代消息差阈值取 1e-3 到 1e-4太小省不了时间SUS串行更新调度按置信度优先更新适合 16 用户大因子图迭代轮数可减半文件里的 MSMPA16_3s2u221i_sorted.m 这类长后缀含义基本是把 Max-Log MPA 应用到一个大的 16 用户场景再用 sorted 排序策略降低迭代轮数。实际改参数时H 的行数等于资源块数列数等于用户数constellation 的长度等于调制阶数。三者之间没有必然数值相等关系唯一必须保证的是H中非零项所在的列必须对应到用户编号且constellation的索引必须能映射回发射比特。3. 变体实现Max-log、阈值剪枝与串行调度在代码中的落地标准 MPA 可以跑通但工程中没人直接拿它做 16 用户实时检测。这个资源包里真正值钱的是变体实现Max-Log MPA 把指数求和改成最大值近似Threshold MPA 在迭代过程中提前终止SUS MPA 改变更新顺序。下面按代码文件逐一拆开看。3.1 Max-Log MPA 与 jacobian 修正项标准 MPA 的消息更新里反复出现log(sum(exp(...)))这种计算在硬件实现中开销很大而且动态范围容易溢出。Max-Log 近似把所有对数域的指数求和直接替换成取最大值算法变成纯比较与加法代价是外信息被高估BER 通常有 0.1 到 0.3 dB 的损失。代码包里 Max_logMPA4.m、MSMPA16_3s2u_sorted.m 都属于这一类。如果不想牺牲精度又希望保留对数域计算的稳定性可以用 jacobian.m 里的修正项。function lse jacobian(a, b) % 计算 log(e^a e^b) 的稳定形式 % 输入 a, b 是对数域消息 % 输出 lse 是它们的对数求和 m max(a, b); lse m log1p(exp(-abs(a - b))); endlog1p在参数接近 0 时比log(1x)更稳定。写代码时我会把jacobian当作一个累积器反复调用因为多个消息求和等价于两两合并。需要特别留意的是当abs(a-b)大于 30 时exp(-30)已经逼近双精度下限的邻域此时直接返回m即可没必要再做额外运算我在调试 MSMPA16 时会把这一步写成if abs(a-b) 30, lse m; return; end能省掉相当一部分浮点计算。3.2 Threshold MPA 的提前终止与迭代管理ThresholdMPA4.m 和 ThresholdMPA16.m 的核心思想很直接相邻两轮迭代的变量节点消息如果几乎没有变化再迭代下去也只是重复计算。具体做法是在每轮更新前保存V_old更新后计算最大绝对差值低于阈值就跳出循环。V_old V; % ... 执行一轮功能节点和变量节点更新 ... delta max(abs(V(:) - V_old(:))); if delta threshold break; end阈值这个参数对结果影响很大。取 1e-2 可能在第 2 轮就跳出BER 出现平台取 1e-6 基本退化成固定迭代次数失去加速意义。我在 16 用户场景下的经验值是 5e-4配合 Max-Log 变体可以把平均迭代轮数从 8 轮压到 4 到 5 轮。另外要注意阈值判断的基准是消息概率还是对数消息概率域里两个 0.01 的差和对数域里同样的差对应的决策显著度完全不同代码里ThresholdSUSMPA16.m是在 SUS 调度基础上叠加阈值判断基准建议跟主算法保持一致。3.3 sorted/unsorted 与 SUS 串行调度的收敛性差异并行 MPA 的每轮迭代里所有功能节点同时更新消息齐步走SUS 串行 MPA 则按某种顺序逐个更新后更新的消息立刻被后续计算使用信息在单轮内传播得更快。sorted 版本会先更新信道增益高、消息置信度大的节点因此收敛更快但也引入了调度顺序对结果的影响。unsorted 版本顺序固定结果可复现适合做性能基线。文件中的 Sus_MPA4.m、Sus_MPA16.m 和 MSMPASUS16_3s2u221i_sorted.m 把两种思路混在一起先排序再串行更新最后用阈值终止。下表列出我在对比中常用的配置变体对应文件更新策略收敛轮数参考适用场景标准 MPAMPA4.m, MPA16.m并行更新6 到 8性能基准小规模验证Max-LogMax_logMPA4.m, MSMPA16_3s2u_sorted.m并行更新最大值近似6 到 8降低计算量硬件友好ThresholdThresholdMPA4.m, ThresholdMPA16.m并行更新提前终止3 到 5中等用户数加速SUSSus_MPA4.m, Sus_MPA16.m串行更新3 到 516 用户大因子图ThresholdSUSThresholdSUSMPA16.m, MSMPASUS16_3s2u221i_sorted.m串行更新阈值终止2 到 4追求低时延的场景表格里的收敛轮数是在信噪比 0 到 10 dB、M4 星座下统计的如果调制阶数升到 16轮数会继续下降因为每个用户的星座点概率分布更尖锐消息收敛更快但单轮复杂度上升。做横向对比时我给的建议是固定最大迭代次数为 10不要单独用阈值终止这样各组之间只差算法本身不会把调度策略的偏差混进去。4. 主脚本与验证从 main4.m 到复杂度计算文件包里除了算法文件还有一批带 CopyforConvergence 和 CopyforComplexityCal 后缀的脚本。这些脚本不是重复代码而是把同一套链路改造成专门做收敛性观测和复杂度统计的版本。先看懂 main4.m 和 main16.m 的数据流再看这两个副本。4.1 发射端、信道与接收端辅助函数sinfoalt.m负责生成符合 SCMA 帧结构的用户信息比特convencoder.m是卷积编码器配合 trellis.m 做网格描述btod.m和dtob.m做比特与十进制整数之间的映射作用是把星座点索引翻译回比特位decoder.m 和 decoder_log.m 分别对应硬判决和软判决的 Viterbi 解码。我建议第一次运行时不要直接改 SNR 范围先用默认参数跑通 main4.m再逐步替换算法文件。4.2 运行主脚本与切换算法在 MATLAB 里把解压目录加入路径后可以直接运行 main4.m。这个脚本的骨架可以简化成下面这样% main4.m 的核心流程示意 N 4; % 用户数 M 4; % 星座点数 SNRdB 0:2:10; % 信噪比扫描范围 nMC 500; % 蒙特卡洛仿真次数 algo ThresholdMPA4; % 切换: MPA4 / Max_logMPA4 / ThresholdMPA4 / Sus_MPA4 ber zeros(size(SNRdB)); for idx 1:length(SNRdB) berSum 0; for trial 1:nMC [txBits, txSym] sinfoalt(N, M); % ... 叠加信道与噪声 ... [rxBits, ~] run_detector(algo, y, H, constellation, sigma2); berSum berSum sum(txBits ~ rxBits); end ber(idx) berSum / (nMC * N * log2(M)); end semilogy(SNRdB, ber, o-); grid on;这里的run_detector是伪代码位置实际文件里直接调用对应的MPA4.m、Max_logMPA4.m、ThresholdMPA4.m或Sus_MPA4.m。切换算法时只需要保证函数签名一致比如输入都是y, H, constellation, sigma2, Niter输出都是软比特或硬比特。这样改是为了把 BER 对比的变量控制在“算法”这一项上避免由于帧长或信道矩阵不一致带来的额外差异。4.3 收敛性与复杂度验证副本的使用main16CopyforConvergence.m 会把每一轮迭代后的消息差异记录下来画成随迭代次数衰减的曲线main16CopyforComplexityCal.m 则用 tic/toc 统计每个信噪比点、每个信息比特平均消耗的仿真时间。实际测量时我会这样做for idx 1:length(SNRdB) runStart tic; [ber(idx), iterAvg(idx)] sim_mpa16_block(SNRdB(idx), algo, nMC); elapsed toc(runStart); complexityPerBit(idx) elapsed / (nMC * N * log2(M)); end其中sim_mpa16_block是 main16 中的内层循环函数iterAvg是平均迭代轮数。复杂度统计建议每比特时间而不是每轮时间因为 Threshold 和 SUS 变体的平均轮数不同只看总时间无法区分“每轮快”还是“轮数少”。下表是一个我在 16 用户、M4、SNR6 dB 时得到的典型结果算法平均迭代轮数每比特耗时相对值BER 相对值MPA167.81.00基准Max-log MPA167.50.62略高 0.1 dBThreshold MPA164.20.55与 Max-log 接近SUS MPA163.10.41高 0.2 dB 左右5. 多用户检测仿真的边界从 MPA4 迁移到 MPA16 的适配技巧最后落在迁移这件事上。直接把 MPA4.m 里的 J/K/M 改成 16 通常跑不动原因不在算法本身而在枚举组合的叠加态。标准 MPA 在功能节点更新里要枚举其他用户所有星座组合J16 时即使每个资源只连两个用户局部组合数也可能超过可接受范围。我迁移时最先改的不是迭代轮数而是把功能节点更新从全枚举换成 Max-Log 截断只保留每个星座点下概率最大的前若干个干扰组合其余直接忽略。这个思路对应到文件里就是 MSMPA16_6s1u_sorted.m 这类实现。第二个适配点是消息初始化。标准 MPA 初始化为 1/M 没有问题但到了 16 用户高信噪比场景概率域很容易出现极值建议直接把消息初始化为log(1/M)并在整个更新过程中保持对数域。这样在迭代后期不会出现概率乘到 0 导致 NaN 的情况也方便直接输出 LLR。代码里已经准备了 decoder_log.m 和 sinfoalt_log.m就是为对数域链路准备的。第三个技巧是加一个配置自检函数把 H 矩阵、星座点数和消息数组维度一次性检查完。function check_mpa_config(H, M, K) % 检查因子图矩阵 H 和星座点数 M 是否与 K 相容 assert(size(H, 1) K, 资源块数不等于 H 的行数); assert(size(H, 2) max(max(H(:)), 1), 用户索引越界); % 再检查每个资源块上的连接数避免全枚举爆炸 conn sum(H ~ 0, 2); if any(conn 5) M 16 warning(检测到高连接数且 M16建议改用 Max-log 或 SUS 变体); end end调用这个函数不需要额外依赖能省掉大部分维度不对齐的排错时间。最后一个值得注意的细节是 SNR 与噪声方差的换算。通信仿真里通常用每比特能量与噪声谱密度之比作为横轴而 MPA 更新公式里出现的 sigma2 是高斯白噪声方差如果代码里星座点没有归一化sigma2 需要乘上发射功率再参与计算。文件里 main.m 和 main16.m 各自保留了独立的噪声生成方式我在对比时会在开头把两个脚本的EbN0转snr公式对齐否则 MPA4 和 MPA16 的 BER 曲线会整体平移看起来像算法差异实际只是噪声基准不一致。迁移规模时不改噪声基准不省略配置自检然后才谈得上比较 Max-log、Threshold 和 SUS 谁更适合你的目标平台。本文还有配套的精品资源点击获取