做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 R2023b，Windows和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。它们的基本参数我不想直接念表格，列成对比更直观：
| 参数 | BG1 | BG2 |
|---|---|---|
| 基图维度 | 46行×68列 | 42行×52列 |
| 信息列数K_b | 22 | 10 |
| 最大信息块K | 8448 | 3840 |
| 设计码率范围 | 约1/3～8/9 | 约1/5～2/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时常用K=1408，此时K_b=22，ceil(1408/22)=64，查表正好Zc=64，而且1408+128=1536小于2264=1408？不对，2264=1408，算出来不满足。这时候要换成Zc=128，填充到2816比特，码率会偏低。为了图省事，最稳妥的做法是让K刚好等于K_b*Zc，比如Zc=64时取K=1408，这样完全没有填充位，仿真逻辑最干净。
提示：如果你想让仿真结果和标准传输块贴近，就必须实现填充和速率匹配；但如果目的是验证译码器性能和算法正确性，取K=K_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)*Z+1 : i*Z, (j-1)*Z+1 : j*Z) = block; end end end H = sparse(double(H)); end实际工程里不会用这种双循环去展开大矩阵，因为太慢，但学习阶段这样最直观。展开后H的大小是m_bZc行×n_bZc列，BG1配Zc=64就是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(end+1) = j; end end if length(pivots) ~= M error('H秩不足，无法构造系统生成矩阵'); end % 这里默认主元列就是前M列，实际情况可能需要列重排 A = Hr(:, 1:M); P = Hr(:, M+1:end); K = N - M; G = [gf(eye(K), 1), gf(P', 1)]; G = double(G.x); end这段代码有个需要注意的地方：rref得到的主元列不一定是前M列，我为了可读性先假设前M列是校验位。真正使用时一定要检查pivots是否为1:M，如果不是，需要按主元列做列置换，否则生成矩阵和校验矩阵会对不上。我在自己项目里会加一段自动列重排，但代码一长反而干扰理解，这里只演示核心思路。
对于Zc=64、BG1的4352比特码长，H的行数是2944，生成矩阵大小是1408×4352，double存储大约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_b=22、可选Zc=64，那就需要填充384个比特到1408。填充位一般置0，编码后再把这些填充位当作打孔位不发送。为了简化仿真，我默认信息块长度取K=K_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)}(end+1) = e; adjV{ci(e)}(end+1) = 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_b=22 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、Zc=64这种参数，跑30帧在每个SNR点，大约是几分钟到十几分钟的量级。我建议先把nFrames设成5跑通流程，确认无误后再加大帧数，避免辛苦调半天结果发现是脚本逻辑错误。
4. 仿真结果与参数调优
4.1 不同码块长度和码率的对比
我仿真了BG1和BG2两组配置。BG1用Zc=64、K=1408，码率约0.3235；BG2用Zc=32、K=320，码率约0.1923。只算核心码字、不做速率匹配时，两条曲线在低信噪比下都能正常下降，但BG2因为码率更低，曲线整体更靠左，也就是在同样Eb/N0下误码率更低，这和理论预期一致：低码率纠错能力更强。
有个值得留意的现象：在BER低于1e-3之后，曲线会明显变陡，这是LDPC码的瀑布区特性。瀑布区位置主要取决于码长和列重分布，BG1因为码长更长，瀑布区通常比BG2更陡。如果你仿真时发现曲线没有出现陡降、而是慢慢滑下去，那大概率不是信道问题，而是译码器实现有bug，或者归一化因子没调好。
4.2 迭代次数与归一化因子如何权衡
迭代次数直接影响吞吐和仿真时间。我在BG1、Zc=64下测试过，maxIter从8加到20，误码率改善不到0.1dB，但仿真时间翻了接近一倍。5G系统里对译码时延非常敏感，工程上默认8到10次就够了，这也是我脚本里maxIter=8的原因。做学术研究时你可以放宽到20次，但不要指望靠加大迭代次数去弥补码本设计的缺陷。
归一化因子alpha的调法更有意思。alpha取1就是纯最小和，性能会差一些；取0.5则消息被压得太狠，收敛变慢甚至不收敛。我实测在BG1、BPSK下alpha=0.75到0.8之间最优，BG2短码时0.8更好。一个比较实用的经验是：先跑低信噪比点，把翻译的CRC或校验失败率记录下来，alpha从0.5到0.9每隔0.05扫一遍，哪个点误码率最低就用哪个。一次扫描可能要多跑几帧，但总比凭感觉猜强。
5. 常见问题与实战排查技巧
5.1 译码完全不收敛
如果你发现无论SNR多高，xhat始终满足不了H*xhat=0，优先检查H矩阵的秩。展开出来的H必须行满秩，也就是行数等于秩数。可以用rank(full(H))验证，BG1展开后秩应该等于2944。不满秩的情况通常出现在基图表抄错、移位值读错、或者expandBaseGraph的循环移位方向不一致。
另一个常见原因是LLR的符号约定反了。BPSK调制里我用的映射是0到+1、1到-1，LLR=2y/sigma^2，此时LLR为正意味着比特大概率是0，硬判决取Ltotal<0判为1。如果你把调制映射写成0到-1、1到+1，那LLR符号就要整体取反，否则校验永远不过。
5.2 性能和理论差太远
排除实现错误后，性能差的头号嫌疑是alpha没调，其次是LLR没有按码率折算噪声方差。还有一个容易被忽略的点：填充位处理。如果K小于K_b*Zc但你统计误码时把填充位也算进去，误码率会被稀释，看起来“变好”了，实际是假象。
另外，如果你用了min(次小值)逻辑写错，校验节点会把自身消息算进去，这会让消息在相邻两次迭代里来回震荡。判断方法很简单：把译码器里的校验节点更新替换成纯BP公式（tanh域）跑一遍，如果纯BP性能正常而NMS异常，问题一定在最小和实现，而不是码本或信道。
5.3 大Zc下仿真速度太慢
Zc=128、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%的“性能差”都不是理论上不懂，而是索引、符号方向、填补位这些小地方出了偏差。我每次改完代码都会先用一个极短码长、比如Zc=2甚至更小的自定义基图把流程跑通，确认H、G、译码器三者的尺寸和秩全部正确，再切到完整5G参数。这套流程看着慢，实际是排错最快的路径。
