简介：面向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-MPA（Product Matrix Message Passing Algorithm，积矩阵消息传递算法）是 SCMA（Sparse 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.m | SCMA 编码器，把用户比特映射成稀疏码字 | 码本维度、非零元素位置 |
PM_MPA.m | PM-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。
参数说明：J=6表示 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 资源的结构，也就是 J=6、K=4、M=4，每个用户只占据 2 个资源块。
| 参数 | 典型值 | 含义 |
|---|---|---|
| J | 6 | 用户数 |
| K | 4 | 资源块数 |
| M | 4 | 每个用户的星座点数量 |
| df | 2 | 每个用户非零资源块数 |
| 过载率 | 150% | J/K，表示频谱资源复用程度 |
| 单资源重叠用户数 | 3 | 每个资源块上叠加的用户数 |
这组参数的含义是：6 个用户的数据挤在 4 个资源块上发射，每个资源块上的信号由 3 个用户的信号叠加而成。df=2意味着每个用户只在 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_exp
4.1 从 MAP 到 MPA 再到 PM-MPA：三次复杂度取舍
接收端要做的事，是从叠加了多用户信号和噪声的 K 维向量里恢复每个用户的比特。最理想的是 MAP 检测，但对 6 用户 4 星座的配置，一次联合遍历就是M^J = 4096种组合，调制阶数再高就完全跑不动。MPA 的做法是在因子图上做消息传递，每个资源块只需遍历M^df种组合。
| 检测方案 | 单资源遍历量 | 复杂度特征 |
|---|---|---|
| 联合 MAP/ML | 4096 | 指数爆炸，不可实际使用 |
| 标准 MPA | 16 | 对每个用户状态遍历邻居码字组合 |
| 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会直接下溢成 0，log 再取就变成-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 或 Inf，BER 曲线出现断崖
原因：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 增大到 10，convTh从1e-3收紧到1e-5，观察曲线尾段是否改善。如果平台还在，检查信道补偿逻辑：H(k,j)在深衰落位置幅度接近 0，直接除会放大噪声，常规做法是在补偿时加一个小常数delta=1e-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 的价值。从那以后，我每次换信道模型或调制阶数，都强制走一遍“信道归一化检查、消息保护、迭代收敛性观察”这三个步骤，短则五分钟长则半小时，能替后面省下大把调参时间。希望这份拆解能让你少踩几个我踩过的坑。
本文还有配套的精品资源，点击获取
