简介面向5G非正交多址技术研究者的MATLAB实现包聚焦SCMA系统下的PM-MPA检测算法。资源基于消息传递与最大后验概率思想提供瑞利信道环境中的完整仿真链路适合通信工程高年级学生、算法工程师及科研人员参考复现。包内共5个文件4个m脚本承担核心功能PM_MPA.m实现迭代消息传递过程simulation.m用于配置系统参数并评估误码性能scmaenc.m完成用户数据到稀疏码字的映射log_sum_exp.m则为概率计算提供数值稳定的对数求和工具另附一个压缩算法包可与软判决MPA思路对照使用。整个资源包仅7KB轻量便携下载后即可快速运行。已有611人学习使用是理解SCMA编码原理、PM-MPA迭代译码流程及非正交多址性能权衡的实用入门资料。1. SCMA 的 PM-MPA 检测器这个 matlab 仿真包把瑞利信道下的多用户链路凑齐了PM-MPAProduct Matrix Message Passing Algorithm积矩阵消息传递算法是 SCMASparse Code Multiple Access稀疏码分多址系统里很常用的一类多用户检测算法。我这次拆解的 matlab 源码包把发射侧的scmaenc.m、检测侧的PM_MPA.m、数值工具log_sum_exp.m和顶层仿真脚本simulation.m串成了一条可以直接跑的完整链路信道模型按瑞利衰落来处理。它适合刚接触 SCMA 仿真、想拿现成编码器和检测器跑 BER 曲线的研究生也适合准备做 MPA 变体算法对比的工程师。下文按“文件怎么拆 → 编码侧怎么实现 → 检测侧怎么迭代 → 哪些参数坑最值得注意”的顺序把这份资源讲透。2. 拆包初看从 simulation.m 到 PM_MPA.m 的调用链2.1 五个文件各自的角色压缩包里出现频率最高的五个文件职责边界其实很清晰我在下面按“谁调谁”的顺序列出来文件角色关键观察点simulation.m顶层仿真脚本负责定义 SNR 扫描点、用户数、迭代次数、信道类型瑞利信道建模方式、BER 统计方式scmaenc.mSCMA 编码器把用户比特映射成稀疏码字码本维度、非零元素位置PM_MPA.mPM-MPA 检测算法核心实现消息初始化、因子节点更新、收敛判断log_sum_exp.m对数域求和工具是否做最大值提取能否防下溢scma-SD-MPA.zip另一个 SCMA 检测算法压缩包做软判决对比用可与 PM-MPA 对照性能simulation.m是入口它生成用户比特后调用scmaenc.m经过瑞利信道加噪后交给PM_MPA.m做检测检测结果再和原始比特比对统计误码率。log_sum_exp.m不单独运行它被PM_MPA.m内部循环调用。scma-SD-MPA.zip则是一个独立对比版本解压后可以参照同样的调用方式替换测试。2.2 simulation.m 的主循环结构我把这类 SCMA 仿真工程里最常见的顶层循环抽出来结构基本如下%% simulation.m 主循环结构与常见 SCMA 工程一致 clear; clc; % ------- 基础参数 ------- J 6; % 用户数 K 4; % 资源块数子载波数 M 4; % 码本星座点数量 maxIter 6; % PM-MPA 最大迭代次数 EbN0dB 0:2:12; % 每比特信噪比扫描 trialNum 1e4; % 每个 SNR 点的蒙特卡洛帧数 % ------- 瑞利信道 ------- % 典型做法每个资源块上信道系数独立 % h (randn 1i*randn) / sqrt(2) % ------- 主循环 ------- for ebnoIdx 1:length(EbN0dB) errorCount 0; bitCount 0; for trial 1:trialNum % 1) 生成随机比特调用 scmaenc 得到 K 维发送向量 % 2) 乘上瑞利衰落系数叠加复高斯白噪声 % 3) 调用 PM_MPA.m 做检测得到估计比特 % 4) 与原始比特比对累计 errorCount / bitCount end ber(ebnoIdx) errorCount / bitCount; end % ------- 画图 ------- semilogy(EbN0dB, ber, -o); grid on;逻辑说明这里用EbN0dB而不是SNR是通信仿真里的习惯因为不同调制阶数下每个符号携带的比特数不同只有折算到每比特能量多条 BER 曲线才有可比性。循环内每帧数据都走一遍“编码 → 信道 → 检测 → 比对”四个步骤最后统计误码数除以总比特数得到 BER。参数说明J6表示 6 个用户共用 4 个资源块过载率 150%这是 SCMA 参考设计中常见的配置maxIter不建议一开始设太大先设 6 跑通流程再逐步加大观察性能变化trialNum在低误码率区域要适当增大否则曲线尾部会抖动得很厉害。2.3 拿到代码后先做三项自查第一次运行这份资源先别急着改参数我一般会做三个快速检查。第一是看当前目录是否已经addpath到所有.m文件所在路径MATLAB 报“未定义函数或变量”八成是路径问题。第二是在simulation.m里搜索码本定义位置这类代码的码本矩阵有时直接写在主脚本里有时单独放在一个codebook.m里务必保证它能被scmaenc.m和PM_MPA.m同时访问。第三是检查log_sum_exp.m文件尾是否有多余测试代码很多下载版代码文件尾部残留调试输出会影响仿真效率。做完这三项基本可以确认代码环境是可复现的状态。3. 编码侧实现scmaenc.m 如何把比特变成稀疏码字3.1 从比特到稀疏码字的两次映射SCMA 的编码过程和传统 CDMA 最大的区别在于“稀疏”两个字。传统 CDMA 每个用户都会扩展占用全部资源SCMA 则让每个用户只占用其中一部分资源留下大量零元素。scmaenc.m做的事情本质上是两次映射第一次把比特组合映射成星座点索引第二次把星座点索引映射成 K 维码字。% scmaenc.m 的核心思路非原包逐行照贴变量名按习惯改写 % codebook: K x M x J 三维矩阵 % dataBits: J x log2(M) 的逻辑比特矩阵 function tx scmaenc(dataBits, codebook) [J, bitsPerSym] size(dataBits); [K, M, ~] size(codebook); tx zeros(K, 1); % 发射向量初始化为 0 for j 1:J % 第一次映射比特 - 符号索引1~M symIdx bi2de(dataBits(j, :), left-msb) 1; % 第二次映射符号索引 - 稀疏码字并叠加到资源上 tx tx codebook(:, symIdx, j); end end逻辑说明codebook(:, symIdx, j)取出第 j 个用户在第symIdx个星座点上的 K 维码字。由于 SCMA 码字是稀疏的这个向量里大部分位置是 0只有少数几个位置非零。循环内做的是多用户信号在同一组资源上的叠加这也是 SCMA “非正交”的直接体现。参数说明bi2de(..., left-msb)的左右 MSB 设置会影响索引顺序如果发射端和接收端用的映射规则不一致BER 曲线会直接崩溃。码本矩阵的第三维是用户序号第二维是星座点序号第一维是资源序号这个维度约定在PM_MPA.m里同样适用改动时两边必须同步。3.2 码本结构与稀疏因子图参数SCMA 的性能很大程度上取决于码本设计。我常见到的资源包里码本采用 6 用户 4 资源的结构也就是 J6、K4、M4每个用户只占据 2 个资源块。参数典型值含义J6用户数K4资源块数M4每个用户的星座点数量df2每个用户非零资源块数过载率150%J/K表示频谱资源复用程度单资源重叠用户数3每个资源块上叠加的用户数这组参数的含义是6 个用户的数据挤在 4 个资源块上发射每个资源块上的信号由 3 个用户的信号叠加而成。df2意味着每个用户只在 2 个资源块上放置能量其余 2 个资源块上为空。这种稀疏性让接收端可以用因子图描述用户和资源之间的关系从而用消息传递算法以较低复杂度完成多用户分离。3.3 验证编码器是否正确的两个自检点下载资源最容易出的问题就是码本文件被误改或复制错位我每次拿到新代码都会跑一下这个自检脚本% 检查第 j 个用户的码字稀疏度是否等于 df j 3; codewords reshape(codebook(:, :, j), K, M); nzCount sum(abs(codewords) 1e-12, 1); % 统计每列非零数 disp(unique(nzCount)); % 期望输出: 2 % 检查每个码字能量是否归一化 energy sum(abs(codewords).^2, 1); disp(energy); % 期望接近 1逻辑说明第一个检查保证码本的稀疏结构与PM_MPA.m里因子图矩阵的假设一致。如果unique(nzCount)输出的是 3 而不是 2说明码本数据错位后续检测算法会把不存在的连接关系当成有连接导致消息更新混乱。参数说明1e-12是判断是否为 0 的阈值因为浮点运算中真正的 0 可能被存成极小的残留值。能量归一化检查则确保每个码字等概率发射避免某个星座点功率异常偏高影响 BER 结果的真实性。4. 检测侧核心PM_MPA.m 里的消息迭代、乘积矩阵与 log_sum_exp4.1 从 MAP 到 MPA 再到 PM-MPA三次复杂度取舍接收端要做的事是从叠加了多用户信号和噪声的 K 维向量里恢复每个用户的比特。最理想的是 MAP 检测但对 6 用户 4 星座的配置一次联合遍历就是M^J 4096种组合调制阶数再高就完全跑不动。MPA 的做法是在因子图上做消息传递每个资源块只需遍历M^df种组合。检测方案单资源遍历量复杂度特征联合 MAP/ML4096指数爆炸不可实际使用标准 MPA16对每个用户状态遍历邻居码字组合PM-MPA约 M²用乘积矩阵缓存减少重复计算PM-MPA 的改进思路在于当某个资源块上重叠了多个用户时标准 MPA 对每个用户都要重新计算一遍邻居用户的联合概率PM-MPA 把这些重复计算整理成一份乘积矩阵缓存更新单个用户时通过“整体乘积”剔除自己那一路。在PM_MPA.m里这个技巧体现为因子节点更新时先算临时累乘再逐用户取值。4.2 PM_MPA.m 的迭代骨架function [bitsEst, llrOut] PM_MPA(y, H, codebook, N0, maxIter, convTh) % y: Kx1 接收向量 % H: KxJ 瑞利信道系数矩阵 % codebook: KxMxJ 码本 % N0: 噪声单边功率谱密度 % maxIter: 最大迭代次数 % convTh: 收敛门限 [K, M, J] size(codebook); % 变量节点消息初始化为等概率 Mv2f ones(M, J) / M; for iter 1:maxIter % ---- 因子节点更新 ---- Mf2v ones(M, J); for k 1:K % 找到占用第 k 个资源的用户集合 userIdx find(abs(codebook(k, 1, :)) 1e-12); % 先算该资源上所有用户消息的乘积矩阵PM 核心 prodMsg ones(M, 1); for jj 1:length(userIdx) prodMsg prodMsg .* Mv2f(:, userIdx(jj)); end % 对每个用户单独生成因子节点消息 for jj 1:length(userIdx) j userIdx(jj); % 剔除自己后与信道、噪声相关的指数项相乘 % 常见实现里会调用 log_sum_exp 做累加 Mf2v(:, j) prodMsg ./ Mv2f(:, j) .* exp(-abs(y(k) - H(k,j) * codebook(k,:,j).).^2 / N0); Mf2v(:, j) Mf2v(:, j) / sum(Mf2v(:, j)); end end % ---- 变量节点更新 ---- Mv2f ones(M, J); for j 1:J resIdx find(abs(codebook(:, 1, j)) 1e-12); % 将自己所占资源上的因子节点消息做乘积 for kk 1:length(resIdx) Mv2f(:, j) Mv2f(:, j) .* Mf2v(:, resIdx(kk), j); end Mv2f(:, j) Mv2f(:, j) / sum(Mv2f(:, j)); % 归一化 end % ---- 收敛判断 ---- if iter 1 max(abs(Mv2f - Mv2f_old), [], all) convTh break; end Mv2f_old Mv2f; end % 硬判决取每列最大概率对应的符号 [~, symIdx] max(Mv2f, [], 1); bitsEst de2bi(symIdx - 1, log2(M), left-msb); llrOut Mv2f; end逻辑说明因子节点更新那段里prodMsg就是 PM-MPA 的乘积矩阵缓存。它先把某个资源上所有用户的消息乘在一起然后对某个用户更新时直接拿总乘积除掉自己省去了逐个重新遍历邻居状态的开销。这个“先整体乘、再逐个除”的操作正是 PM-MPA 相比标准 MPA 的核心区别。参数说明convTh是收敛门限典型值在1e-3到1e-5之间。设太大迭代提前终止BER 变差设太小迭代次数拉满复杂度优势消失。N0必须和simulation.m里的噪声功率保持一致否则检测器内部指数项的权重是错的高 SNR 区域的表现会非常奇怪。4.3 log_sum_exp 在消息更新中的用法MPA 类算法的消息更新里经常出现形如log(exp(a) exp(b))的运算而 matlab 原生log(sum(exp(x)))在 x 取较大负值时exp会直接下溢成 0log 再取就变成-Inf。资源包里单独放一个log_sum_exp.m就是为了解决这个数值问题。function y log_sum_exp(x, dim) % 对数域求和计算 log(sum(exp(x))) 的数值稳定版本 if nargin 2 dim 1; end maxVal max(x, [], dim); y maxVal log(sum(exp(x - maxVal), dim)); end逻辑说明核心技巧是先减去最大值再求指数把指数函数的自变量整体平移到非正区间。这样exp(x - maxVal)的最大值是 1不会上溢最小值受精度限制也不会轻易下溢成 0。参数说明dim指定求和维度默认是 1。在实际使用中如果消息矩阵是M x J维想按用户维度求和就把dim设为 2。需要注意的是这个函数只接受实数输入复数相加必须先拆成实部虚部分别处理。5. 避坑排查瑞利信道归一化、log_sum_exp 溢出与迭代不收敛5.1 现象BER 曲线比文献差好几个 dB平躺下不去原因最常见的是simulation.m里瑞利信道系数没有归一化。(randn 1i*randn)产生的信道功率是 1但很多人会忘记除sqrt(2)导致信道增益偏大等效噪声被低估。另一个原因是把EbN0和SNR混用SCMA 多用户叠加后每个资源上的符号能量不等于单个用户的比特能量。解决检查信道生成代码确认h (randn 1i*randn) / sqrt(2)再看N0的计算是否用10^(-EbN0dB/10)并除以每比特对应资源数。我一般会在simulation.m里加一行disp(norm(h, fro)^2 / K)理想值接近 1。5.2 现象log_sum_exp 输出 NaN 或 InfBER 曲线出现断崖原因log_sum_exp输入里出现NaN通常是消息矩阵中出现了零概率值归一化时除数为 0。另一个可能是exp(x - maxVal)里x - maxVal全部为非常大负数sum 之后为 0取 log 得到-Inf而不是有效数值。解决在PM_MPA.m的变量节点更新后加一个保护Mv2f max(Mv2f, eps);然后再归一化。同时检查收敛判断时是否用了abs(Mv2f - Mv2f_old)如果消息矩阵本来就是概率值diff量级很小直接用maxDiff convTh即可不需要再取对数。5.3 现象高信噪比区域误码率下降变缓出现错误平台原因PM-MPA 是近似算法迭代次数固定为 6 时高 SNR 区域残留误差主要来自消息近似而非噪声。另外一个容易被忽视的点是收敛门限convTh设得太大算法在还没收敛时就提前退出。解决把maxIter从 6 增大到 10convTh从1e-3收紧到1e-5观察曲线尾段是否改善。如果平台还在检查信道补偿逻辑H(k,j)在深衰落位置幅度接近 0直接除会放大噪声常规做法是在补偿时加一个小常数delta1e-6保护分母。5.4 现象PM-MPA 和标准 MPA 性能完全重合复杂度优势体现不出来原因代码里没有真正缓存乘积矩阵只是在循环内重复计算本质还是标准 MPA。很多版本的PM_MPA.m只是把标准 MPA 的变量名改成 PM内部逻辑没有体现“先乘整体、再除自己”的操作。解决检查因子节点更新部分是否在for jj循环外先算了prodMsg。如果没有参考本文 4.2 节的写法把乘积矩阵提到内层循环外面。验证方式很简单在PM_MPA.m里记录每次迭代的乘法次数和标准 MPA 版本对比差距不明显就说明缓存逻辑没生效。6. 进阶把 PM-MPA 和 SD-MPA 画到同一张 BER 图上做对比6.1 复用一份仿真脚本的快速做法拿到scma-SD-MPA.zip后不需要另写仿真框架。把simulation.m里调用的检测函数封装一层用函数句柄切换算法是最省事的做法% 统一检测接口 detectMPA (y, H, codebook, N0) PM_MPA(y, H, codebook, N0, 6, 1e-4); detectSDMPA (y, H, codebook, N0) SD_MPA(y, H, codebook, N0, 6, 1e-4); % 在 SNR 循环里调用其余代码完全复用 [bitsEst, ~] detectFunc(y, H, codebook, N0);对比时建议固定同一组随机种子确保两种算法经受完全相同的信道和噪声样本。这样画出来的 BER 曲线差异才能只反映算法本身的能力。6.2 判读对比结果时看什么性能上SD-MPA软判决通常比 PM-MPA 有零点几个 dB 的增益尤其在低迭代次数下。复杂度上PM-MPA 在中低 SNR 区域迭代 3 到 4 次就能收敛而标准 MPA 往往需要 6 次以上。如果两条曲线完全重合要检查迭代次数是否被固定成相同值PM-MPA 的优势场景是“相同迭代次数下更快收敛”。我习惯在代码里加tic/toc记录每个 SNR 点的平均耗时这条时间曲线往往比 BER 曲线更能说明 PM-MPA 的价值。从那以后我每次换信道模型或调制阶数都强制走一遍“信道归一化检查、消息保护、迭代收敛性观察”这三个步骤短则五分钟长则半小时能替后面省下大把调参时间。希望这份拆解能让你少踩几个我踩过的坑。本文还有配套的精品资源点击获取
