简介：围绕稀疏码分多址（SCMA）与DS-MPA检测算法，这份MATLAB仿真资源为5G非正交多址技术研究提供了直观的代码实现。包内4个.m脚本分别负责系统仿真主程序、SCMA编码映射、DS-MPA多用户检测及数值稳定辅助计算，覆盖从发送端码字生成、瑞利衰落信道传输到接收端迭代检测的完整链路，并支持通过调整迭代次数与信道参数分析误码率、误符号率。DS-MPA算法采用分解级联思想，将复杂的联合多用户检测拆分为单用户问题的迭代求解，兼顾性能与复杂度；瑞利信道建模则贴近实际多径传播环境，便于评估算法鲁棒性。资源整体仅4KB，结构紧凑、易于阅读和修改，已有289人学习下载，适合通信工程专业学生或无线通信研究人员快速开展SCMA相关仿真实验与课题验证。
1. 从 MPA 到 DS-MPA：这份 SCMA 仿真资源到底能帮你省多少事
如果你正在搭 SCMA 链路仿真，大概率会被多用户检测卡住：6 个用户共享 4 个资源块，每个用户 4 个码字，最优联合检测要枚举 4^6=4096 种码字组合，每个资源块上还得全部遍历一遍，没跑几步计算量就爆了。工程上没人这么干，都用 MPA 这类消息传递算法在因子图上做近似。而 DS-MPA 更进一步，把每个资源块的联合检测拆成单用户级联处理，复杂度又降了一截，迭代 6 次就能收敛到接近最优的误码率。这份资源包里的四个文件——scmaenc.m、DS_MPA.m、log_sum_exp.m、simulation.m——覆盖了从比特编码、稀疏码字叠加、瑞利信道传输到对数域检测的完整链路。适合正在写 SCMA 相关论文、或者搭 5G 非正交多址仿真、需要一个能直接跑通的基线系统的人。
2. DS-MPA 检测原理：它凭什么把复杂度从指数级拉下来
2.1 因子图上的消息传递到底在传什么
要理解 DS-MPA，先得抓住 MPA 的动机。SCMA 码字是稀疏的：J 个用户共享 K 个资源块，每个用户的码字只有 N 个非零元素（N<K），所以每个资源块上只叠加了 d_f 个用户的信号。最优 MAP 检测要在每个资源块上枚举这 d_f 个用户的所有码字组合，计算量是 M^d_f。如果 d_f=3、M=4，单个资源块只要 64 次乘加，看起来不多，但乘以资源块数 K 和迭代轮数之后，整个系统的计算量会迅速膨胀，用户数再多几个就直接不可行了。
MPA 的做法是把多用户检测问题画成因子图：左侧是用户节点，右侧是资源节点，用户在某资源块上有非零码字，就画一条边。每次迭代分为两步。第一步，资源节点根据接收信号和从其他用户节点传来的消息，更新自己发给某个用户节点的消息；第二步，用户节点把从各个资源节点收到的消息乘起来，更新自己发出去的码字概率。迭代若干轮后，消息逐步收敛到用户码字的近似后验概率。
我对这套机制最大的体会是：MPA 每一步都在做局部计算，单个资源节点只需要考虑和自己相连的那 d_f 个用户，这正是稀疏码字带来的红利。如果码字不稀疏、因子图全连接，局部计算策略也救不回来。所以 SCMA 码本设计里，稀疏模式永远排在第一优先级，星座图设计反而在其次。DS-MPA 就是在 MPA 这个骨架上，把「局部计算」这一步再拆得更细。
2.2 Decomposition + Successive：DS-MPA 的两个关键动作
DS-MPA 的 DS 是 Decomposition 和 Successive 两个词的缩写。Decomposition 指的是把一个资源块上的联合检测任务拆成 d_f 个单用户子任务：处理某个用户时，把其他用户信号当作干扰，用当前已有的估计算出这个用户每个码字的概率；处理完这个用户后，立刻把它的估计从接收信号里消掉，再处理下一个用户。后面的消掉再处理，就是 Successive——串行级联。
这和并行 MPA 的差异是本质性的。并行 MPA 每轮迭代内所有用户消息同时更新，大家看到的都是上一轮的消息；DS-MPA 同一轮内后处理的用户能看到前面用户刚算出来的新消息，信息更新鲜，收敛自然更快。这也是资源包默认迭代次数只设 6 次的原因——实测里 DS-MPA 用一半轮数就能追平并行 MPA 跑 12 轮的效果。代价是串行更新不好直接做并行实现，硬件加速时得靠流水线拆。
检测器最后要给每个用户输出对数似然比 LLR，而 LLR 计算里必然出现 log(sum(exp(…))) 形态的运算，这就是 log_sum_exp.m 存在的意义。下面先把它的实现说清楚，再讲 DS_MPA 的更新顺序。
2.3 log_sum_exp：对数域计算的后悔药
function y = log_sum_exp(x) % 输入 x 是列向量，输出 log(sum(exp(x))) % 直接 exp 再求和，x 超过 700 左右就溢出成 Inf xmax = max(x); if isinf(xmax) y = xmax; else y = xmax + log(sum(exp(x - xmax))); end逻辑说明：先把输入向量里的最大值 xmax 提出来，所有元素减去它再做 exp，这样最大的项映射成 exp(0)=1，不会触发 Inf；最后结果再加上 xmax，数学上等价于原始的 log(sum(exp(x)))。如果输入里已经出现 Inf，说明上游检测出了 NaN 或溢出，这时的安全做法是直接把 xmax 当作结果返回，避免 0/0 或 Inf-Inf 这类未定义运算。
参数说明：x 是某个资源节点传给用户节点的概率对数向量，长度等于用户码字数 M，在这个资源包里 M=4。在实际调用中，x 里各元素就是「该用户取第 m 个码字」的对数概率，log_sum_exp 算的是归一化因子。少了这一步，直接 exp 求和，高信噪比下概率对数稍大一点就会爆 Inf，这是很多初版 MPA 代码跑出 NaN 的头号原因。
2.4 DS_MPA.m 的串行更新核心
DS_MPA.m 的循环结构大致是这样的：先初始化所有用户码字的对数先验，全部设为 -log(M)，然后对每个资源块、每个用户做串行消息更新，更新完一个用户就把它对相关资源块的贡献刷新一遍。核心代码片段如下：
% 串行更新：按用户顺序逐个处理 for iter = 1:max_iter for j = 1:J for k = find(F(:, j))' % 取出资源块 k 上除了用户 j 之外的其他用户 others = find(F(k, :)); others(others == j) = []; % 遍历用户 j 的 M 个码字，计算条件对数概率 for m = 1:M % 用当前其他用户的消息估计干扰项，累加得到候选概率 msg = channel_llr(k, j, m, others, Msg, H, y); Msg{k, j}(m) = msg; end % 归一化：对数域用 log_sum_exp 而不是直接求 exp Msg{k, j} = Msg{k, j} - log_sum_exp(Msg{k, j}); end end end逻辑说明：最外层循环控制迭代轮数，资源包默认 max_iter=6。中间两层分别遍历用户和该用户占用的资源块，注意串行的关键就在这——用户 j 的消息一算完，立刻写回 Msg 结构，后面用户再算自己的消息时，拿到的就是用户 j 本轮刚更新的值，而不是上一轮的旧值。这就是 Successive 的含义。
参数说明：Msg{k,j} 存的是资源块 k 传给用户 j 的消息，长度是 M=4，存的是对数概率，不是概率本身。find(F(:,j))' 取的是用户 j 的稀疏模式里非零位置对应的资源块索引。归一化那行用了 log_sum_exp 而不是 sum(exp(...))，这一步能保证 Msg 里的数值始终是对数量级，后面算 LLR 时直接做差即可，不会溢出也不会下溢成零。
3. scmaenc.m 编码链路：稀疏码字是怎么从比特里长出来的
3.1 码本结构和稀疏模式：先看非零位置怎么排
SCMA 编码不是常规的「比特→星座符号→乘扩频序列」，而是比特直接映射到一个 K 维复数码字。这个码字的特点是 K 维里只有 N 个位置非零，而且同一资源块上叠加的不同用户，非零位置是错开的，这样多个用户的信号在资源块上叠加后，才不会出现完全重叠的干扰。
资源包常见的配置是 J=6 个用户、K=4 个资源块、N=2 个非零位置、M=4 个码字。6 个用户从 4 个资源块里选 2 个，正好是 C(4,2)=6 种组合，可以实现一对一的稀疏模式映射。实践中用一张 4×6 的 0/1 矩阵 F 来表示：每行代表一个资源块，每列代表一个用户，每列恰好两个 1，每行恰好三个 1——三个 1 意味着每个资源块上同时叠加 d_f=3 个用户的信号。
这张 F 矩阵是编码和检测共用的核心数据结构。scmaenc.m 根据 F 把用户的码字放到对应的资源块上，DS_MPA.m 根据 F 决定哪些用户需要参与某个资源块的消息更新。两端引用同一个 F，改码本时就不容易两边对不上。我在拆这个包时，第一步就是把 F 打出来看，确认稀疏模式符合预期，再往下读代码。
3.2 从比特到码字：scmaenc.m 的实现骨架
function tx = scmaenc(bit, Codebook, F) % bit: J x (log2(M)*frames) 的 0/1 矩阵，J 个用户的二进制数据 % Codebook: 1xJ 的元胞数组，每个元胞是 K x M 的复矩阵 % F: K x J 稀疏模式矩阵 % tx: K x frames 的叠加发送信号 [J, total_bits] = size(bit); K = size(F, 1); M = size(Codebook{1}, 2); logM = log2(M); frames = total_bits / logM; tx = zeros(K, frames); for j = 1:J bits_j = bit(j, :); idx = bi2de(reshape(bits_j, logM, frames).', 'left-msb') + 1; for f = 1:frames tx(:, f) = tx(:, f) + Codebook{j}(:, idx(f)); end end逻辑说明：外层循环按用户处理，每个用户每 log2(M) 个比特转换成一个十进制码字索引，然后去码本 Codebook{j} 里取第 idx 列。这个列本身就是 K 维复数向量，非零位置与 F 的第 j 列一致。把该用户的码字加到发送信号 tx 对应帧上，所有用户处理完后，tx 就是叠加后的 SCMA 发送信号。
参数说明：bi2de 是 Communications Toolbox 的函数，没有这个工具箱时可以自己实现：idx = bits_j(1)*2 + bits_j(2) + 1（假设 M=4、logM=2）。帧数 frames 由总比特数除以每符号比特数得到，仿真时一般设成几千到几万帧，用来累积统计误码率。Codebook{j}(:, idx) 取的是复数码字，不是实数，用户叠加也是复数域叠加，后面瑞利信道衰落的实部虚部都要一起处理。
3.3 码本功率归一化：第一眼看不出来、跑起来才头疼的坑
码本里每个码字的功率必须提前归一化。具体做法是让每个码字列向量的非零元素模平方和等于同一个常数，常见做法是归一化到 1。如果不做这一步，不同用户的码字幅度参差不齐，仿真里的信噪比定义就会全乱套——某个用户统计 BER 时等效 SNR 比别人偏了 1~2 dB，最后的平均误码率曲线会被拉偏，跟已发表论文对不上。
我一般会在读码本时顺手算一下每列的 norm，打印出来看是不是一致。公开代码里拿来的码本，十个有八个没做这步。如果你发现仿真曲线整体偏移、或者高 SNR 区出现诡异平台，优先检查码本功率，别急着改迭代次数。
4. simulation.m 完整仿真：瑞利信道下把 BER 曲线跑出来
4.1 主循环骨架
simulation.m 做的事情可以概括成一段话：设定系统参数，初始化码本和稀疏模式，然后在一系列 SNR 点上，反复执行「生成随机比特→编码→过瑞利信道→加高斯噪声→DS-MPA 检测→统计误码」这个循环。核心结构如下：
J = 6; K = 4; M = 4; max_iter = 6; EbN0dB = 0:2:16; num_frames = 5000; for snr_idx = 1:length(EbN0dB) EbN0 = 10^(EbN0dB(snr_idx)/10); N0 = 1 / (log2(M) * EbN0); noise_var = N0 / 2; total_bit_errors = 0; total_bits = 0; for frame = 1:num_frames bit = randi([0 1], J, log2(M)); tx = scmaenc(bit, Codebook, F); h = (randn(K, J) + 1j*randn(K, J)) / sqrt(2); rx = sum(h .* tx, 2) + sqrt(noise_var) * (randn(K,1) + 1j*randn(K,1)); llr = DS_MPA(rx, h, Codebook, F, max_iter); % 硬判决 + 统计误码 est_bits = llr < 0; total_bit_errors = total_bit_errors + sum(bit(:) ~= est_bits(:)); total_bits = total_bits + numel(bit); end BER(snr_idx) = total_bit_errors / total_bits; end逻辑说明：外层循环遍历 SNR 点，内层循环跑固定帧数。每帧内先随机生成 6 个用户的比特，编码成发送信号 tx；随后生成瑞利信道矩阵 h，尺寸是 K×J，表示每个用户在每个资源块上的信道增益；接收端把信道作用后的信号叠加，再加上复高斯噪声。DS_MPA 检测完输出 LLR，硬判决后统计误码。
参数说明：num_frames 设 5000 是折中——帧太少 BER 抖动大，帧太多仿真时间翻倍。高 SNR 区误码率低，想看到 1e-4 量级至少要几万帧，我一般会先跑 2000 帧验证流程，确认无误后再加大到 20000 帧出图。
4.2 瑞利信道系数：给每个资源块一个复衰落
瑞利信道的核心是生成服从复高斯分布的衰落系数：h = (randn + 1j*randn) / sqrt(2)。这样生成的 h，模平方的期望是 1，也就是信道平均不放大也不衰减信号，方便后面按 Eb/N0 换算噪声功率。
这里的细节是 h 的维度。SCMA 上行链路里，每个用户的信号只落在自己码字的非零资源块上，而 F 矩阵记录了这个对应关系。simulation.m 里直接生成 K×J 的满矩阵，让 DS_MPA 在更新消息时根据 F 取用对应位置的 h，这个做法实现起来简单，但要注意：F 里为零的位置，h 就算生成了也不该参与计算，否则会把不存在的信道路径当成有效信号处理。
4.3 Eb/N0 与噪声功率换算：整个仿真最容易算错的公式
噪声功率的换算我单独拿出来说，因为这里翻车概率最高。发送端码本每列归一化到 1，那么每个用户每符号能量是 1，每比特能量 Eb = 1/log2(M)。给定 Eb/N0，噪声单边功率谱密度 N0 = 1/(log2(M)*EbN0)。复噪声是二维的，实部虚部各占 N0/2 的方差，所以加噪时要用sqrt(N0/2)*(randn+1j*randn)，而不是直接乘sqrt(N0)。
常见的错误是把复噪声总方差当成 N0，结果等效 SNR 高了 3 dB，BER 曲线比正确结果偏优，跟论文对不上时第一个怀疑的就是这里。另一个经常被忽略的点是：Eb/N0 定义里 Eb 是每用户每比特能量，多用户叠加后的总接收功率是 J 倍，但噪声功率不变，所以用户数变化时 BER 曲线位置会移动，拿不同用户数的结果放同一张图对比时，必须要确认两边用的 Eb/N0 口径一致。
5. 避坑指南：跑这套 DS-MPA 代码最容易翻车的五个地方
5.1 现象：高信噪比下 BER 突然变成 NaN
原因：LLR 计算里出现了 log(0) 或者 0/0。高 SNR 下某个码字后验概率会极端集中，正常码字的概率可能下溢成 0，取对数就是 -Inf，再参与运算就变成 NaN。
解决：所有概率运算在取对数之前加一个极小量，比如log(p + 1e-12)；对数域里则用 log_sum_exp 统一做归一化，避免直接构造概率再取对数。我检查代码时还会顺手把中间变量里的 NaN 位置打出来，定位到具体是哪个资源块哪次迭代出的问题。
5.2 现象：迭代次数从 6 加到 12，BER 纹丝不动
原因：消息初始化写错了。把先验对数概率初始化为 0，等于认为每个码字的概率是 exp(0)=1，归一化之前概率总和是 4 而不是 1，消息更新的归一化因子全错，迭代收敛到了错误的不动点。
解决：初始化改为-log(M)，也就是四个码字均匀分布的对数概率。改完这行，BER 曲线通常立刻恢复正常下降趋势。这个坑我见过不止一次，症状很有迷惑性——低 SNR 时影响不明显，高 SNR 时曲线就是压不下去。
5.3 现象：DS-MPA 跑出来的结果和普通并行 MPA 完全一样
原因：看起来串行的代码实际还是并行更新。常见写法是先把所有用户的新消息算到临时变量里，一轮结束后统一写回 Msg 结构，这就把串行退化成了并行。
解决：检查用户循环内部是否直接对 Msg{k,j} 就地赋值，后一个用户读取 Msg 时是否包含前面用户刚更新的值。DS-MPA 的性能优势全靠这个「新消息即时可见」，写成临时变量中转就等于丢掉了 Successive 的核心。
5.4 现象：瑞利信道下 BER 曲线比已发表论文差 2~3 dB
原因：瑞利信道系数没归一化。直接randn + 1j*randn生成的复高斯随机变量功率期望是 2，会让信号被平均放大 3 dB，但噪声没变，等效 SNR 偏低，曲线自然右移。
解决：除以 sqrt(2)，让 E(|h|^2)=1。我建议在仿真里打印一次mean(abs(h(:)).^2)确认这个值接近 1，每次改信道模型都顺手验一下。
5.5 现象：换一个码本后 BER 性能大变，甚至不如随机星座
原因：码本非零位置和 F 矩阵不匹配，或者码本功率没归一化。SCMA 码本设计本身就很敏感，非零位置的选取直接决定因子图结构，换码本时 F 必须同步更新。
解决：码本放进 Codebook cell 前，先逐列检查非零位置和 F 的第 j 列是否一致、每列能量是否相等、星座各维度是否在合理范围内。我一般会在 scmaenc.m 里加一个断言，发现不匹配立即报错，免得仿真跑了半小时才发现在用错码本。
6. 进阶：用迭代扫描验证 DS-MPA 的收敛边界，别把 6 次迭代当玄学
资源包默认迭代次数是 6，但 6 不是从天上掉下来的。不同信道条件、不同用户数下，DS-MPA 收敛需要的迭代次数会变。我在跑正式仿真前，通常会做一次收敛性扫描：
EbN0_fixed = 10; % 固定在中等信噪比观察收敛趋势 num_frames = 10000; for iter_test = 1:10 % 固定随机种子，保证每个迭代次数跑的是同一批数据 rng(42); % 完整仿真循环，detector 内部 max_iter 换成 iter_test BER_iter(iter_test) = run_simulation(EbN0_fixed, num_frames, iter_test); end plot(1:10, log10(BER_iter), '-o');固定随机种子这一步是很多人忽略的关键：帧数有限时 BER 本身有随机抖动，不固定种子，不同迭代次数的 BER 差异里混入了噪声，很难看清真正的收敛趋势。固定 rng(42) 之后，BER 随迭代次数的变化就纯粹来自算法本身。脚本跑完，看曲线在哪个迭代次数之后不再明显下降，那个点就是当前配置的收敛边界。如果 6 次时曲线还在陡降，就要把默认迭代次数调上去；如果 3 次就平了，6 次纯属浪费算力。
我在一次改信道模型时遇到过高 SNR 下 BER 平台降不下去的情况，直觉以为迭代不够，扫了一圈发现是信道系数没归一化。从那以后，我每次调迭代次数都强制先做一遍收敛扫描，顺手把 h 的功率也打印出来。这套习惯帮我省掉了大量排查时间。希望帮到你。
本文还有配套的精品资源，点击获取
