简介一套面向SCMA稀疏码分多址系统研究者的PM-MPA检测算法MATLAB源码包适合5G通信、非正交多址方向的学生与工程师快速上手算法仿真。资源共5个文件以m脚本为主4个m文件1个zip核心代码PM_MPA.m实现了产品矩阵消息传递迭代逻辑simulation.m可配置用户数、星座与信道参数scmaenc.m完成稀疏码字映射log_sum_exp.m用于数值稳定的对数求和另附带SCMA-SD-MPA的zip压缩包便于对比不同检测器性能。整个压缩包仅7KB结构简洁无冗余文件。已有611人学习下载。借助这套程序读者可深入理解消息传递机制在瑞利信道下的检测流程并直接进行性能评估与算法改进实验是研究SCMA低复杂度检测算法值得参考的实用工具。1. 同样是 SCMA 仿真为何你的 MPA 慢得让人等不下去我碰到过不少做物理层仿真的朋友把 SCMA 的检测器从 QPSK 的 MAP 检测换成 MPA 之后第一反应是“终于跑通了”第二反应是“这跑得也太慢了”。PM-MPAPartially Marginalized Message Passing Algorithm部分边缘化消息传递算法不是一套新的编解码方案而是对标准 MPA 的检测循环做了一次工程修正把迭代中最先收敛的一部分变量节点从消息更新里“摘”出去让剩余的迭代只处理真正还在变化的消息。对 6 用户 4 资源、母矩阵度数为 3 的经典 SCMA 配置PM-MPA 通常能砍掉三成左右的计算量BER 损失可以控制在 0.2 dB 以内。这篇文章就用 MATLAB 把 PM-MPA 完整走一遍从因子图搭建到边缘化集合选择再到那些让我曾经翻车的数值坑。2. PM-MPA 在加速什么因子图上的消息传递与部分边缘化2.1 标准 MPA 的消息传递机制与复杂度来源SCMA 的检测问题本质上是一个多用户联合估计问题。J 个用户各自从 M 个码字里选一个叠加在 K 个资源上传输接收端拿到的是 y_k Σ h_kj·c_j(m_j)_k n_k其中 c_j(m_j)_k 是第 j 个用户第 m_j 个码字在资源 k 上的分量。穷举所有用户的码字组合 M^J在 J6、M4 时是 4096 种直接做 MAP 在演示场景还能忍一旦接入信道编码或提高调制阶数就完全不可行。MPA 把联合检测拆到因子图上。因子图有 J 个变量节点对应 J 个用户的码字选择和 K 个功能节点对应 K 个接收资源消息在两类节点之间来回传递。功能节点向变量节点传的消息反映的是“结合本资源上其他用户的分布当前用户选某个码字的可信度”计算时需要枚举这个功能节点连接的所有 d_f 个用户的码字组合变量节点向功能节点传的消息则把来自其他功能节点的信息乘起来再与自己的先验相乘。这段乘积在实现时通常转成对数域求和避免码字概率越乘越小导致数值下溢。复杂度的大头在功能节点更新。每个功能节点连接的 d_f 个用户每个用户有 M 个码字所以要遍历 M^d_f 个组合。对 S(4,6) 母矩阵K4、d_f3、M4那就是每一轮迭代 4×64256 次指数项计算。如果跑 6 轮迭代就是 1536 次。仿真 100 万帧时这个量级在 MATLAB 里用循环实现压测耗时非常难看。这也是为什么 PM-MPA 这类“降低单轮迭代开销”的策略值得专门做一次代码级拆解。2.2 部分边缘化不是近似是主动剪枝PM-MPA 的出发点非常朴素迭代中消息收敛的速度并不均匀。信道条件好的用户或者因子图上度数低的用户往往一两轮迭代后后验概率就稳定了。此时继续让这些变量节点参与每一轮的消息更新等于反复计算一堆已经不再变化的量。部分边缘化做的就是把这类节点“冻结”住。通常分三步走先做少量预热迭代让所有消息在因子图上传播开对每个变量节点计算当前的后验边缘概率按一定准则挑出 N_marg 个“已基本确定”的用户后续迭代中这些用户的变量节点消息不再更新功能节点更新时也不再枚举它们的所有码字只使用固定的后验或者直接用后验最大的码字参与组合计算。严格讲这只是一种降低迭代计算量的策略而不是对 MPA 消息公式本身的近似。被选中的节点后验已经收敛硬固定不改变因子图上联合后验分布的形状实际引入的性能损失来自“被固定节点的真实码字判断错误”这种情况。因此边缘化集合的选取直接决定了 BER 损失大小这也是后面第 4 章要单独展开的原因。2.3 边缘化前后复杂度对比从 4^3 到 4^d_f用经典配置算一笔账。母矩阵 F 如下每行 3 个 1每列 2 个 1用户 1 占资源 1、2用户 2 占资源 1、3其余用户以此类推。标准 MPA 下功能节点更新的组合遍历总数是 K × M^df 4×64 256 次指数项计算每轮。PM-MPA 选了不同数量的用户边缘化之后剩余组合数取决于被冻结节点在因子图上的分布配置每轮 FN 组合总数6 轮迭代总指数计算量相对标准 MPA标准 MPAdf325615361.0PM-MPA边缘化 1 个用户 116011520.75PM-MPA边缘化 2 个用户 1 和用户 6647680.50这里假设预热 2 轮按标准 MPA 计算正式迭代 4 轮使用剪枝后的组合遍历。边缘化 1 个用户时资源 2 和资源 3 各少一个自由维度变成 M^216边缘化 2 个用户时如果选用户 1 和用户 6四个资源各少一个自由维度全部变成 16每轮组合总数就是 4×1664。这解释了为什么边缘化集合的分布比“选几个”更重要选得均匀复杂度收益才明显。3. 用 MATLAB 搭 PM-MPA 最小检测器数据结构与核心循环3.1 系统参数与因子图矩阵MATLAB 里搭 SCMA 系统我习惯把参数集中在一个脚本头部后面无论是发端、收端还是仿真循环都直接引用这些量。下面是一份最小可用的参数配置% scma_config.m K 4; % 资源节点数 J 6; % 用户数 M 4; % 码本大小每用户码字个数 niter 6; % 标准 MPA 迭代轮数 Warmup 2; % PM-MPA 预热迭代轮数 SNRdB 10; % 仿真信噪比 % 因子图母矩阵 S(4,6)行资源列用户 F logical([1 1 1 0 0 0; 1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1]);这里 F 的行重都是 3列重都是 2对应每个资源上叠加 3 个用户、每个用户占用 2 个资源。码本我临时生成一份简化版从 QPSK 星座出发给每个用户一个不同的相位旋转再把非零分量放到该用户占用的资源位置上qpsk [11i, 1-1i, -11i, -1-1i] / sqrt(2); CB cell(J,1); for j 1:J cb zeros(M, K); nz find(F(:,j)); % 该用户占用的资源 cb(:, nz) qpsk * exp(1i * (j-1) * pi/8); CB{j} cb; end这份码本只保证稀疏结构与相位可分性距离标准 SCMA 码本还有差距。做正式误码率实验时建议换成文献中的低互相关 SCMA 码本否则某些用户对会被相同相位旋转的星座点拉近边缘化时更容易选错。发射侧把码字按因子图叠加加高斯白噪声H ones(K, J); % 先做单径理想信道后续再换成平坦衰落 bits randi([0 1], J, log2(M)); mapIdx bi2de(bits) 1; x zeros(K, 1); for j 1:J x x H(:,j) .* CB{j}(mapIdx(j), :).; end y x sqrt(10^(-SNRdB/10)/2) * (randn(K,1) 1i*randn(K,1));这里的 H 先用全 1 矩阵相当于先把算法逻辑调通再上衰落信道。SNR 转到噪声方差的公式是 N0 10^(-SNRdB/10)符号能量归一化到 1所以实部虚部各加 sqrt(N0/2) 的高斯噪声。3.2 对数域 MPA 的实现骨架消息传递最容易踩的坑是概率域下溢。几个用户的消息乘几轮之后最小的概率项可能变成 0取对数就是 -Inf。所以建议直接用对数域消息更新全部用加法合并消息用 log-sum-exp 而不是直接乘。下面是一份能跑通标准 MPA 的骨架函数适用于列重为 2、行重为 3 的 S(4,6) 母矩阵function [bits, llr] mpa_detector(y, H, CB, F, niter, N0) % 输入: y(Kx1), H(KxJ), CB(Jx1 cell, 每用户 MxK), F(KxJ logical), N0 K size(y, 1); J size(F, 2); M size(CB{1}, 1); % 对数域消息初始化, 用 -1e8 而不是 -inf, 防止 logsumexp 出 NaN v2f -log(M) * ones(M, J, K); f2v -1e8 * ones(M, J, K); % 预生成 df3 时的组合矩阵 [A, B, C] ndgrid(1:M, 1:M, 1:M); combos [A(:), B(:), C(:)]; % M^3 x 3 for iter 1:niter % ---- 功能节点更新 ---- for k 1:K u find(F(:,k)); % 本资源关联的 3 个用户 for c 1:size(combos,1) ms combos(c,:); sd 0; for d 1:3 sd sd H(k,u(d)) * CB{u(d)}(ms(d), k); end e -abs(y(k)-sd)^2 / N0; % 该组合对每个用户的 f2v 消息做 log-sum-exp 累加 for d 1:3 others setdiff(1:3, d); acc e v2f(ms(others(1)), u(others(1)), k) ... v2f(ms(others(2)), u(others(2)), k); f2v(ms(d), u(d), k) logsumexp([f2v(ms(d), u(d), k), acc]); end end end % ---- 变量节点更新 ---- for j 1:J r find(F(j,:)); % 该用户占用的 2 个资源 for k r others setdiff(r, k); % 发往资源k, 取来自其他资源的f2v之和 v2f(:,j,k) sum(f2v(:,j,others), 2); end end end % 后验边缘概率与 LLR llr zeros(J, log2(M)); for j 1:J r find(F(j,:)); logPost sum(f2v(:,j,r), 2); % 均匀先验, 比值时被约掉 logPost logPost - max(logPost); post exp(logPost); post post / sum(post); for b 1:log2(M) bit0 sum(post(mod(0:M-1, 2)0 1)); % 示例映射, 按实际码本映射表替换 bit1 1 - bit0; llr(j,b) log(bit1/bit0 eps); end end bits (llr 0); end代码里的两个关键点要说明一下。第一logsumexp 是 MATLAB 没有内置的函数我一般自己写取出最大值 maxVal返回 maxVal log(sum(exp(x - maxVal)))。初值用 -1e8 而不是 -Inf这样 logsumexp 不会在全 -Inf 的行上返回 NaN。第二变量节点更新时每个用户只占 2 个资源所以发往资源 k 的消息就是从另一个资源来的 f2v 消息本身如果换成列重不是 2 的母矩阵这段要改成通用累加形式。3.3 加入部分边缘化的最小改动标准 MPA 跑通后PM-MPA 的最小改动不是重写循环而是在迭代流程里加一个“冻结掩码”。先做预热迭代并计算后验确定哪些用户进入边缘化集合随后在正式迭代中按掩码跳过这些用户的变量节点更新同时让功能节点更新时使用固定码字% PM-MPA 控制流程 margSet [2 5]; % 选择边缘化的用户索引 margMask false(J, 1); margMask(margSet) true; % 预热先跑 Warmup 轮标准 MPA for iter 1:Warmup [v2f, f2v] mpa_one_iter(v2f, f2v, y, H, CB, F, N0); end % 计算候选用户的后验, 确定冻结消息 for j margSet r find(F(j,:)); logPost sum(f2v(:,j,r), 2); logPost logPost - max(logPost); fixedPost(:,j) exp(logPost); % 保存, 后续不再更新 [~, fixedIdx(j)] max(fixedPost(:,j)); % 硬判决码字索引 end % 正式迭代 for iter Warmup1:niter for k 1:K u find(F(:,k)); margInRes margMask(u); % 本资源是否连接冻结用户 if any(margInRes) % 组合枚举时, 冻结用户的码字直接取 fixedIdx % 示例: df3 且第3个用户被冻结时, 只剩两层循环 for m1 1:M for m2 1:M ms [m1 m2 fixedIdx(u(3))]; % 后续指数项计算与标准版本一致 end end else % 原始标准 MPA 更新 end end end这段是控制流骨架省略了“计算指数项并累加 f2v 消息”的重复代码因为那部分与标准版本逐字相同。重点在于PM-MPA 的工程侵入性很低主体检测器仍然复用标准 MPA 的每一轮更新逻辑边缘化只是一个外层的掩码与分支控制。这个设计也解释了为什么很多开源仿真里 PM-MPA 只比标准 MPA 多几十行代码。3.4 代码运行结果确认LLR 与 BER 的初步 sanity check第一次跑通别急着上 BER 曲线先做三件事确认代码没有方向性错误。第一把 SNR 设成 30 dB检测器的 LLR 符号应该完全正确后验最大值对应的码字索引应该与发射码字索引一致第二把 y 的实部虚部互换或把星座旋转 90 度观察 LLR 是否成比例变化排除索引错位第三固定随机种子对比标准 MPA 和 PM-MPA 在相同 2000 帧上的 BER如果两者差异小于 0.1 dB说明边缘化没有破坏检测能力。我一般把第三个对比脚本固定成目录下的独立文件每次换码本、换信道都先跑它再决定要不要继续调参。4. 决定 PM-MPA 性能的 5 个关键参数迭代次数、边缘化集合与外环调度4.1 内环迭代次数看收敛曲线在哪变平先看一个实验习惯不确定迭代设多少时不看 BER看对数域消息的变化量。第 t 轮和第 t-1 轮之间所有 v2f 消息的绝对值变化均值如果小于 1e-3多数情况下再往下迭代已经很难改变硬判决结果。标准 MPA 在 df3、M4 的 SCMA 上迭代 4 到 6 轮后 BER 曲线基本贴着 MAP 走。PM-MPA 因为冻结了一部分节点剩余还在更新的消息变化空间变小通常 3 到 4 轮足够。经验值是预热 2 轮加正式 3 轮在边缘化集合选得准的前提下就能和标准 MPA 6 轮持平。如果正式迭代只给 1 到 2 轮冻结节点的后验可能还没真正稳定BER 损失会明显放大。这里要留意的是MPA 每轮迭代是 FN 更新和 VN 更新交替一次我统计“轮数”时按完整的 FNVN 算不单独算 FN 半轮。4.2 边缘化集合怎么选度数、信道还是后验熵边缘化集合的选择是 PM-MPA 里最“玄学”也最敏感的部分我试过三种准则选择准则做法效果与风险固定度数选列重最小的用户直接进集合实现最简单但没利用信道信息性能损失最不确定信道幅度选 |H(:,j)| 最大的用户信道强则后验可信稳定BER 损失小但每帧信道变化时集合要重新算后验熵预热后计算每用户后验熵选熵最小的用户最符合“已收敛”的语义推荐代价只是多算一次熵我推荐后验熵准则。后验熵小意味着该用户当前的后验已经很尖就算继续迭代也不会翻转。MATLAB 里实现就是% 计算每用户后验熵示例用 f2v 消息做近似 post zeros(M, J); for j 1:J r find(F(:,j)); logPost sum(f2v(:,j,r), 2); logPost logPost - max(logPost); post(:,j) exp(logPost); post(:,j) post(:,j) / sum(post(:,j)); end entropy -sum(post .* log(post eps), 1); [~, sortedIdx] sort(entropy, ascend); margSet sortedIdx(1:Nmarg);注意别把熵最小的节点全选到同一个资源上。比如用户 1 和用户 2 的熵都小但它们都连在资源 1 上选进去只会让资源 1 的组合维度大幅下降其他资源仍然保持 64 的遍历量整体复杂度收益被摊薄。常见做法是加一个约束每个资源最多选 1 个边缘化用户选的时候按熵排序逐个检查。4.3 冻结消息的两种姿态软保持与硬判决冻结消息有两种做法。软保持是把固定后的后验作为 v2f 消息继续参与后续迭代功能节点更新时对冻结用户的所有码字做概率加权复杂度省得不彻底但性能损失更小。硬判决是直接选定后验最大的码字功能节点组合遍历时少一层循环复杂度下降最猛错一个就是错一簇。我的习惯是预热迭代少1 到 2 轮时用软保持兜底迭代较充裕时用硬判决。两者在总迭代 4 到 6 轮时的 BER 差距通常在 0.1 到 0.3 dB复杂度差距却能到 20% 左右。追求实时性或高吞吐仿真就往硬判决靠做性能上限评估就用软保持。4.4 外环调度边缘化集合要不要每轮更新这里有个很多论文没明说的细节边缘化集合是静态的还是动态更新的。静态集合在预热后固定实现简单适合帧长较短、信道慢变的场景动态集合每隔几轮重新评估一次熵把后验变差、出现翻转的节点解冻适合信道质量中等的场景因为后验熵会随迭代加深而变化之前熵小的节点不一定一直可信。动态更新的代价是每轮要额外算一次熵和排序在 MATLAB 里是几十微秒级的开销完全可接受。实际工程中我每 2 轮重新评估一次而不是每轮都做。这样既照顾了消息演变又避免了集合频繁抖动导致的不稳定。注意动态更新时“冻结”的语义会弱化一个节点可能被解冻后重新参与迭代这时它的 v2f 消息要从当前 f2v 重新合成不能沿用冻结前的旧值。4.5 预热轮数边缘化之前给多少传播时间预热轮数决定了消息在因子图上传播的充分程度。预热太少比如只跑 1 轮后验熵还没收敛选出来的边缘化节点往往是“假尖峰”预热太多比如 4 轮以上复杂度优势就提前被消耗掉了。S(4,6) 母矩阵的消息传播半径很短2 轮预热基本能让所有功能节点的影响传递到相邻变量节点。我通常把预热轮数和正式迭代轮数一起调总迭代预算固定时预热占三分之一左右。比如总迭代 6 轮预热 2 轮总迭代 4 轮预热 1 轮。如果用了第 4.4 节的动态更新预热还可以再省一点因为后续每隔 2 轮重新评估集合本身就是一种纠错机制。5. 避坑与常见问题排查从 NaN 发散到错误加速5.1 现象运行几轮迭代后出现 NaNBER 曲线断裂原因对数域消息里出现了 -Inf。初始化时我用了 -log(M) 而不是 -1e8组合枚举时某个组合在所有消息项上都是 -Inflogsumexp 里 exp(-Inf - (-Inf)) 变成 NaN。这条坑在概率域 MPA 里也常见只是转对数域后表现形式从“下溢为 0”变成“NaN 发散”。解决logsumexp 的输入先做一次 max 检查如果所有输入都是 -Inf直接返回 -1e8同时消息数组初始化统一用 -1e8 代替 -Inf。我把这条规则写进所有 MPA 实现函数的头部注释里属于血泪经验。5.2 现象PM-MPA 的 BER 比标准 MPA 差一个量级原因边缘化集合选错了用户。最常见的是选了信道最差的用户其后验熵虽然也小但不是因为“收敛”而是因为该用户信号被强用户压住后验被错误压缩成假尖峰。熵小不等于可信这是 PM-MPA 最容易踩的坑。解决先看预热结束后每个用户的后验熵是“真尖”还是“假尖”。简单判别法是看该用户在最近两轮迭代中硬判决结果是否翻转若翻转次数大于 0就不该边缘化。代码里加一个判定条件marg_flag (entropy thr) (flipCount 0)。5.3 现象静态信道下 BER 正常换成衰落信道后 PM-MPA 完全不能用原因每帧信道都在变边缘化集合如果沿用上一帧的结果就等于用老信道的冻结节点处理新信道的接收信号后面的解码全错。解决每帧重新计算边缘化集合。实现上把“预热 选集合 正式迭代”封装成独立函数帧与帧之间不保留任何边缘化状态。这个生命周期问题很容易被忽略因为静态信道下整个仿真流程跑得太顺很多人根本不会想到集合是“一帧一算”还是“一次算完”。5.4 现象MATLAB 向量化 MPA 循环后速度反而更慢原因向量化版本会把 M^df 的组合矩阵一次性铺开配合 cellfun 或 arrayfun 做扩展。df3、M4 时 64 的组合矩阵很轻松但结合 J6、K4、迭代 6 轮后中间张量会膨胀到上万维内存带宽和函数调用开销已经超过循环本身。解决这个规模下普通 for 循环加预分配数组在 R2021 之后的 MATLAB 里通过 JIT 加速并不慢。真正要优化的是数值运算本身把 exp 计算里的平方项改写成 abs(y - sd)^2 的复用变量避免每次循环都新建复数临时数组。这条排查思路同样适用于其他算法类仿真。5.5 现象硬判决边缘化后LLR 出现波浪形偏差原因硬判决把冻结用户近似为确定码字但该码字对应的星座点与实际发送信号之间存在残差这些残差被当作加性噪声处理导致冻结用户所在资源的 LLR 来回摆动尤其在高 SNR 区域更明显。解决改用软保持或者把硬判决码字对应的残余偏差计入噪声方差N0_eff N0 E|h|^2·|c_fixed - c_true|^2。由于 c_true 未知用该用户后验概率加权的期望误差来近似。这个修正能让硬判决版本在高 SNR 下的 LLR 输出稳定很多。6. 让 PM-MPA 从跑通到可信与标准 MPA 的基线对比习惯PM-MPA 的性能验证我坚持一个死规矩任何一组参数先跑标准 MPA 作为基线再跑 PM-MPA两者用同一个随机种子、同一帧数据。对比内容不只是 BER 曲线还有两个附加量。第一是平均迭代轮数的实测。很多实现把迭代次数写死但 PM-MPA 的收益恰恰体现在提前收敛上。可以在正式迭代里加一个停止条件如果所有未冻结用户的后验熵都低于阈值就提前跳出。实际统计时平均迭代轮数往往比标称值低 15% 到 25%这个数字应该写在实验记录里而不是只报告 BER。第二是复杂度实测。用 timeit 包住检测函数分别记录标准 MPA 和 PM-MPA 在同一组 2000 帧上的耗时再和第 2.3 节的组合数对比确认复杂度下降没有白干。如果耗时下降远小于组合数下降多半是代码里多了不必要的矩阵复制或者 logsumexp 实现写得过重。把“理论组合数”和“实测耗时下降比例”放在同一张表里是排查性能瓶颈最直接的方式。最后别忘了把边缘化集合的选取准则写进实验记录。同一份 PM-MPA 代码用后验熵和用信道幅度选出来的节点完全不同BER 结果也可能差一个量级。我见过不少仿真报告只写“PM-MPA 损失 0.x dB”却不说边缘化怎么选的别人复现时根本对不上数。如果你打算在自己的项目里试用 PM-MPA我的建议是从静态信道 后验熵选取 Nmarg2 起步先确认能复现标准 MPA 的基线再加硬判决和动态集合做优化。这个渐进路线能帮你避开大部分坑也方便回头排查是哪一步引入了性能损失。希望帮到你。本文还有配套的精品资源点击获取