ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

SCMA的MPA译码Matlab仿真:从码本设计到误码率曲线

SCMA的MPA译码Matlab仿真:从码本设计到误码率曲线

简介:基于MPA的SCMA的Matlab仿真项目是一份面向5G非正交多址接入技术研究的可运行代码包,适合通信工程专业学生、算法研究人员与系统工程师深入理解稀疏码多址原理和消息传递迭代译码实现。压缩包共包含11个文件,其中9个m脚本覆盖了SCMA编码、稀疏码本设计、多用户接入、MPA译码及误码率统计等核心环节,并配有一份md格式说明文档和一份txt说明文件,方便读者快速梳理仿真流程与函数调用关系。整个资源包仅11KB,非常轻量,在Matlab环境下即可直接运行和二次修改。该仿真项目目前已有1675人学习使用。通过研读源码与运行结果,读者可以掌握MPA因子图消息更新的具体编程路径,观察不同信噪比下的误码性能,并可进一步调整迭代次数、码本尺寸等参数,为优化SCMA系统或开展多址接入技术对比提供实验支撑,是研究和教学场景中一份实用的参考实现。

1. 为什么还要自己拉一遍「基于MPA的SCMA的Matlab仿真」

「基于MPA的SCMA的Matlab仿真」这串词,看起来像通信课的一道大作业,真正动手把它跑通的人都知道,里面的坑比想象中多不少。SCMA(稀疏码分多址)靠码本稀疏性换来过载能力,MPA(消息传递算法)则靠因子图迭代把这笔非正交的账算回来。这篇内容要做的事很具体:从系统参数、码本、发射机到MPA迭代译码,搭出一条能出误码率曲线的完整仿真链路。适合通信方向做链路级仿真的研究生、想做非正交多址对比的工程师,以及所有听过SCMA但没亲手写过消息传递代码的人。跑通之后你会发现,真正吃时间的不是编码,而是那些藏在信噪比换算和因子图索引里的玄学。

2. 仿真骨架:SCMA链路怎么搭,先让发射端动起来

2.1 系统模型与参数表:先定过载率、码本大小和稀疏模式

SCMA的系统参数一般围绕用户数 J、资源块数 K、码本大小 M、每个码字的非零资源数 df 来定。业内做这类仿真最常见的配置是 J=6、K=4、M=4、df=2,过载率 J/K=150%。这个配置的复杂度正好卡在「能让MPA跑得动」的范围:每个资源块上只有 df·J/K=3 个用户叠加,功能节点更新时的组合遍历量是 M²=16,代码写起来不痛苦,跑一帧也就毫秒级。

参数先列成表,后面所有代码都围绕这张表展开:

参数含义演示取值
J用户层数6
K资源块数4
M码本大小(每用户候选码字数)4
df每个码字中的非零符号数2
dc每个资源块上叠加的用户数3
λ过载率150%
EbN0dB每比特信噪比扫描范围0:2:12

稀疏模式决定了因子图的长相。演示系统里每个资源块恰好被3个用户占用,每个用户占用2个资源块,我用一个 K×J 的0/1矩阵 smap 来表达,行是资源块、列是用户。这个矩阵在发射端和译码端都要用到,建议在一开始就定义成全局唯一的变量。

J = 6; K = 4; M = 4; df = 2; EbN0dB = 0:2:12; % 稀疏模式: 行=资源块, 列=用户 % 用户1占资源1,2; 用户2占资源1,3; ... 每列恰好2个1 smap = [1 1 1 0 0 0; 1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1];

这段代码里的 smap 不是随便写的,它保证了每一列恰好有2个非零、每一行恰好有3个非零,也就是 df 和 dc 都均衡。用这个矩阵去画因子图,就是6个变量节点、4个功能节点、12条边。注意 smap 的列顺序就是用户编号顺序,译码端如果换了用户顺序,消息传递会整个错位,后面专门讲这个坑。一般做法是把 smap 和 J、K、M 一起放在一个结构体里,编码和译码都从同一个结构体取参数,避免两头维护。

2.2 演示码本与SCMA编码:把比特映射成稀疏码字

SCMA的码本是一个 K×M×J 的三维复数数组,cb(:, m, j) 表示用户 j 的第 m 个码字。每个码字在 df 个资源上有非零符号,其余位置是0。码本设计是SCMA性能的关键,直接关系到同一资源块上叠加星座的可分性。正式研究中码本通常要做欧氏距离优化,但做链路验证时可以用一个固定结构的演示码本:基础QPSK星座加每用户相位旋转。

% make_demo_codebook.m % 生成可复现的演示码本: 基准QPSK + 每用户相位旋转 % 注意: 该码本未做距离优化, 仅用于链路调试 function cb = make_demo_codebook() J = 6; K = 4; M = 4; base = [1+1i; -1+1i; -1-1i; 1-1i] / sqrt(2); % 平均符号能量=1 cb = zeros(K, M, J); for j = 1:J theta = (j - 1) * pi / 6; % 每用户一个旋转角 sym = base * exp(1i * theta); % 旋转后的符号集 pos = find(smap(:, j)); % 用户j的两个非零资源 cb(pos(1), :, j) = sym.'; cb(pos(2), :, j) = sym.'; end end

逻辑说明:每个用户在两个非零资源块上放置同一组旋转后的QPSK符号,旋转角度随用户索引递增。这样做的目的是让不同用户在同一资源块上的星座不再完全重合,MPA才有足够的区分度;如果所有用户都用同一组QPSK符号,叠加后会出现大量完全相同的候选组合,译码性能会很差。参数上,读者可以把 pi/6 改成任意步进,但要注意角度差异太小会导致星座距离变近、高SNR段出现错误平层。

编码过程就是每个用户按自己的比特选一个码字,然后按资源块叠加。

% scma_transmit.m % 输入: cb (KxMxJ), bitIdx (Jx1, 取值1..M) % 输出: y (Kx1), 叠加后的无噪声发射向量 function y = scma_transmit(cb, bitIdx) [K, ~, J] = size(cb); y = zeros(K, 1); for j = 1:J y = y + cb(:, bitIdx(j), j); end end

代码里 bitIdx 是每个用户选中的码字索引,由 log2(M)=2 个比特映射而来。发射端的映射我用的是直接二进制:bits 两列第一列是低位。实际项目中常见做法是把直接映射换成格雷映射,能让星座点相邻比特只差1位,从而降低高SNR下的误比特率,但演示链路用直接映射不影响MPA本身的功能验证。

2.3 信道与接收向量:AWGN怎么加,EbN0怎么进sigma2

信道模型先按最简单的平坦AWGN来:y_received = transmit_vector + noise。MPA译码需要两个关键量:接收向量和噪声方差 sigma2。噪声方差的来源是EbN0,而EbN0到sigma2的换算是整个仿真里最容易出错的一步。正确做法是先由码本统计出平均码字能量,再除以 log2(M) 得到每比特能量,最后除以线性EbN0。

% 主仿真片段: 一帧的发射与加噪 bits = randi([0 1], J, 2); % 每用户2比特 bitIdx = bits * [1; 2] + 1; % 00->1, 01->2, 10->3, 11->4 y = scma_transmit(cb, bitIdx); % 由码本统计平均码字能量 Es = mean(mean(sum(abs(cb).^2, 1), 2), 3); % 标量 Eb = Es / log2(M); % 每比特能量 sigma2 = Eb / (10^(EbN0dB / 10)); % 复噪声总方差 n = sqrt(sigma2/2) * (randn(K,1) + 1i*randn(K,1)); r = y + n;

这一段里 sigma2 是复数噪声的总方差,所以拆成实部虚部各占一半。Eb 的计算用 mean 对三个维度求平均,不依赖某个具体用户或码字,这样即便码本被整体缩放,Eb 也会跟着变化,信噪比定义始终一致。换到频率选择性信道时,常见做法是把 h 参数从标量扩展成 K×J 矩阵,接收向量变成逐资源块乘系数再叠加,后面的MPA代码需要同步改。

获得 r 之后,发射端这一侧就完整了。接下来要面对的是整个仿真的核心:到底怎么从叠加的信号里把6个用户的比特逐个抽出来,也就是MPA译码部分。

3. MPA译码:因子图上的消息迭代,从公式到Matlab

3.1 因子图与消息定义:MPA在log域里到底存了什么

MPA把多用户检测建模成因子图上的消息传递。因子图里有两类节点:变量节点对应6个用户,功能节点对应4个资源块。用户 j 与资源块 k 之间有一条边,当且仅当 smap(k,j)=1。迭代过程中每条边上会跑两类消息:从用户到资源块的 V2F 消息,以及从资源块到用户的 F2V 消息。

初学者最容易搞混的一点:这里的消息不是定义在某个资源符号上的,而是定义在码字索引 m 上的,每条消息都是一个 M×1 的向量。原因是SCMA里一个码字同时决定用户在多个资源上的符号,如果把消息按符号维度拆开,码字内部的关联就丢了。所以 V2F{k,j}(m) 表示用户 j 告诉资源块 k「我的码字是 m」的对数置信度,F2V{k,j}(m) 表示资源块 k 根据接收信号和其他用户的消息推断「用户 j 的码字是 m」的对数置信度。

所有消息都在对数域里运算。初始化时没有任何先验信息,每条边的 V2F 都置为全0向量,表示所有码字等概率。0是log域里的均匀分布,因为它对应概率1,而消息之间通常只关心相对大小,常数偏移在后面的更新中会自动抵消。这个初始化看似随意,却直接影响迭代是否稳定,第四章会专门讲。

MPA的迭代更新分两步。功能节点更新(F2V)的max-log近似形式是:

L(k→j)(m) = max over others的码字组合 [ -|y_k - Σ h_k,i x_k,i|^2 / sigma2 + Σ V2F消息 ]

变量节点更新(V2F)则是把该用户连接的其他资源块发来的 F2V 消息逐码字相加。

3.2 功能节点更新:把每个资源块上的16种组合算干净

以资源块1为例,它连接用户1、2、3。要给用户1更新 F2V,需要固定用户1的码字索引 m,然后遍历用户2和用户3的所有码字组合,总共4×4=16种。对每种组合计算一个度量:信道似然项加上来自用户2、用户3的 V2F 消息,取最大值作为用户1的消息值。dc=3 时组合数只有16,这就是SCMA选择稀疏结构的核心原因——如果用传统叠加编码,组合数会是 M^dc,资源块上用户再一多就彻底算不动了。

% scma_mpa_decode.m 核心片段: 功能节点更新 % cb: KxMxJ, smap: KxJ, y: Kx1, h: Kx1 平坦信道增益 % V2F/F2V: KxJ cell, 每个元素是 Mx1 log域消息 for k = 1:K nbrs = find(smap(k, :)); % 资源k上的用户集合, dc=3 for j = nbrs others = setdiff(nbrs, j); % 其余两个用户 fMsg = zeros(M, 1); for m = 1:M xj = cb(k, m, j); % 用户j在资源k上的符号 best = -inf; for idxA = 1:M xA = cb(k, idxA, others(1)); vA = V2F{k, others(1)}(idxA); for idxB = 1:M xB = cb(k, idxB, others(2)); vB = V2F{k, others(2)}(idxB); metric = -abs(y(k) - h(k)*(xj + xA + xB))^2 / sigma2; best = max(best, metric + vA + vB); end end fMsg(m) = best; end F2V{k, j} = fMsg; end end

这段代码用两层循环固定遍历 others(1) 和 others(2),是因为演示配置 dc=3。如果以后换参数让 dc 变大,两层循环要改成动态组合遍历,常见做法是用 ndgrid 生成所有组合索引矩阵再批量计算。metric 里的 h(k) 目前是标量平坦增益,换成频率选择性信道时要把 h 扩展成逐用户系数,令 xj = h(k,j)*cb(k,m,j)。vA 和 vB 来自上一次迭代的 V2F 消息,因此迭代顺序必须是先更新全部 F2V,再更新全部 V2F,不能同一个迭代次数里混着用新一轮和旧一轮的消息,否则会引入时序偏差。

3.3 变量节点更新与硬判决:迭代怎么收尾,比特怎么出

变量节点更新比功能节点简单得多。用户 j 发给资源块 k 的 V2F 消息,等于用户 j 从其他资源块收到的 F2V 消息逐码字相加。这个更新的依据是:用户 j 只连接两个资源块,在贝叶斯意义下后验概率正比于各分支似然的乘积,取对数后就是求和。

% 变量节点更新 for j = 1:J resNbrs = find(smap(:, j)); % 用户j的资源集合 for k = resNbrs vMsg = zeros(M, 1); for l = setdiff(resNbrs, k) vMsg = vMsg + F2V{l, j}; end V2F{k, j} = vMsg; end end

迭代结束后,每个用户把所有 F2V 消息累加得到后验对数消息,取最大值的索引作为码字估计,再映射回比特。

% 迭代结束后判决 LLR = zeros(M, J); for j = 1:J resNbrs = find(smap(:, j)); for k = resNbrs LLR(:, j) = LLR(:, j) + F2V{k, j}; end end [~, idxEst] = max(LLR, [], 1); bitsEst = zeros(J, 2); bitsEst(:, 1) = mod(idxEst - 1, 2); bitsEst(:, 2) = floor((idxEst - 1) / 2);

这里比特映射和发射端保持一致:发射端 bitIdx = bits(:,1) + 2*bits(:,2) + 1,译码端反过来拆就行。需要提醒的是,SCMA的因子图是有环的,MPA本质是在环图上做近似推断,所以迭代次数不是越多越好。常见做法是固定5到8次,既保证消息在环路里充分传播,又避免过度自信。迭代次数的坑在下一章详细展开。

4. MPA仿真避坑:五个让曲线翻车的细节

4.1 功率归一化和EbN0换算:曲线整体平移的元凶

坑1:EbN0被当成SNR直接用。现象是误码率曲线形状完全正常,但整体比预期差3dB左右,而且不同码本配置下偏移量还不一样。原因在于SNR是信号功率对噪声功率的比,而EbN0是每比特能量对噪声功率谱密度的比,中间差了码字能量和比特数的换算。如果你的仿真直接把 EbN0dB 塞进 sigma2 = 10^(-EbN0dB/10),那曲线一定会平移。解决方法是严格走一遍:先用码本统计平均码字能量 Es,再算 Eb = Es/log2(M),最后 sigma2 = Eb / EbN0_linear。注意这里 Es 要用 mean 对码本三个维度求平均,而不是拿某一个用户的码字能量代替。

坑2:码本符号没做能量归一化。现象是不同用户之间发射功率不一致,误码曲线在高SNR段出现地板,某个用户永远比别的用户差一截。原因是演示码本里如果直接拼接QPSK符号,旋转后符号幅度虽然不变,但码字总能量可能因为非零位置分布不同而产生差异,尤其当 df 不固定时更明显。解决方法是统计每个用户的平均码字能量后做归一化:cb = cb ./ sqrt(Es_per_user),并更新 Eb 的计算。整体缩放不会改变MPA译码的相对判决,只是让信噪比定义变公平。实际项目里还见过有人忘了把噪声方差跟信号一起缩放,导致曲线随码本缩放左右乱飘,这是同一个坑的两面。

4.2 迭代相关的坑:不收敛、震荡和max-log近似

坑3:迭代次数越多,误码率反而越差。现象是迭代次数从5加到20,BER先降后升,甚至在某个迭代次数出现周期性波动。原因有两个:一是SCMA因子图有环,多轮迭代会让消息在环路上自我强化,最终偏向某个局部最优;二是log域消息没有做归一化,绝对值随迭代不断增大,数值上也开始不稳定。解决方法是固定迭代次数在5到8次,同时在每次变量节点更新后对消息减去其最大值。减去最大值不改变码字间的相对关系,纯粹是为了防止数值溢出,这样消息值域就维持在0附近。

坑4:max-log近似在高SNR段出现假平层。现象是曲线低SNR段趋势正常,越往高SNR走误差越大,甚至出现一段斜率变缓的伪平层。原因是功能节点更新时用 max 替代了 log-sum-exp,这一近似在消息可靠性低时误差较小,但消息差异变大后,被忽略的 exp 项占比上升,判决边界被扭曲。解决方法是高SNR段改用完整LogSumExp,Matlab里可以写成 max(a) + log(sum(exp(a - max(a))))。如果嫌复杂度高,另一个做法是只在最后两次迭代用LogSumExp,前几次用max-log,性能损失很小。

4.3 索引与结构:稀疏模式对不齐这类玄学

坑5:发射端和译码端稀疏模式对不齐,BER永远在0.3以上。现象是无论SNR多高,误码率都降不下去,看起来像译码器完全失效。原因多半是编码侧和译码侧各写了一份 smap,其中一份转置了,或者用户顺序换了,因子图上的边连错对象。这类问题特别隐蔽,因为代码不报错,只是结果荒谬。解决方法是把 smap、码本、用户顺序全部打包进一个参数结构体,编码和译码共用同一份数据,不在任何函数里重新声明。血泪经验是:换参数文件或重排用户顺序时,先跑一帧无误码的端到端测试,把所有译码结果打出来人工核对,再开始大面积扫SNR。

5. 从跑通到可信:用收敛性检查和误码统计给仿真结果上保险

5.1 动态停止迭代:监视LLR变化,别数固定圈数

固定迭代次数5到8次虽然稳定,但不同SNR段需要的迭代数不一样。低SNR时消息收敛慢,高SNR时两三轮就够了。常见做法是监视相邻两次迭代的变量节点LLR变化量,变化小于阈值就提前退出,既能省时间又能避免过迭代。阈值我一般取1e-3,按所有用户所有码字索引的平均绝对差来算。

% 动态停止: 计算相邻迭代的LLR平均变化量 delta = mean(abs(LLR_new(:) - LLR_old(:))); if delta < 1e-3 break; end

这个阈值灵敏度很高,实测中1e-2到1e-4之间对BER影响不大,主要是省迭代。注意 LLR_new 要在每轮迭代结束时先算出来,和上一轮的 LLR_old 比较后再进入下一轮。

5.2 误码统计至少攒够多少:关于BER曲线最后一段的可信度

误码率曲线最后一段往往只有几十个误码,波动极大。按置信区间的说法,二项分布标准差是 sqrt(p(1-p)/n),想置信带宽在20%以内,至少得攒够100个误码。也就是说,如果某一帧SNR段只有10个误码,要继续跑到至少100个误码再记下BER,否则曲线最后一段就是在画噪声。

我现在的习惯是:每次仿真把随机种子、码本文件名、迭代次数、信噪比扫描范围一起存进结果文件,方便事后定位问题。曾有一次调了个参数曲线明显变好,后来发现是漏了噪声归一化,从此所有仿真都先跑一发端到端无误码测试再信参数。收敛检查、误码计数、参数存档这三件事做到位,曲线才是能给别人看的东西。希望帮到你。

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进