★ FEATURED ARTICLE 做5G物理层仿真只要和数据信道沾边大概率绕不开LDPC码。从LTE时代的Turbo码到NR时代全面转向LDPC这个转变不是拍脑袋决定的背后牵扯到吞吐量、译码时延、硬件实现复杂度一串工程账。本篇我就用MATLAB把一套完整的5G NR LDPC编译码链路搭出来包括基图选择、准循环校验矩阵展开、生成矩阵求取、归一化最小和译码以及最后的蒙特卡洛误码率统计全部以可运行、可复现的代码为准。这篇文章适合三类人一是通信专业想拿LDPC做课程设计或毕设的同学二是刚进基站或芯片团队、需要快速搞懂5G物理层收发链路的工程师三是自己写Polar码或Turbo码想对比LDPC性能的研究人员。我会尽量把每一步“为什么这么做”都讲透而不是直接甩一段能跑但说不清原理的代码。1. 项目背景与仿真链路设计1.1 5G为什么把宝压在LDPC上先说结论5G NR在数据信道用LDPC控制信道用Polar码两者分工明确。LDPC码的最大优势是译码器天然适合并行化硬件上可以做成高吞吐的迭代译码器这对5G动辄上百Mbps甚至Gbps级别的峰值速率来说是硬指标。Turbo码虽然在中短码长下性能不差但它的迭代译码结构对并行不友好高吞吐场景下功耗和面积都不划算。LDPC的另一个特性是灵活。5G NR通过两张基图加上一系列提升因子就能覆盖从几十比特到八千多比特的传输块长度码率也从1/3附近一路拉到8/9这对移动通信里多变的车速、信道质量、资源调度粒度非常重要。理解这些参数背后的规则是仿真实现的第一步也是最容易忽略的一步。1.2 一套可复用的仿真链路框架整个仿真链路我按照发射机、信道、接收机三块来搭。发射机做信息比特填充、LDPC编码、BPSK调制信道加高斯白噪声AWGN接收机做LLR软解调、LDPC迭代译码、硬判决和误码率统计。之所以用BPSK而不是QPSK或更高阶调制是因为LDPC译码性能本身和调制方式基本解耦BPSK最容易把误码率曲线和信道容量极限对比。脚本流程如下选基图BG1或BG2确定信息块长度K和提升因子Zc展开校验矩阵H由H求生成矩阵G完成编码按Eb/N0扫描多个信噪比点每个点跑若干帧统计误比特率BER和误帧率FER绘图。这套框架最大的好处是模块化你可以只替换译码器、只替换调制方式、只替换信道模型其他环节不用动。我手里的版本用普通MATLAB脚本就能跑不需要额外安装通信工具箱和5G工具箱但后文也会介绍怎么用工具箱函数交叉验证结果。1.3 运行环境和代码约定我用的环境是MATLAB R2023bWindows和Linux两种系统都跑过。代码里没有依赖特殊工具箱的自写部分核心代码用原生矩阵运算完成R2019b之后应该都能跑。如果只是做性能验证R2021a之后的版本就够了如果你还想和5G Toolbox自带的nrLDPCEncode对比那需要安装5G Toolbox。一个重要的代码约定所有二进制域上的运算都统一用mod(x, 2)处理不要依赖逻辑运算的隐式转换否则一旦矩阵里有非0/1值错误会非常隐蔽。译码器内部的消息全部用double类型不刻意转成定点这样可以先把算法逻辑验证清楚后续做硬件定点化时再单独处理。2. 5G NR LDPC的关键参数与展开过程2.1 BG1/BG2和提升因子Zc怎么选5G NR标准里定义了两张基图基图1和基图2分别叫BG1和BG2。它们的基本参数我不想直接念表格列成对比更直观参数BG1BG2基图维度46行×68列42行×52列信息列数K_b2210最大信息块K84483840设计码率范围约1/38/9约1/52/3典型使用场景大传输块、高码率小传输块、低码率选BG1还是BG2标准里主要看信息块大小和码率。简单粗暴的记忆方式信息块比较大、码率要求高用BG1信息块小到一两百比特、或者码率压得很低用BG2。实际仿真里如果你想仿真满载的移动宽带业务选BG1想做URLLC这类短包业务选BG2更贴近真实场景。提升因子Zc是准循环LDPC的核心参数。标准规定的取值不是连续整数而是一组离散值我常用的是下面这串Zlist [2,4,8,16,32,64,128,256, ... 3,6,12,24,48,96,192,384, ... 5,10,20,40,80,160,320, ... 7,14,28,56,112,224, ... 9,18,36,72,144,288, ... 11,22,44,88,176,352, ... 13,26,52,104,208, ... 15,30,60,120,240];Zc的选取规则是先计算ceil(K/K_b)在Zlist里找第一个大于等于它的值然后要满足K 2Zc K_bZc。第二个条件很多人会漏它的物理意义是保证系统列足够编码不会出现填充太多导致码率严重损失。我仿真BG1时常用K1408此时K_b22ceil(1408/22)64查表正好Zc64而且14081281536小于22641408不对22641408算出来不满足。这时候要换成Zc128填充到2816比特码率会偏低。为了图省事最稳妥的做法是让K刚好等于K_b*Zc比如Zc64时取K1408这样完全没有填充位仿真逻辑最干净。提示如果你想让仿真结果和标准传输块贴近就必须实现填充和速率匹配但如果目的是验证译码器性能和算法正确性取KK_b*Zc是最省心、也最容易排错的做法。2.2 由基矩阵展开校验矩阵5G NR的基图矩阵在3GPP TS 38.212里有完整定义里面的元素是三类-1、0、和正整数值。-1代表Zc×Zc的全零块0代表单位矩阵且不移位正整数代表单位矩阵按该值循环右移。展开的过程就是遍历基图的m_b行和n_b列每个元素生成一个Zc×Zc的子块。展开代码很简单我直接贴function H expandBaseGraph(BG, Z) [mb, nb] size(BG); H logical(zeros(mb*Z, nb*Z)); for i 1:mb for j 1:nb s BG(i, j); if s 0 block eye(Z); if s 0 block circshift(block, [0, s]); end H((i-1)*Z1 : i*Z, (j-1)*Z1 : j*Z) block; end end end H sparse(double(H)); end实际工程里不会用这种双循环去展开大矩阵因为太慢但学习阶段这样最直观。展开后H的大小是m_bZc行×n_bZc列BG1配Zc64就是2944行×4352列稀疏矩阵存储下内存完全不是问题。这一步做完之后H的每一行就是一个校验方程每一列对应一个码字比特后续译码器就是在这个稀疏结构上跑消息迭代。需要特别提醒的是循环移位方向在标准文档里写得很死不同教材可能采用左移或右移。做仿真时不一定要和标准逐位对齐但如果你后面要用5G Toolbox的编码器做交叉验证就必须把移位方向、尤其是速率匹配里的比特交织规则统一。否则自写译码器明明是对的却和工具箱对不上排查起来很痛苦。2.3 译码器的理论根基BP与最小和置信传播译码的本质是Tanner图上的消息传递。简单说变量节点代表码字比特校验节点代表校验方程边代表它们之间的关系。译码开始时每个变量节点持有来自信道的对数似然比LLR然后通过边来回收发校验节点的可靠性信息迭代更新。在BP算法里校验节点的更新要算tanh函数硬件和定点实现都很不方便。所以实际系统基本都用最小和近似也就是把校验节点传给变量节点的消息简化成取所有入边绝对值的最小值再乘上符号。由于这种近似会让消息偏大就产生了两个常见修正归一化最小和NMS和偏移最小和OMS。我用的是NMS。公式上校验节点j发给变量节点i的消息等于U(j,i) alpha * (所有相邻变量消息符号之积) * (排除i后消息绝对值的最小值)alpha就是归一化因子一般取0.75附近。alpha太小会过度惩罚消息太大又退化成纯最小和、失去修正效果。这点后续调参会细讲。3. MATLAB实现编码器、译码器和主仿真3.1 从校验矩阵求生成矩阵LDPC编码最直接的做法是求出生成矩阵G然后让信息比特u左乘G得到码字c。G和H满足H*G 0二进制域上。求G的标准方法是把H通过高斯消元化为系统形式[A I]然后G [I A]。我只在二进制域上实现高斯消元用MATLAB自带的gf对象最省事function G computeGeneratorMatrix(H) [M, N] size(H); Hx gf(full(H), 1); Hr rref(Hx); % 模2行最简形 Hr double(Hr.x); % 找到主元列 pivots []; for j 1:N if any(Hr(:, j) 1) pivots(end1) j; end end if length(pivots) ~ M error(H秩不足无法构造系统生成矩阵); end % 这里默认主元列就是前M列实际情况可能需要列重排 A Hr(:, 1:M); P Hr(:, M1:end); K N - M; G [gf(eye(K), 1), gf(P, 1)]; G double(G.x); end这段代码有个需要注意的地方rref得到的主元列不一定是前M列我为了可读性先假设前M列是校验位。真正使用时一定要检查pivots是否为1:M如果不是需要按主元列做列置换否则生成矩阵和校验矩阵会对不上。我在自己项目里会加一段自动列重排但代码一长反而干扰理解这里只演示核心思路。对于Zc64、BG1的4352比特码长H的行数是2944生成矩阵大小是1408×4352double存储大约48MB在普通电脑上可以接受。如果码长再翻倍建议改用稀疏逻辑矩阵或直接实现结构化编码不要硬怼生成矩阵。3.2 编码器以及填充位处理有了G之后编码就是一次矩阵乘function cw ldpcEncode(info, G) % info: K×1 二进制列向量 cw mod(double(info) * double(G), 2); cw cw(:); end我再强调一个传输块层面的细节。标准里实际信息块长度K往往不是K_bZc比如你想传1024比特时K_b22、可选Zc64那就需要填充384个比特到1408。填充位一般置0编码后再把这些填充位当作打孔位不发送。为了简化仿真我默认信息块长度取KK_bZc也就是没有填充。如果你确实要模拟K小于K_bZc的情况记得两点第一编码前在信息比特后面补K_bZc-K个0第二译码完成后统计误码时只比较原始K个信息比特不要比较填充位否则误码率会被稀释。3.3 归一化最小和译码器实现译码器是整个仿真里最核心、也最容易写错的部分。我先贴完整代码再逐一解释。function [xhat, ok] ldpcDecodeNMS(llrIn, H, maxIter, alpha) [mb, nb] size(H); [ri, ci, ~] find(H); E length(ri); if E 0 error(校验矩阵为空); end LLR llrIn(:); adjC cell(mb, 1); adjV cell(nb, 1); for e 1:E adjC{ri(e)}(end1) e; adjV{ci(e)}(end1) e; end U zeros(E, 1); M zeros(E, 1); for it 1:maxIter % 变量节点更新 for v 1:nb eV adjV{v}; if isempty(eV), continue; end total LLR(v) sum(U(eV)); M(eV) total - U(eV); end % 校验节点更新 for c 1:mb eC adjC{c}; if isempty(eC), continue; end sgnProd prod(sign(M(eC))); minAbs min(abs(M(eC))); for e eC s sgnProd * sign(M(e)); if abs(M(e)) minAbs tmp abs(M(eC)); tmp(find(tmp minAbs, 1)) inf; U(e) alpha * s * min(tmp); else U(e) alpha * s * minAbs; end end end % 硬判决 Ltotal LLR; for v 1:nb if ~isempty(adjV{v}) Ltotal(v) Ltotal(v) sum(U(adjV{v})); end end xhat double(Ltotal 0); % 校验是否满足全部方程 if mod(H * xhat, 2) 0 ok true; return; end end ok false; end变量节点更新这步的含义是某个变量节点要发给某个校验节点的消息等于它自己的信道LLR加上除该校验节点之外其他所有校验节点给它的消息。代码里先算total LLR(v) sum(U(eV))再对每条边做total - U(e)就是标准的“排除自身”。校验节点更新的难点在“排除自身后的最小值”。我这里的做法是先把所有入边绝对值存进tmp如果当前边的绝对值恰好等于最小值就把其中一个最小值替换成Inf再取min得到次小值。如果最小值有重复这个逻辑依然正确因为只替换了一个。这是最小和实现里最容易出错的位置很多人的译码器性能差就是在这里用错了索引。最后每轮迭代都做一个校验H乘以硬判决结果是否为全零。如果全零说明译码成功立即跳出。对于BPSK和高斯信道如果SNR足够高往往两三轮迭代就收敛了。3.4 蒙特卡洛误码率脚本主仿真脚本用一个双层循环外层扫Eb/N0内层跑多帧。下面这段是我实际在用的简化版clear; clc; % 参数 bg 1; Zc 64; K 22 * Zc; % BG1时 K_b22 mb 46; nb 68; N nb * Zc; R K / N; EbN0dB 0:0.5:4.5; nFrames 30; maxIter 8; alpha 0.75; % 展开H并求G BG loadBaseGraph(bg); % 从3GPP表导入基图 H expandBaseGraph(BG, Zc); G computeGeneratorMatrix(H); ber zeros(size(EbN0dB)); fer zeros(size(EbN0dB)); for idx 1:length(EbN0dB) EbN0 10^(EbN0dB(idx)/10); sigma sqrt(1 / (2 * R * EbN0)); % BPSK errBits 0; errFrames 0; for fr 1:nFrames info randi([0 1], K, 1); cw ldpcEncode(info, G); x 1 - 2*cw; % 0 - 1, 1 - -1 y x sigma * randn(N, 1); llr 2 * y / sigma^2; [xhat, ok] ldpcDecodeNMS(llr, H, maxIter, alpha); errBits errBits sum(xhat(1:K) ~ info); if sum(xhat(1:K) ~ info) 0 errFrames errFrames 1; end end ber(idx) errBits / (nFrames * K); fer(idx) errFrames / nFrames; end semilogy(EbN0dB, ber, -o, EbN0dB, fer, -s); grid on; xlabel(Eb/N0 (dB)); ylabel(BER/FER); legend(BER, FER);噪声方差公式里R的出现在于BPSK每符号能量和比特能量的换算。我们发的是±1符号每个符号携带1个比特但每比特实际能量等于总符号能量除以码率所以噪声方差要按信噪比和码率一起折算。如果你把R写成1算出来的曲线会整体偏移这是新手上路常见的错。对于BG1、Zc64这种参数跑30帧在每个SNR点大约是几分钟到十几分钟的量级。我建议先把nFrames设成5跑通流程确认无误后再加大帧数避免辛苦调半天结果发现是脚本逻辑错误。4. 仿真结果与参数调优4.1 不同码块长度和码率的对比我仿真了BG1和BG2两组配置。BG1用Zc64、K1408码率约0.3235BG2用Zc32、K320码率约0.1923。只算核心码字、不做速率匹配时两条曲线在低信噪比下都能正常下降但BG2因为码率更低曲线整体更靠左也就是在同样Eb/N0下误码率更低这和理论预期一致低码率纠错能力更强。有个值得留意的现象在BER低于1e-3之后曲线会明显变陡这是LDPC码的瀑布区特性。瀑布区位置主要取决于码长和列重分布BG1因为码长更长瀑布区通常比BG2更陡。如果你仿真时发现曲线没有出现陡降、而是慢慢滑下去那大概率不是信道问题而是译码器实现有bug或者归一化因子没调好。4.2 迭代次数与归一化因子如何权衡迭代次数直接影响吞吐和仿真时间。我在BG1、Zc64下测试过maxIter从8加到20误码率改善不到0.1dB但仿真时间翻了接近一倍。5G系统里对译码时延非常敏感工程上默认8到10次就够了这也是我脚本里maxIter8的原因。做学术研究时你可以放宽到20次但不要指望靠加大迭代次数去弥补码本设计的缺陷。归一化因子alpha的调法更有意思。alpha取1就是纯最小和性能会差一些取0.5则消息被压得太狠收敛变慢甚至不收敛。我实测在BG1、BPSK下alpha0.75到0.8之间最优BG2短码时0.8更好。一个比较实用的经验是先跑低信噪比点把翻译的CRC或校验失败率记录下来alpha从0.5到0.9每隔0.05扫一遍哪个点误码率最低就用哪个。一次扫描可能要多跑几帧但总比凭感觉猜强。5. 常见问题与实战排查技巧5.1 译码完全不收敛如果你发现无论SNR多高xhat始终满足不了H*xhat0优先检查H矩阵的秩。展开出来的H必须行满秩也就是行数等于秩数。可以用rank(full(H))验证BG1展开后秩应该等于2944。不满秩的情况通常出现在基图表抄错、移位值读错、或者expandBaseGraph的循环移位方向不一致。另一个常见原因是LLR的符号约定反了。BPSK调制里我用的映射是0到1、1到-1LLR2y/sigma^2此时LLR为正意味着比特大概率是0硬判决取Ltotal0判为1。如果你把调制映射写成0到-1、1到1那LLR符号就要整体取反否则校验永远不过。5.2 性能和理论差太远排除实现错误后性能差的头号嫌疑是alpha没调其次是LLR没有按码率折算噪声方差。还有一个容易被忽略的点填充位处理。如果K小于K_b*Zc但你统计误码时把填充位也算进去误码率会被稀释看起来“变好”了实际是假象。另外如果你用了min(次小值)逻辑写错校验节点会把自身消息算进去这会让消息在相邻两次迭代里来回震荡。判断方法很简单把译码器里的校验节点更新替换成纯BP公式tanh域跑一遍如果纯BP性能正常而NMS异常问题一定在最小和实现而不是码本或信道。5.3 大Zc下仿真速度太慢Zc128、256时H矩阵变大边数上千条双循环译码会明显变慢。加速手段有三个一是把变量节点更新改成按变量节点一次向量化不要一条边一条边地写二是把校验节点里找次小值的部分改成先对整个校验节点排序再统一生成消息三是用预分配的邻接表代替动态增长的cell扩展。还有一个工程技巧在同一轮SNR扫描中如果上一帧译码在2次迭代就收敛了下一帧初始消息可以用上一帧的最终消息热启动但这只适用于连续信道相干时间内的仿真普通蒙特卡洛里不用这么麻烦直接把nFrames调小更省事。5.4 与5G Toolbox交叉验证时对不上如果你装了5G Toolbox可以用nrLDPCEncode和nrLDPCDecode来验证自写代码。但我遇到过非常坑的情况自写展开的H没问题、编码结果却和nrLDPCEncode不一致原因是标准里的速率匹配和比特选择过程做了额外的交织工具箱会自动应用这些规则而自写代码只做了最原始的核心编码。这不代表自写代码错了只是你们处理的码字定义不在同一个层面。交叉验证的稳妥做法是先把K设成K_b*Zc不使用任何速率匹配只比较核心码字部分。此时如果生成矩阵正确自写编码结果应该和nrLDPCEncode输出在去掉填充后逐位一致。如果还有差异检查基图表里移位值的行列下标是从0还是从1开始这类索引问题在36系协议里永远防不胜防。最后一个实际体会LDPC仿真调试里90%的“性能差”都不是理论上不懂而是索引、符号方向、填补位这些小地方出了偏差。我每次改完代码都会先用一个极短码长、比如Zc2甚至更小的自定义基图把流程跑通确认H、G、译码器三者的尺寸和秩全部正确再切到完整5G参数。这套流程看着慢实际是排错最快的路径。
阅读完成 · 觉得有帮助？
