简介:面向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 次方。拿最常见的配置举例:K=6、M=4 时,组合数是 4 的 6 次方等于 4096;码本扩到 M=8,组合数直接变成 262144。这还只是一个资源块的量级,实际仿真里一个 frame 有几十上百个符号,穷举检测从计算量上就不可能落地。
MPA(消息传递算法)把问题搬到因子图上:N 个资源节点和 K 个变量节点之间,只有码字非零位置对应的节点才连边。消息在两类节点之间交替传递,每个资源节点只需要枚举连到它身上的 df 个用户的码字组合,单节点复杂度从 M^K 降到 M^df。在 SCMA 的典型设计里,df 远小于 K,比如 K=6、N=4 时 df=3,单节点枚举量只有 4 的 3 次方等于 64 种组合,这是 MPA 能落地的根本原因。
DS-MPA 是在 MPA 基础上再做一层分解。它的思路是:当资源节点的度数 df 仍然偏大时,把连到同一资源节点的用户集合拆成两个或多个子集,子集内部先做部分消息计算,再把子集之间的结果级联起来,作为该资源节点向外传递的消息。这样单次更新的枚举规模从 M^df 降到大约 M^(df/2+1) 的量级,代价是消息不再严格精确,会引入少量近似误差,但后续迭代会让误差逐步被修正。
常见做法是:分解点选在 df 最大的那类资源节点上,用一个 split 向量记录每个子集分到几个用户。比如 df=4 拆成 2+2,df=6 拆成 3+3,拆完每个子集内部的枚举量从 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 次,背后是收敛曲线的形状决定的。在我跑过的瑞利信道仿真里,SNR=10 dB、K=6 满负载条件下,4 次迭代的 BER 和 6 次迭代大约差半个数量级,而 6 次和 10 次的曲线几乎重叠——继续加迭代只是线性增加计算量,性能上不再有可观测的改善。如果你把 max_iter 从 6 改成 10,仿真时间大约增加三分之二,BER 可能只改善一两个百分点。
迭代次数的选择其实和系统负载有关:K=4 的轻负载场景 3 次迭代就收敛了,K=8 的重负载场景 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。K=6、M=4 的配置下,每个符号周期需要枚举的组合数虽然已经通过分解降下来,但一帧 64 个符号、每个 SNR 点 1000 帧,累积起来仍然可观。我一般先把 n_frames 设为 100 跑一版粗结果,确认 BER 曲线趋势正常后,再加到 1000 以上做正式统计。这套源码的典型参数配置如下:
| 参数 | 取值 | 说明 |
|---|---|---|
| K | 6 | 用户数 |
| N | 4 | 资源元素数 |
| M | 4 | 每个用户的码本大小 |
| dv | 2 | 码字非零元素数 |
| max_iter | 6 | DS-MPA 迭代次数 |
| SNR 范围 | 0~14 dB | 步进 2 dB |
| n_frames | 1000 | 每个 SNR 点的帧数 |
| n_symbols | 64 | 每帧符号数 |
4. 避坑与常见问题排查:数值溢出、码本错位与信道归一化
4.1 log_sum_exp 相关:LLR 出现 NaN 或 Inf
现象:仿真跑到高 SNR 段,LLR 矩阵里出现 NaN 或 Inf,BER 曲线在中高信噪比处突然跳回 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 没有改善;二是码本本身的最小欧氏距离太小,某些码字对在无噪声的情况下也会混淆。
解决:先固定 SNR=14 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 曲线和 K=1 时的单用户曲线放在同一张图上。单用户没有多用户干扰,检测退化为星座点硬判决,这条曲线就是任何多用户检测器的性能下界。用同一套码本跑下来,满负载 K=6 的 DS-MPA 和单用户下界的差距通常在 2~3 dB 附近——如果实测超过 3 dB,多半是码本设计或信道归一化出了问题,回到第 4 章的排查清单过一遍。
从那以后,我每次跑完新版仿真都会强制走一遍这个流程:先画收敛曲线确认迭代次数够用,再和单用户下界对比确认差距合理,最后才把帧数拉满跑正式结果。这套流程帮我挡掉了至少三次「看起来没什么问题、实际信道归一化写错」的翻车现场,希望帮到你。
本文还有配套的精品资源,点击获取