多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

SCMA仿真核心:MPA译码完整链路与Matlab实现

SCMA仿真核心:MPA译码完整链路与Matlab实现 简介基于MPA的SCMA的Matlab仿真压缩包是一份面向第五代移动通信系统的非正交多址接入技术仿真代码适合通信专业高年级学生、科研人员及无线通信算法工程师在学习和项目中参考。包内共11个文件其中9个M文件源码完整覆盖稀疏编码设计、星座映射、多用户叠加传输、信道加噪以及MPA迭代译码等关键环节另附Markdown说明和TXT文档介绍运行环境与参数设置整个压缩包仅11KB紧凑且易读。目前已有1674人下载学习。代码除标准MPA译码实现外还提供了多种改进型MPA变体与仿真主程序支持BPSK、QPSK、16QAM等调制方式切换并可调整用户数、信噪比、迭代次数等参数直接观察误码率及收敛特性变化从而深入理解SCMA系统在不同调制与负载条件下的性能差异也为后续优化编码矩阵、改进译码策略提供了可直接复用的实验平台。1. 为什么SCMA仿真绕不开MPA一份Matlab仿真zip背后的完整链路经常有人拿到“基于MPA的SCMA的Matlab仿真”这类zip包解压后看到一堆码本矩阵、mpa函数和嵌套for循环不知道先看哪个文件也不知道BER曲线画出来该跟谁对比。其实这条链路特别短发射端把每个用户的比特映射成稀疏码字在资源块上故意重叠接收端用MPA消息传递算法把叠加信号分离出来。整套仿真不需要真实信道环境Matlab脚本加上AWGN就能把误码率曲线跑出来目标是用数据说清楚“SCMA在高过载下为什么比OFDMA更省频谱”。这篇笔记适合正在做5G多址接入课程设计、需要复现SCMA仿真结果、或者想真正弄懂MPA译码流程的人。下面按“码本与因子图 → MPA译码公式与代码 → 完整仿真链路 → 踩坑记录 → 码本微调”的顺序讲代码可以直接抄进脚本跑。2. 先把SCMA链路拆开码本、因子图与6用户4资源块的选型理由2.1 SCMA编码的本质比特到高维稀疏星座点的映射SCMA和传统CDMA最大的区别是CDMA用扩频序列区分用户SCMA用“稀疏码字”区分。每个用户拿一本码本M个候选码字每个码字是一个K维复数列向量但向量里只有少数几个非零元素。发射时用户从自己码本里挑一个码字多个用户的码字在K个资源块上叠加同一个资源块上会有好几个用户信号挤在一起。最标准的演示配置是6用户4资源块J6个用户共享K4个资源块每用户占2个资源块行重2每资源块承载3个用户列重3码本大小M4。过载率是J/K1.5意味着比传统正交多址多承载50%的用户。接收端要做的事就是把这3个用户在一个资源块上的叠加信号拆开。码本生成的代码可以这样写注意每列能量必须归一化% 生成6用户4资源块SCMA码本 % 因子图矩阵6行(用户) x 4列(资源块) F [1 1 0 0; 1 0 1 0; 1 0 0 1; 0 1 1 0; 0 1 0 1; 0 0 1 1]; M 4; % 码本大小每用户4个候选码字 K 4; % 资源块数 J 6; % 用户数 base [1, -1, 1i, -1i]; % 基础星座点QPSK codebook cell(J, 1); for u 1:J res_idx find(F(u, :)); % 该用户占用的资源块位置 codewords zeros(K, M); for r 1:length(res_idx) % 不同资源块用不同相位旋转避免两个维度星座完全重叠 theta (u-1) * 2*pi/J (r-1) * pi/M; codewords(res_idx(r), :) base * exp(1i*theta); end % 每个码字能量归一化到1这是BER曲线横轴不漂移的前提 codewords codewords ./ sqrt(sum(abs(codewords).^2, 1)); codebook{u} codewords; end逻辑说明F矩阵每行两个1表示每个用户占用两个资源块这是SCMA获得频率分集的关键每列三个1表示每个资源块上正好叠加三个用户MPA译码时的状态数就是4的三次方等于64种。相位旋转的伪随机设计让不同用户的星座在同一个资源块上错开降低叠加后星座点碰撞的概率是码本设计里最简单有效的一招。参数说明M4对应每个码字携带log2(4)2比特一个SCMA符号在6用户下总共有12比特。r循环里的( r-1)*pi/M让同一用户在不同资源块的星座旋转不同角度避免两个非零维度发射完全相同的星座。归一化那行是按列求能量再除漏掉这一行后面所有SNR换算都会错。参数取值含义M4每用户4个候选码字每个码字承载2比特K4时频资源分成4个资源块J66个用户共享4个资源块过载率1.5行重2每用户占2个资源块提供频率分集列重3每资源块叠加3个用户MPA状态数4³642.2 为什么要用MPA最优联合检测的复杂度被困死了接收端理论上可以走最大似然检测把6个用户的所有码字组合全部枚举算概率选最大。M4时全局组合数是4的6次方等于4096种一个符号枚举一次还能忍仿真几千个符号就非常吃力扩展到12用户8资源块、M16时组合数直接变成16的12次方任何暴力枚举都不可行。MPA做的事就是利用码字稀疏性把全局枚举拆成每个资源块上的局部枚举。对上述6用户4资源块配置每个资源块只有3个用户叠加状态数4³644个资源块加起来每轮迭代算4×64256个状态和全局4096相比复杂度降了一个数量级。代价是需要迭代若干轮让消息在因子图上多传播几次才能收敛。这就是标题里“基于MPA的SCMA仿真”的核心套路MPA不是对联合检测的粗糙近似而是把联合检测分解成因子图上的局部计算再用消息传递把局部后验汇总回每个用户。因子图节点关系由F矩阵决定F(u,k)1表示用户u占用资源块k。后面译码函数里的更新顺序是“先资源块节点后变量节点”每一轮把某个资源块算出的置信度传给相邻变量节点变量节点汇总后再传回其他资源块如此反复。在过载率只有1.5的配置下MPA收敛后的性能和最大似然基本贴在一起差距通常在0.1dB以内。选这个配置做仿真的好处是既能看出MPA相对暴力枚举的复杂度优势又不至于把MPA本身的迭代误差混进对码本设计的判断里。新手做SCMA仿真第一站几乎都是这个标准配置和你手里那份zip里的主体结构一致。3. MPA译码从公式到Matlab消息传递、Log域与迭代实现3.1 消息传递的三条线初始化、资源节点更新、变量节点更新MPA的核心只有三个公式。初始化时没有先验信息每个用户取任意码字的概率均匀分布Log域里就是全0消息。资源节点更新第k个资源块的接收值y_k已知时要计算“用户u取了第m个码字、其他用户取各种组合”的联合似然公式是对资源块k上叠加的用户集合N(k)每种码字组合都算一个欧氏距离度量再把其他用户传来的先验消息乘进去最后把同一个用户同一个码字的所有组合累加。对6用户4资源块配置N(k)里有3个用户每种组合的枚举量是4³64。变量节点更新用户u把自己的信息传给资源块k时要收集除k以外所有相邻资源块传来的消息取乘积。因为每个用户只占2个资源块这个更新简化为“把另一个资源块传来的整条消息搬过来”。正是这一步让消息绕过因子图的环迭代才有意义。最后一步是符号判决。每个用户的总后验等于它两个资源块消息的乘积转成比特LLR后统计误码。这个流程里最容易出问题的就是Log域实现时消息没有清零或归一化下面这段代码是完整可运行的实现。3.2 Log域MPA的最小可运行代码直接按概率写会因为连乘下溢成NaN工程上都转Log域。概率相乘变成指数相加欧氏距离项变成-|y-hx|²/(2σ²)省掉大量exp计算。下面函数输入接收向量、噪声方差、码本和因子图输出每个用户的比特级软信息function LLR mpadecode_logdomain(y, sigma2, codebook, F, M, J, K, maxIter) % y : K*1 接收向量一个SCMA符号 % sigma2 : 复噪声每维方差 % codebook : J*1 cell每个元素是 K*M 复数矩阵 % F : J*K 因子图矩阵 % LLR : J*log2(M) 软比特输出 logM log2(M); % 预计算每个资源块上的活跃用户列表和组合总数 users_res cell(K, 1); num_comb zeros(K, 1); for k 1:K users_res{k} find(F(:, k)); num_comb(k) M^(length(users_res{k})); end % Log域初始消息置0等价于均匀先验 MsgV2F zeros(J, M, K); % 变量节点到资源节点 MsgF2V zeros(J, M, K); % 资源节点回变量节点 for iter 1:maxIter % 每一轮迭代开始必须清空F2V重新累加否则消息会跨轮累积出错 MsgF2V zeros(J, M, K); % ---- 资源节点更新 ---- for k 1:K us users_res{k}; dr length(us); % 该资源块上叠加的用户数这里是3 for s 0:num_comb(k)-1 % 把s拆成dr个M进制位得到一组码字索引组合 tmp s; state zeros(1, dr); for i dr:-1:1 state(i) mod(tmp, M) 1; tmp floor(tmp / M); end % 计算该组合的欧氏距离度量 对数域先验 tx_sum 0; prior_sum 0; for i 1:dr u us(i); tx_sum tx_sum codebook{u}(k, state(i)); prior_sum prior_sum MsgV2F(u, state(i), k); end metric -abs(y(k) - tx_sum)^2 / (2*sigma2) prior_sum; % 把这条路径的贡献累加到对应用户的F2V消息上 for i 1:dr u us(i); MsgF2V(u, state(i), k) logsumexp(MsgF2V(u, state(i), k), metric); end end % 归一化防止消息无限增长或溢出 for i 1:dr u us(i); MsgF2V(u, :, k) MsgF2V(u, :, k) - max(MsgF2V(u, :, k)); end end % ---- 变量节点更新 ---- for u 1:J res_list find(F(u, :)); for k res_list acc zeros(1, M); for g res_list if g ~ k acc acc MsgF2V(u, :, g); end end MsgV2F(u, :, k) acc; end end end % ---- 软信息输出 ---- LLR zeros(J, logM); for u 1:J post zeros(1, M); for k find(F(u, :)) post post MsgF2V(u, :, k); end post post - max(post); % 防exp溢出 prob_m exp(post); prob_m prob_m / sum(prob_m); % 归一化成概率 % 码字索引到比特的映射用格雷码按位算LLR for b 1:logM mask bitget(0:M-1, b); LLR(u, b) log(sum(prob_m(mask 1)) 1e-12) - ... log(sum(prob_m(mask 0)) 1e-12); end end end function v logsumexp(a, b) m max(a, b); v m log(exp(a - m) exp(b - m)); end逻辑说明资源节点更新里s从0遍历到63每次拆成三个M进制位就得到3个用户的码字索引组合metric中的-abs(y-tx_sum)²/(2σ²)是高斯信道下的对数似然prior_sum把上一轮变量节点传来的消息加进来这就是“消息传递”的实体现。MsgF2V每轮迭代开头必须清零这是最容易写漏的一行漏掉后消息会跨轮累积BER曲线直接失效。参数说明maxIter在6用户4资源块配置下建议12到20轮。sigma2是复噪声的每维方差定义成sigma2N0/2时距离项分母2σ²正好等于N0这是最容易写错的地方。归一化那三行不能删删掉后Log域数值会逐渐漂移最终在某个SNR点爆出NaN。想提速可以把logsumexp替换成max那就是Max-Log近似速度提升明显代价是BER曲线整体偏移0.2到0.3dB。注意如果你手里那份zip里的mpa函数没有“每轮清零F2V”这一步先把它加上。这是我在整理常见实现时发现最高频的代码缺陷往往表现为“迭代次数增加但BER不下降”。3.3 迭代参数怎么设12次还是20次迭代次数是MPA仿真里最微妙的参数。6用户4资源块的因子图有环消息需要多轮传播才能稳定。经验值是12轮起步20轮可以认为是收敛参照。判断方法很简单把maxIter分别设成8、12、20、30各跑一次同一个SNR点的BER如果12次和20次差距小于0.05dB就说明12次已经够用。用Max-Log近似时收敛会慢一些建议直接取20。如果BER曲线在高SNR段出现“平台期”先别怀疑码本极大概率是迭代次数不够。把迭代次数翻倍看曲线是否继续下降是排查这个问题的标准操作。4. 完整链路从比特生成到BER曲线的Matlab仿真搭建4.1 发射端码本索引生成与稀疏叠加发射端按用户独立生成随机比特每log2(M)2个比特映射成一个码字索引。码字索引到比特的映射要固定同一张表接收端LLR计算和它对称。然后每个用户在自己的两个非零资源块上放置码字其余位置为0最后所有用户叠加成K维发送向量% 发射端生成一个SCMA符号并完成稀疏叠加 M 4; logM 2; J 6; K 4; % F和codebook沿用第2节变量 bits_u randi([0 1], J, logM); % 每用户2比特 idx_u bi2de(bits_u, left-msb) 1; % 码字索引范围1~M tx_sym zeros(K, J); for u 1:J tx_sym(:, u) codebook{u}(:, idx_u(u)); % 取出该用户的码字 end tx sum(tx_sym, 2); % 稀疏叠加K*1发送向量逻辑说明bi2de把每行的2比特转成0到3的整数加1对齐Matlab列索引。tx_sym每一列是用户u的码字只有两个非零位置按行求和后每个资源块上正好是3个用户码字对应位置的叠加。发送向量每个元素的期望功率等于列重3因为每个资源块上有3个能量为1的码字分量。参数说明bi2de的left-msb要和接收端mask的bitget顺序一致。写完F矩阵后建议打印sum(F,2)和sum(F,1)核对行重2、列重3这个检查三十秒能省掉后面排查半天。4.2 接收端AWGN信道与MPA译码调用AWGN下接收向量ytxnoise噪声按EbN0换算。关键是噪声方差设定很多仿真曲线横移都是栽在这里。换算关系是一个SCMA符号总比特数为J·logM12资源块数K4频谱效率R12/43 bit/s/Hz。每资源块期望信号能量Es等于列重3。% 接收端加噪 MPA译码 误码统计 EbN0_dB 6; R J * logM / K; % 频谱效率 3 Es 3; % 每资源块平均信号能量等于列重 sigma2 Es / (2 * R * 10^(EbN0_dB/10)); % 复噪声每维方差 noise sqrt(sigma2) * (randn(K, 1) 1i*randn(K, 1)); y tx noise; LLR mpadecode_logdomain(y, sigma2, codebook, F, M, J, K, 20); % 硬判决并与发射比特对比 bit_est double(LLR 0); ber_seg mean(bit_est(:) ~ bits_u(:));逻辑说明噪声生成用randn实部和虚部各乘sqrt(sigma2)得到每维方差都是sigma2的复高斯噪声。LLR0表示该比特为1的后验概率更大直接把软信息硬判决成比特跟发射端比对。单个符号的BER没有统计意义下面要循环几千帧再平均。参数说明sigma2公式里的Es用3而不是1。分子取3的原因是一个资源块上有3个能量为1的码字叠加实际接收能量就是3。用错会整条BER曲线平移约4.8dB且斜率看起来完全正常这是最隐蔽的横轴错误。提示跑全曲线前先固定一个低SNR点用BPSK理论BER做参照校准sigma2比直接跑完再怀疑横轴靠谱得多。4.3 跑完看什么BER曲线与理论界的对比主脚本套一层EbN0循环从0dB扫到14dB每个点跑5000到10000个SCMA符号最后semilogy画图EbN0_list 0:2:14; BER zeros(size(EbN0_list)); nSymbols 5000; for i 1:length(EbN0_list) bit_err 0; bit_tot 0; EbN0_dB EbN0_list(i); R J * logM / K; Es 3; sigma2 Es / (2 * R * 10^(EbN0_dB/10)); for s 1:nSymbols % 生成比特 - 发射 - 加噪 - MPA译码同前面逻辑略 bit_err bit_err sum(bit_est(:) ~ bits_u(:)); bit_tot bit_tot J * logM; end BER(i) bit_err / bit_tot; end semilogy(EbN0_list, BER, o-); grid on; xlabel(EbN0 (dB)); ylabel(BER);跑完看两个参照低信噪比下BER应该接近高斯信道下单用户QPSK的性能因为叠加干扰在低SNR下被噪声主导高信噪比下曲线斜率反映码本的最小欧氏距离特性斜率越陡码本越好。如果曲线在10dB以上出现地板效应先查迭代次数和码本归一化不要急着怀疑MPA本身。高SNR点BER抖动是正常的把该点的nSymbols从5000提到20000能压住抖动。想要更快出图低SNR降到2000帧也够用。5. SCMAMPA仿真避坑6条血泪经验5.1 现象BER曲线整体向右平移约4.8dB高信噪比贴着单用户QPSK下界原因sigma2公式里的Es取错。按每个用户码字能量1去算实际每个资源块上叠加3个用户总信号能量是3而不是1。整条曲线平移了10log10(3)4.77dB但曲线斜率看着完全正常特别容易漏掉。解决永远从实测能量反推。发射端先发几千个随机符号统计mean(abs(tx).^2)打印出来用这个实测值代回sigma2公式。这个方法同样能校验码本能量没归一化的情况一石二鸟。5.2 现象Log域MPA第一次跑就出现NaNBER曲线在0.5附近横着原因消息归一化没做exp之后数值爆炸或者logsumexp实现有误两个负无穷相加变成-Inf。另一个高频原因是发射端的bi2de和接收端LLR的bitget位顺序不一致导致LLR符号整体反了。解决资源节点更新完立即对当前节点消息减max输出LLR前再减一次。调试时把中间消息打印出来检查是否每轮迭代后落在[-30, 0]区间。位顺序问题就固定随机种子rng(42)逐比特比对发射端映射表一次能定位。5.3 现象BER曲线在中间出现平台期迭代8次和20次结果一样原因迭代次数太少消息尚未在因子图上充分传播。6用户4资源块配置有环信息绕一圈需要好几轮。Max-Log近似会进一步放慢收敛有时12次和20次仍有0.1dB的差距。解决把maxIter先设成30跑一次拿30次的结果当基准再分别用8、12、20次画出对比曲线。选迭代次数时找“再翻倍性能变化小于0.05dB”的最小值。高SNR段一旦出现“假地板”优先怀疑这里。5.4 现象MPA性能比最大似然差1dB以上怀疑码本设计有问题原因大部分情况是消息更新里漏乘了信道系数。AWGN仿真里信道幅度虽然恒为1但实现时不保留h_{kv}这一项以后换平坦衰落信道就会翻车。另一种常见原因是Max-Log近似后消息没重新归一化导致软信息幅度整体偏大或偏小。解决在资源节点更新里显式乘上信道系数h默认设成全1矩阵。验证时写一个4^64096状态枚举的ML解调函数在同一个噪声序列上比较两者的LLR方向一致率。如果MPA的LLR大量反转去查变量节点更新是否把两个邻居资源块的消息都乘上了。5.5 现象仿真跑得很慢一个SNR点要等十几分钟原因嵌套循环太多而且一些本可预计算的量在EbN0循环里反复生成。6用户4资源块每资源块64种组合4个资源块每轮256次度量计算如果每个符号都重新生成组合索引表和码本时间全耗在重复劳动上。解决码本、每资源块的用户列表、组合索引表全部提到EbN0循环外预计算。资源块枚举用位拆解方式一次生成组合表缓存不要在迭代里重建。再进一步把logsumexp替换成max做Max-Log近似速度再快约一半。先保证正确再优化避免边优化边引入NaN最后分不清是优化错了还是译码出错了。5.6 现象matlab中文注释乱码代码读起来像黑匣子原因老脚本是GBK编码保存新版matlab如matlab 2023b默认按UTF-8打开中文注释全乱。这是编码历史问题不是代码坏了。解决编辑器中直接另存为UTF-8编码。整个zip里多个脚本都乱码时用记事本打开确认原编码再统一批量转换。MPA的状态更新逻辑本来就绕中文注释能省很多理解成本别让乱码影响判断。6. 码本微调与验证技巧把BER曲线再往下压一点基础链路跑通后最有价值的微调方向是码本的相位旋转。第2节里theta用了均匀分布但这个均匀分布远不是最优解。做法是把相位因子改成待搜索变量以“任意多用户码字组合叠加后的最小欧氏距离最大”为目标在[-π, π]上扫一组theta。搜索在MPA之前单独做把4个资源块上所有可能叠加的星座组合都枚举一遍计算最小距离选最大的一组固定进码本。这一步能直观看到BER曲线在10⁻³量级附近明显下移。另一个必须养成的验证习惯是做ML对照。虽然6用户4资源块的ML状态数有4096单帧枚举不快但只在高SNR点跑2000个符号做对照完全可行。同一个噪声序列上ML的BER和MPA的BER之间的差距就是MPA的近似损失。如果这个损失小于0.2dB说明你后续所有码本改进的结论都建立在可信的译码器之上。我后来把配置扩到12用户8资源块时每次改码本都先跑小规模ML对照再跑完整MPA这个习惯帮我排掉了至少三个“码本好像变好了、其实是译码器没收敛”的假阳性结论。最后建议把参数集中管理。把F矩阵、M、J、K、maxIter、sigma2计算全部写进一个config脚本里发射端和接收端都从它读取。我最早做这个仿真时参数散落在两三个脚本里每次调参都怀疑自己改到了别处后来集中到一个文件才踏实。如果你从zip包里拿到的版本比较老按这个方式重构一遍会比在旧脚本上打补丁省心得多。希望帮到你。本文还有配套的精品资源点击获取
返回列表