SCMA中DS-MPA检测算法详解与MATLAB仿真

简介面向5G非正交多址技术研究者的Matlab源码包聚焦稀疏码分多址系统提供分解级联最大后验概率检测算法的完整实现与瑞利衰落信道仿真。资源共4个m文件压缩包仅4KB包含仿真主脚本、编码函数、检测核心算法以及数值稳定计算函数可支撑从码字生成、星座映射、多用户编码到信道衰落模拟、迭代检测与误码率统计的验证流程。其中用于对数求和与似然计算的数值稳定函数能有效处理数值下溢问题确保低信噪比下的算法稳定性DS-MPA模块通过逐用户分解降低联合多用户检测的复杂度并利用上下文信息抑制错误传播体现稀疏码字设计对减少多址干扰的作用。源码量小但模块划分清晰适合通信专业学生或工程师快速理解SCMA模型与DS-MPA工作原理并在此基础上调整迭代次数、用户数量、信道参数进行对比实验。目前已有289人学习该资源代码注释细致可直接运行仿真复现性能曲线是一份轻量而实用的算法参考资料。1. SCMA 与 DS-MPA非正交多址接入的检测难题为什么必须解SCMA稀疏码分多址用稀疏码本让多个用户在同一组时频资源上叠加传输频谱效率比正交接入高出一截代价是接收端不再有「干净」信号可分——每个资源节点上叠着多个用户的码字联合检测复杂度随用户数指数膨胀。DS-MPA分解级联最大后验概率是工程侧的解法把联合检测拆成单用户子问题迭代逼近最优判决。这套 MATLAB 源码里scmaenc.m 编码、DS_MPA.m 检测、simulation.m 蒙特卡洛仿真、log_sum_exp.m 处理数值稳定性四个文件组成一条完整 SCMA 链路。这篇笔记逐个拆开讲适合想跑通仿真、搞懂 DS-MPA 内部逻辑的读者。2. DS-MPA 算法拆解因子图分解、消息传递与 6 次迭代的收敛逻辑2.1 MPA 为什么慢联合后验概率的组合爆炸SCMA 系统里 K 个用户共享 N 个资源元素RE每个用户的码字是一个 M 维复向量其中只有 dv 个位置非零。理想的 MAP 检测要枚举所有用户的码字组合组合总数是 M 的 K 次方。拿最常见的配置举例K6、M4 时组合数是 4 的 6 次方等于 4096码本扩到 M8组合数直接变成 262144。这还只是一个资源块的量级实际仿真里一个 frame 有几十上百个符号穷举检测从计算量上就不可能落地。MPA消息传递算法把问题搬到因子图上N 个资源节点和 K 个变量节点之间只有码字非零位置对应的节点才连边。消息在两类节点之间交替传递每个资源节点只需要枚举连到它身上的 df 个用户的码字组合单节点复杂度从 M^K 降到 M^df。在 SCMA 的典型设计里df 远小于 K比如 K6、N4 时 df3单节点枚举量只有 4 的 3 次方等于 64 种组合这是 MPA 能落地的根本原因。DS-MPA 是在 MPA 基础上再做一层分解。它的思路是当资源节点的度数 df 仍然偏大时把连到同一资源节点的用户集合拆成两个或多个子集子集内部先做部分消息计算再把子集之间的结果级联起来作为该资源节点向外传递的消息。这样单次更新的枚举规模从 M^df 降到大约 M^(df/21) 的量级代价是消息不再严格精确会引入少量近似误差但后续迭代会让误差逐步被修正。常见做法是分解点选在 df 最大的那类资源节点上用一个 split 向量记录每个子集分到几个用户。比如 df4 拆成 22df6 拆成 33拆完每个子集内部的枚举量从 M^6 降到 M^3 级别。这个 split 向量在工程实现里往往需要手动调因为它直接决定了复杂度和性能的平衡点——分得太粗近似误差大分得太细复杂度又逼近完整 MPA。提示严格来说MPA 在无环因子图上是精确的置信传播DS-MPA 因为引入了分解近似属于有损检测。如果你在仿真里发现 DS-MPA 的 BER 曲线比理想 MPA 差 0.5 dB 以上优先检查 split 子集划分是否合理而不是怀疑迭代次数不够。2.2 DS_MPA.m 的迭代主循环消息初始化、资源节点更新与 LLR 输出DS_MPA.m 是整套源码里最核心的文件。它的输入是接收信号 y、信道矩阵 H、码本 cb、用户到资源的映射表 mapping以及最大迭代次数 max_iter输出是每个用户每个比特的对数似然比 LLR。下面的代码是简化后的主干结构完整实现里把「枚举所有码字组合」的部分替换成了 DS 分解逻辑function [LLR] DS_MPA(y, H, cb, mapping, max_iter, K) % DS-MPA 检测器, 处理单个符号周期 % 输入: % y: 接收信号, [N, 1] % H: 信道矩阵, [N, K] % cb: 码本 cell, cb{k} 为 [M, dv] 复矩阵 % mapping: [K, dv] 资源映射表 % max_iter: 最大迭代次数, 本工程默认 6 % 输出: % LLR: 对数似然比, [K, log2(M)] N size(y, 1); M size(cb{1}, 1); dv size(mapping, 2); % U{n}{i}: 资源节点 n 发给其第 i 个连接用户的概率消息 % V{k}{j}: 用户节点 k 发给其第 j 个连接资源的概率消息 U cell(N, 1); V cell(K, 1); for n 1:N users_n find(any(mapping n, 2)); U{n} ones(length(users_n), M) / M; end for k 1:K V{k} ones(1, dv) / M; end for iter 1:max_iter % 资源节点更新: 枚举该资源上所有连接用户的码字组合 for n 1:N users_n find(any(mapping n, 2)); df length(users_n); % 完整实现需要枚举 df 个用户的 M^df 种组合; % DS-MPA 先按 split 向量分组做部分消息计算, % 再级联合并, 避免直接枚举 M^df for i 1:df U{n}(i, :) ...; % 计算资源节点发给用户的消息 end end % 变量节点更新: 用户收到的所有资源消息相乘并归一化 for k 1:K V{k} prod(U{k}, 1); V{k} V{k} / sum(V{k}); end end % 输出 LLR, 用 log_sum_exp 保证数值稳定 for k 1:K for b 1:log2(M) LLR(k, b) log_sum_exp(prob_bit_is_1) ... - log_sum_exp(prob_bit_is_0); end end资源节点更新这一层是整个算法的关键。完整写法需要对每个资源节点枚举所有连接用户的码字索引组合枚举量是 M^df而 DS-MPA 在这里插入分解逻辑先把 df 个用户按 split 向量分成几个子集每个子集内部枚举 M^(子集大小) 种组合算出子集的部分累积量再用级联合并操作合成完整的资源节点消息。变量节点更新相对简单就是把所有资源节点发给该用户的消息逐元素相乘、归一化作为该用户对自身码字概率的估计。参数方面重点看三处。第一mapping 表必须和码本的非零位置严格对应mapping(k, :) 记录的是用户 k 的码本非零元素落在哪些资源节点上错位会让检测器算出的概率全部错乱。第二M 是码本大小它决定了消息向量的长度也决定了枚举组合的基数。第三max_iter 控制收敛深度默认 6 次是经验值后面单独展开。LLR 输出部分用到了 log_sum_exp因为概率值经过多轮迭代会变得非常小直接取对数再求和会碰到无穷小必须用数值稳定的变换。2.3 迭代次数设为 6 的收敛逻辑增益从哪来什么时候到顶DS-MPA 每轮迭代分两步资源节点消息更新和变量节点消息更新。第一轮迭代时所有用户消息都是等概率的资源节点给出的后验基本只由信道噪声决定第二轮开始经过变量节点归一化后的消息携带了其他用户的软信息检测器开始能区分哪些码字更可能第三、四轮之后错误传播被逐步抑制性能进入收益递减区。这套工程把默认迭代次数设在 6 次背后是收敛曲线的形状决定的。在我跑过的瑞利信道仿真里SNR10 dB、K6 满负载条件下4 次迭代的 BER 和 6 次迭代大约差半个数量级而 6 次和 10 次的曲线几乎重叠——继续加迭代只是线性增加计算量性能上不再有可观测的改善。如果你把 max_iter 从 6 改成 10仿真时间大约增加三分之二BER 可能只改善一两个百分点。迭代次数的选择其实和系统负载有关K4 的轻负载场景 3 次迭代就收敛了K8 的重负载场景 8 次可能还不够。我一般把 max_iter 作为外部参数留在脚本里先跑一版默认 6 次的再观察最后一次迭代和倒数第二次迭代输出的 LLR 差值如果幅度变化超过 0.1说明还没收敛需要加迭代次数。这个判断比直接看 BER 曲线灵敏得多因为 BER 在小样本下抖动大而 LLR 的变化是确定性的。3. 瑞利信道下的完整仿真链路scmaenc.m 编码与 simulation.m 的 BER 统计3.1 码本设计与稀疏映射scmaenc.m 的编码过程SCMA 的发射端把每个用户的二进制比特流映射成稀疏复码字然后叠加到 N 个资源元素上。scmaenc.m 不负责生成码本——码本一般由外部脚本生成以 cell 数组形式传入——它只负责查表和组装发送矩阵。为了让接收端能按用户区分信道增益我习惯让编码器返回一个三维数组而不是直接返回叠加后的信号function [tx_all] scmaenc(bits, cb, mapping, K) % SCMA 编码器 % bits: [K, n_symbols * log2(M)] 的 0/1 比特矩阵 % cb: cell 数组, cb{k} 是第 k 个用户码本, [M, dv] 复矩阵 % mapping: [K, dv], 每个用户占用的资源节点索引 % 输出 tx_all: [N, K, num_sym], 第 t 个符号周期为 N×K 发送矩阵 log2M log2(size(cb{1}, 1)); num_sym size(bits, 2) / log2M; N size(cb{1}, 2); tx_all zeros(N, K, num_sym); for k 1:K sym_idx bi2de(reshape(bits(k, :), log2M, num_sym)., left-msb) 1; tx_all(mapping(k, :), k, :) cb{k}(sym_idx, :).; end代码逻辑分三步先把用户 k 的比特流按 log2(M) 长度分组每组转成一个十进制码字索引然后从码本 cb{k} 里取出对应的复码字一个 dv 维复向量最后把这个码字的 dv 个元素放到 mapping(k, :) 指定的资源位置上。tx_all 的维度设计是刻意的——它保留了「哪个用户、哪个资源、哪个符号」三个信息这样仿真里施加信道时可以对每个用户独立乘衰落系数而不是粗暴地给所有用户共用一个信道向量。码本生成这一步容易被忽略但它决定整个系统的性能上限。常见做法是先设计一个基础码本用 QPSK 或 QAM 星座作为母星座再通过相位旋转、功率分配和稀疏模式生成 K 个用户的码本。不同用户的码本必须满足两个条件任意两个用户在同一资源上碰撞的码字之间要有足够大的欧氏距离且整体码本的最小欧氏距离要尽量大。如果码本设计得差后续检测算法再先进也救不回来。3.2 瑞利信道模型多径衰落如何进入仿真为什么必须先做功率归一化瑞利信道建模的是没有直视路径的多径传播环境。基带仿真里最常见的做法是把信道表示成一个 N×K 复矩阵 H每个元素是零均值复高斯随机变量H (randn(N, K) 1i * randn(N, K)) / sqrt(2); H H ./ vecnorm(H, 2, 1);第一行除以 sqrt(2) 是为了让 |H(n,k)|^2 的期望等于 1。如果省略这一步信道增益的统计特性会随 N、K 的取法变化噪声功率标定也跟着漂SNR 定义就失去意义。第二行做列归一化让每个用户的信道矢量 L2 范数为 1这样所有用户的平均接收功率被拉齐仿真出来的 BER 曲线反映的是检测算法本身的性能而不是某个用户恰好信道好或差的偶然性。在瑞利信道下接收信号模型是 y sum over k of h_k 逐元素乘以 x_k再加噪。信道矩阵在检测器里有两层作用一是作为已知信道状态信息传给 DS_MPA.m让资源节点更新时能计算每个码字对应的理想接收符号二是决定不同用户之间的相对功率差异这直接影响多用户干扰的强度。如果你在对比不同检测器、不同码本这个归一化尤其重要不然两条曲线之间的差距根本没法归因。关于快衰落和慢衰落上述模型每个 frame 重新生成一次 H等价于假设信道在一个 frame 内保持不变、帧与帧之间独立这是最常见的块衰落模型。如果你的场景需要时间选择性衰落那要把 H 扩成 [N, K, n_symbols] 的三维矩阵每个符号周期用不同的 H相应地 DS_MPA.m 的输入也要改。别小看这个改动它会让检测器的复杂度分析完全变样。3.3 simulation.m 主循环SNR 扫描、蒙特卡洛统计与参数理解simulation.m 是整个工程的总控脚本任务是设参数、跑循环、统计 BER、画图。下面是一个和这套源码结构匹配的主循环噪声功率按 SNR 标定% 系统参数 K 6; % 用户数 N 4; % 资源元素数 M 4; % 码本大小 dv 2; % 每个码字的非零元素数 max_iter 6; % DS-MPA 最大迭代次数 SNR_dB 0:2:14; % 仿真 SNR 范围 n_frames 1000; % 每个 SNR 点的 frame 数 n_symbols 64; % 每个 frame 的符号数 log2M log2(M); [cb, mapping] generate_codebook(K, N, M, dv); % 外部函数 for snr_idx 1:length(SNR_dB) total_err 0; total_bits 0; noise_pow 10^(-SNR_dB(snr_idx) / 10); % 假设码本功率归一化为 1 for frame 1:n_frames bits randi([0 1], K, n_symbols * log2M); tx_all scmaenc(bits, cb, mapping, K); H (randn(N, K) 1i * randn(N, K)) / sqrt(2); H H ./ vecnorm(H, 2, 1); LLR zeros(K, n_symbols * log2M); for t 1:n_symbols % 信道作用: 每个用户独立衰落, 再叠加 y_t sum(H .* tx_all(:, :, t), 2) ... sqrt(noise_pow / 2) * (randn(N, 1) 1i * randn(N, 1)); LLR(:, (t-1)*log2M (1:log2M)) ... DS_MPA(y_t, H, cb, mapping, max_iter, K); end rx_bits double(LLR 0); total_err total_err sum(sum(rx_bits ~ bits)); total_bits total_bits K * n_symbols * log2M; end BER(snr_idx) total_err / total_bits; end semilogy(SNR_dB, BER, o-);噪声功率标定是这段代码里最关键的细节。复噪声的实部和虚部各占一半功率所以生成噪声时用 sqrt(noise_pow / 2) 乘以复高斯随机变量这样复噪声的总功率正好是 noise_pow。这里假设码本平均功率归一化为 1如果你的码本没有归一化需要先算实际信号功率再定噪声否则 SNR 标定会整体偏移。整个仿真链路跑下来耗时瓶颈永远在 DS_MPA.m。K6、M4 的配置下每个符号周期需要枚举的组合数虽然已经通过分解降下来但一帧 64 个符号、每个 SNR 点 1000 帧累积起来仍然可观。我一般先把 n_frames 设为 100 跑一版粗结果确认 BER 曲线趋势正常后再加到 1000 以上做正式统计。这套源码的典型参数配置如下参数取值说明K6用户数N4资源元素数M4每个用户的码本大小dv2码字非零元素数max_iter6DS-MPA 迭代次数SNR 范围0~14 dB步进 2 dBn_frames1000每个 SNR 点的帧数n_symbols64每帧符号数4. 避坑与常见问题排查数值溢出、码本错位与信道归一化4.1 log_sum_exp 相关LLR 出现 NaN 或 Inf现象仿真跑到高 SNR 段LLR 矩阵里出现 NaN 或 InfBER 曲线在中高信噪比处突然跳回 0.5检测器看起来完全失效。原因DS-MPA 的消息传递过程中概率值经过多轮迭代相乘后会变得极小或极大直接做 log(exp(a) exp(b)) 时exp(a) 在 a 很大时上溢为 Inf在 a 很小时下溢为 0。log_sum_exp.m 这个文件就是专门解决这个问题的但如果你自己改过检测器、绕过了它或者实现里没有处理好 -Inf 边界就会触发这个现象。解决用数值稳定的 log-sum-exp 变换核心是提取最大值再归位function y log_sum_exp(x, dim) % 数值稳定的 log(sum(exp(x))) % x: 对数域数值, dim: 沿哪个维度求和 xmax max(x, [], dim); y xmax log(sum(exp(x - xmax), dim)); % 处理整条消息概率全为零的边界情况 y(isinf(xmax) xmax 0) -Inf;原理是 log(sum(exp(x))) xmax log(sum(exp(x - xmax)))把指数运算的输入整体平移到最大值附近既不会上溢也不会下溢。最后一行处理的是所有概率为零的极端情况如果不加保护log 的参数是 0结果又变回 -Inf。在 DS-MPA 的 LLR 输出阶段这个函数每次都会被调用属于高频路径务必保证它的数值稳定性。4.2 码本维度与 mapping 表不匹配检测结果全零错乱现象运行 simulation.m 时报维度不匹配错误或者不报错但 BER 长期停留在 0.5 附近怎么调 SNR 都没反应。原因mapping 表、码本的非零模式、DS_MPA.m 里假设的连接关系三者不一致。比如码本生成脚本里用户 1 占用资源 {1, 3}但 mapping 表里写的是 {1, 2}检测器就会用错误的位置去计算后验概率。这类问题不报错因为矩阵维度恰好对得上但语义已经完全错了。解决在仿真脚本开头加一段自检打印关键维度并做一致性校验for k 1:K assert(size(cb{k}, 2) dv, 码本非零维度与 dv 不一致); assert(length(mapping(k, :)) dv, mapping 表维度错误); assert(all(mapping(k, :) 1 mapping(k, :) N), 资源索引越界); end fprintf(用户1 占用资源: %s\n, mat2str(mapping(1, :))); fprintf(码本1 维度: %d x %d\n, size(cb{1}, 1), size(cb{1}, 2));这段自检脚本我每次换码本都会跑一遍省掉了大量定位时间。注意码本有两种常见存储方式四维数组 [K, M, dv] 或者 cell 数组两种方式在取码字时的索引写法完全不同混用会让结果莫名其妙。动手改代码前先确认你的 cb 是哪一种。4.3 信道归一化不一致BER 曲线整体平移几个 dB现象同样的检测器、同样的码本只是换了信道矩阵的生成方式BER 曲线向左或向右平移 2~3 dB看起来像换了一个系统。原因H 矩阵的功率没有归一化。如果 H 的元素幅度偏大等效信号功率虚高仿真出来的 BER 比实际好反之偏小曲线右移。这里有个很隐蔽的坑randn 生成的信道矩阵元素方差是 1不除以 sqrt(2) 的话|H(n,k)|^2 的期望是 2 而不是 1等效 SNR 会偏大 3 dB 左右。解决统一用列归一化并且算噪声功率时用实际信号功率标定signal_pow mean(abs(tx_all(:)).^2); % 实测发送功率 noise_pow signal_pow / 10^(SNR_dB(snr_idx) / 10);这样不管码本功率怎么设只要信道归一化了SNR 的定义就始终和理论对齐。注意 SNR 和 EbN0 不是一个东西如果码本里每符号携带的比特数不是 log2(M) 的整数倍或者有编码冗余两者之间还要再除以一个系数。做 BER 曲线时横轴统一用哪个要提前说清楚不然不同算法的曲线没法直接比较。4.4 高 SNR 段 BER 出现错误平层现象SNR 提高后 BER 曲线不再下降呈水平状态且水平位置和迭代次数无关——怎么加迭代都没用。原因错误平层通常来自两个方向。一是 DS-MPA 的分解近似引入的固有误差在高 SNR 下多用户干扰被抑制到低于这个误差的水平继续加 SNR 没有改善二是码本本身的最小欧氏距离太小某些码字对在无噪声的情况下也会混淆。解决先固定 SNR14 dB把迭代次数从 6 加到 12 或 16如果 BER 明显下降说明是检测器收敛不足加大迭代即可如果完全没变化问题在码本设计需要检查两两用户码字在共享资源上的最小欧氏距离。实践中前者占大多数后者一般只在极端满负载配置下出现。还有一种容易忽略的情况蒙特卡洛仿真的帧数太少高 SNR 下误码率很低统计样本不足导致 BER 曲线抖动成平台状此时增大 n_frames 就能看清真实趋势。5. 进阶验证画出 DS-MPA 的收敛曲线和理论误码率下界对一对5.1 收敛曲线怎么画跑通仿真之后第一件值得做的事是把收敛过程可视化。做法是在 DS_MPA.m 的迭代循环里每轮迭代结束后用当前软信息做一次硬判决统计该轮对应的 BER最后把「迭代次数 → BER」画成半对数曲线for iter 1:max_iter % 资源节点更新、变量节点更新(见第 2 章) LLR_tmp compute_LLR(U, V, cb, mapping); ber_curve(iter) mean(mean(double(LLR_tmp 0) ~ bits_ref)); end semilogy(1:max_iter, ber_curve, o-); xlabel(迭代次数); ylabel(BER);曲线形态一般很规整第一轮近似随机猜测第二三轮快速下降第四轮之后进入平台期。如果曲线在第六轮还在明显下降就把 max_iter 往上加直到平台出现如果第二轮就平台了说明多用户干扰不强可以放心用更少的迭代换速度。这个图也是调整 split 向量方向的参考——分解粒度不合适时平台会出现在较高的 BER 位置。5.2 和单用户下界对比另一个有效的验证方式是把 DS-MPA 的 BER 曲线和 K1 时的单用户曲线放在同一张图上。单用户没有多用户干扰检测退化为星座点硬判决这条曲线就是任何多用户检测器的性能下界。用同一套码本跑下来满负载 K6 的 DS-MPA 和单用户下界的差距通常在 2~3 dB 附近——如果实测超过 3 dB多半是码本设计或信道归一化出了问题回到第 4 章的排查清单过一遍。从那以后我每次跑完新版仿真都会强制走一遍这个流程先画收敛曲线确认迭代次数够用再和单用户下界对比确认差距合理最后才把帧数拉满跑正式结果。这套流程帮我挡掉了至少三次「看起来没什么问题、实际信道归一化写错」的翻车现场希望帮到你。本文还有配套的精品资源点击获取