LDPC码的编码实现是很多刚开始做信道编码仿真的同学绕不过去的一个坎。LDPC码本身性能逼近香农极限在5G新空口、Wi-Fi 6、DVB-S2这些标准里都是主力信道编码方案。但真正动手写Matlab仿真的时候你会发现译码端调通还算容易编码端反而最卡人。直接用生成矩阵G去编码H矩阵稀疏的优势全部丢掉一个码长上千的矩阵乘法在普通笔记本上能跑到你怀疑人生。这篇文章就把我用Matlab实现LDPC码近似下三角编码的完整过程讲清楚包括原理推导、代码实现、误码率曲线绘制还有我在仿真里踩过的那些坑。适合正在做LDPC仿真、或者准备把LDPC写进自己通信系统里的朋友参考。1. LDPC编码到底难在哪里1.1 LDPC为什么在5G时代“翻红”LDPC码是上世纪60年代提出的之后沉寂了很久直到90年代末才被重新发掘。它的核心优势在于校验矩阵H是稀疏的配合置信传播BP迭代译码性能可以非常接近香农限。相比Turbo码LDPC的译码并行度更高不需要交织器错误平层也更低所以在5G增强移动宽带场景里被选为数据信道的编码方案。不过大家平时接触LDPC基本都是从译码开始的。译码器拿到对数似然比LLR在Tanner图上做消息传递迭代几次以后硬判决得到码字。这个过程只需要H矩阵参与运算H矩阵的稀疏性让每次迭代的消息更新都只涉及少量边复杂度可控。但编码端不一样。编码的本质是要找到满足Hc^T 0的码字c这意味着你得先有办法从信息位s构造出校验位p。这就是LDPC编码比译码难的地方。1.2 从H到G是一条“死路”最直观的编码方式是把校验矩阵H通过高斯消元化成系统形式 [I P] 或者 [P I]然后得到生成矩阵G编码时直接算 c sG。这在理论上是完全没有问题的但实际做仿真的时候这一套基本跑不起来。原因很简单你的LDPC码为了保证性能H矩阵设计成稀疏的但H经高斯消元化成系统形式以后得到的 [I P] 中P部分往往是稠密的。稠密矩阵乘法的复杂度是O(n²)n是码长。算一笔账码长n 2048码率1/2一次编码需要做大约1000×1000量级的乘加运算。仿真时要跑成千上万帧还要在多个信噪比点下重复这个计算量是灾难性的。要是码长再上到4096或者8192普通的Matlab脚本根本没法在合理时间内完成。很多人第一次做LDPC仿真就卡在这里译码器已经写好了BER曲线却跑不出来因为编码端太慢。我当时也是这样后来才意识到必须绕开生成矩阵G直接从稀疏的H矩阵出发做编码。1.3 从H矩阵直接编码的思路既然不能显式构造G那就想办法直接在H上做文章。校验方程 Hc^T 0 等价于让码字分成信息位和校验位两组然后通过H的分块结构解出校验位。关键在于我们希望这个过程借助H的稀疏结构以接近线性的复杂度完成。如果H能写成一个近似下三角的形式例如H [ A B T ; C D E ]其中T是一个下三角矩阵那么编码过程就可以先从信息位推算出第一段校验位再用回代方式解出第二段校验位。这个过程的复杂度主要取决于T的大小和矩阵运算方式整体可以控制在O(n)到O(n²)之间远好于直接乘稠密G。这就是近似下三角编码方法的基本出发点。2. 近似下三角编码的原理与推导2.1 什么是近似下三角形式近似下三角编码也常被称作RU编码方法核心思路不是把H完全化成系统形式而是只做有限的行列置换让H变成一个带“缺口”的下三角结构。具体来说经过行列置换H被分成六块H [ A B T ; C D E ]这里矩阵的总列数是n总行数是m n - k。T是一个(m - g) × (m - g)的下三角矩阵g称为gap缺口大小表示右下角没有完全下三角化的列数。A的大小是(m - g) × kB是(m - g) × gC是g × kD是g × gE是g × (m - g)。如果g 0T本身是满秩的下三角矩阵编码会非常干净。但实际中随机生成的稀疏LDPC校验矩阵经过简单的行列置换不一定能让g 0通常需要允许一个小的gap比如g 16或者32。这个gap越大编码复杂度越高所以预处理的目标就是尽量让gap变小。2.2 关键公式推导从分块方程解出校验位设码字向量c分为三段c [s, p1, p2]其中s是信息位长度为kp1是第一部分校验位长度为gp2是第二部分校验位长度为m - g。把H和c代到校验方程Hc^T 0里按分块展开得到两个方程A s^T B p1^T T p2^T 0 C s^T D p1^T E p2^T 0先看第一个方程。因为T是下三角矩阵所以可以从下往上做回代用s和p1表达p2p2^T −T^{-1} (A s^T B p1^T)注意这里并不会真的计算T的逆矩阵而是利用T的下三角结构做前向代入复杂度只有O((m - g)²)量级。把这个p2代回第二个方程整理后得到(C − E T^{-1} A) s^T (D − E T^{-1} B) p1^T 0令 φ D − E T^{-1} B则有p1^T −φ^{-1} (C − E T^{-1} A) s^T也就是说只要φ矩阵可逆p1就能直接由信息位s算出来p1出来以后p2用一次回代就能解出。整个编码过程就完成了。2.3 gap值的选择与计算量影响从上面的公式可以看出gap值g直接决定复杂度。φ是一个g×g矩阵求φ的逆的复杂度是O(g³)。回代求解p2的复杂度是O((m - g)²)。整体编码复杂度大约是O(n g³)。如果g很小比如g 16那涉及φ的那些矩阵运算几乎可以忽略但如果g变得很大比如g 128φ的逆运算就会成为瓶颈。所以预处理阶段的核心目标就是在保持H稀疏性的前提下把g压缩到尽量小。我在实际仿真中一般接受g ≤ 32这个量级。对于码长256到2048的仿真这个规模下的编码时间完全可接受。2.4 为什么不用完全下三角你可能想问既然近似下三角这么绕直接想办法把H化成完全下三角不就行了吗理论上可以实际很难。完全下三角要求g 0这意味着校验矩阵本身就具备一种非常规律的结构。常规的规则LDPC码随机生成出来经过行列置换很难同时满足下三角和稀疏性两个约束。为了强行让g 0往往要做大量的行变换这个过程中H矩阵的稀疏性会被严重破坏到头来编码快了但矩阵变稠密存储和计算又变成新的瓶颈。近似下三角编码的意义就在这里稍微容忍一点点gap换来的代价其实很小但保留了H的稀疏性整体编码过程依然高效。用生活里的话说这是“花小钱办大事”的思路。3. Matlab代码实现把公式变成能跑的仿真3.1 生成稀疏校验矩阵H仿真第一步得先有一副“牌”可用。我常用的是规则码列重dv 3行重dc 6码率1/2。生成方式很多这里给一种最容易理解的下标随机交织构造法它生成的是一个行重列重都固定的正则H矩阵。function H generateH(n, dc, dv) % n: 码长 % dc: 行重, dv: 列重 % 返回 m x n 的 GF(2) 稀疏校验矩阵 m n * dv / dc; % 校验方程数 if mod(m, 1) ~ 0 error(参数必须满足 n*dv 能被 dc 整除); end H zeros(m, n); % 逐列放置1的位置保证每列有dv个1每行有dc个1 for j 1:n colIdx zeros(1, dv); for t 1:dv % 随机选择当前行重小于dc的行 candidates find(sum(H, 2) dc); colIdx(t) candidates(randi(length(candidates))); end H(colIdx, j) 1; end % 以防万一清理偶环简单策略接受现有矩阵后续译码验证 H sparse(H); end这段代码的思路是逐列确定1的位置每列随机放到行重还没满的行上。生成结果可能带有少量长度为4的环对于教学仿真影响不大。你如果追求更优的性能可以在这个基础上加环检测和重排逻辑后面再聊。3.2 预处理把H变成近似下三角结构接下来是核心的预处理环节。我直接给一个简化版的高斯消元找下三角块的方法。它的思想是选一部分列作为“主元列”对这些列做GF(2)高斯消元尽量把它们约化成下三角结构记录行置换和列置换gap就是消元过程中没能消成主元的那几列。function [Htilde, rowPerm, colPerm, gap] ru_preprocess(H) % RU预处理将H变换为近似下三角形式 % 返回 Htilde H(rowPerm, colPerm) % 使得 Htilde [A B T; C D E]T为下三角 H mod(H, 2); [m, n] size(H); rowPerm 1:m; colPerm 1:n; Htmp H; % 选主元列从右往左构造下三角 % 这里使用贪心策略优先选择当前可用列 pivCols []; % 主元列位置在当前Htmp中的列号 for pivotRow m:-1:1 % 在当前剩余矩阵中找一列在pivotRow行有1且在这行之前尽量简单 colCand find(Htmp(pivotRow, :) 1); if isempty(colCand) continue; end % 简单贪心取第一个候选列 selCol colCand(1); pivCols [pivCols, selCol]; % 消去其他行在该列上的1 rowsHaveOne find(Htmp(:, selCol) 1); rowsHaveOne(rowsHaveOne pivotRow) []; for r rowsHaveOne(:) Htmp(r, :) mod(Htmp(r, :) Htmp(pivotRow, :), 2); end end % 将主元列移到矩阵后端非主元列放到前端作为信息位和gap部分 nonPiv setdiff(1:n, pivCols); newColOrder [nonPiv, fliplr(pivCols)]; % pivCols从右下到右上fliplr保持顺序 colPerm newColOrder; % 行也重新排序将pivotRow序列恢复到自然顺序 [~, rowOrder] sort(pivCols); % 这个细节需要结合循环内记录优化 rowPerm rowOrder; Htilde H(rowPerm, colPerm); gap n - length(pivCols); % 非主元列数量这里作为近似gap end这段代码在我实际使用中做了一个关键简化它没有记录每步消元对应的行置换细节所以如果你直接复制运行行置换部分可能需要再处理一下。我更推荐的做法是把预处理视为“离线操作”在Matlab里跑一次把rowPerm、colPerm和gap保存成.mat文件后续每次仿真直接加载不必重复计算。我在仿真中就是这么组织的省了很多时间。真的要写一个严格正确的RU预处理核心是维护一个“未消元行”的集合每次选择一个合适的列把它和某一行交换位置然后对矩阵的其余部分做消元。因为涉及稀疏矩阵的频繁动态操作在Matlab里用全矩阵会更简单但码长超过2000以后会变慢。我觉得对于教学仿真直接预处理一次然后保存结果就够了没必要把预处理本身优化成生产级代码。3.3 编码函数信息位到完整码字有了近似下三角矩阵后编码函数可以直接按公式实现。这里要注意一个特别容易出错的地方预处理时对列做了置换等价于改变了码字比特的排列顺序。所以编码时信息位要先按照colPerm对应的顺序摆放编码完成后再把码字恢复到原始顺序。function [c, p1, p2] ldpc_encode_ru(s, H, rowPerm, colPerm, gap) % s: 原始信息位向量 (1 x k) % 返回完整的码字c原始比特顺序 [m, n] size(H); k n - m; % 根据colPerm把信息位放到置换后的对应位置 % 在列置换后前面部分是信息位 gap校验位后面部分是T对应的校验位 % 这里假设colPerm前k列为信息位后m列中gap个为p1其余为p2 Htilde H(rowPerm, colPerm); % 分割矩阵 A Htilde(1:end-gap, 1:k); B Htilde(1:end-gap, k1:kgap); T Htilde(1:end-gap, kgap1:end); C Htilde(end-gap1:end, 1:k); D Htilde(end-gap1:end, k1:kgap); E Htilde(end-gap1:end, kgap1:end); % 计算phi D - E * inv(T) * B实际用回代求解 X solve_tri(T, B); % X T^-1 B phi mod(D - E * X, 2); if rcond(double(phi)) 1e-10 error(phi矩阵接近奇异请重新生成H或者增大gap); end Y solve_tri(T, A); % Y T^-1 A % p1 -inv(phi) * (C * s - E * Y * s) s_perm s; % 简化处理信息位按原顺序实际应遵循colPerm p1 mod(-phi \ (mod(C * s_perm., 2) - E * mod(Y * s_perm., 2)), 2).; p1 mod(p1, 2); % p2 -inv(T) * (A * s B * p1) rhs mod(A * s_perm. B * p1., 2); p2 mod(-solve_tri(T, rhs), 2).; % 组装置换后的码字 c_perm [s_perm, p1, p2]; % 恢复原始顺序 c zeros(1, n); c(colPerm) c_perm; c mod(c, 2); end function X solve_tri(T, B) % 利用T的下三角结构解 X T^-1 B % 这里直接写满矩阵回代教学用 [r, c] size(T); X zeros(r, size(B, 2)); % 简单方式直接使用inv但代码可读性好 % 实际工程中应用向前代入 X T \ B; X mod(X, 2); end在上面的代码中关于信息位顺序的问题我做了转换在预处理时前三列已经假定是信息位集合。实际应用中你需要在预处理导出colPerm时就确定哪些列是信息位哪些列是校验位。一个比较方便的做法是选择一个初始列顺序让前k列是信息位后m列是校验位做完RU预处理后把前k列的原始列号存下来。这样编码时代码里s_perm s(:, colPerm(1:k))逻辑更清晰。为了避免混淆我这里省略了严格置换恢复但这部分值得你花时间仔细设计一次因为它是整个编码实现的“暗坑”之一。3.4 调制、AWGN信道与软信息计算编码完成后下一步是BPSK调制。我规定码字比特0映射到11映射到-1。这样接收端y x noise。噪声方差按照Eb/N0换算公式是sigma² 1 / (2 * R * 10^(EbN0/10))其中R k/n是码率。收到y后进入译码器的信道LLR初始值为Lch 2 * y / sigma²。这个公式的来源是BPSK下高斯信道LLR推导记得符号要对齐。function llr get_channel_llr(y, sigma2) llr 2 * y / sigma2; end3.5 LLR-BP译码器实现译码器我用置信传播的全信息版本。变量节点更新、校验节点更新轮流迭代。校验节点更新用tanh规则公式简洁但计算量偏大实际仿真时我喜欢先用Min-Sum快速验证链路然后再切回标准BP出性能曲线。下面给一个标准的对数域BP实现。function [decoded, iter] ldpc_bp_decode(llr, H, maxIter) % llr: 信道初始对数似然比 (1 x n) % H: 校验矩阵稀疏 [m, n] size(H); [rows, cols] find(H); % 校验矩阵非零位置 numEdges length(rows); % 初始化校验节点到变量节点的消息为0 msg_c2v zeros(1, numEdges); % 变量节点到校验节点消息初始为信道信息 msg_v2c llr; % 为方便消息传递建立索引映射表 % 这里用cell数组存每条边对应的变量节点和校验节点 v2cMsg zeros(numEdges, 1); for iter 1:maxIter % 校验节点更新 for e 1:numEdges % 实际Matlab应向量化教学演示用循环 varIdx cols(e); % 该边连接的变量节点 checkIdx rows(e); % 该边连接的校验节点 % 收集除当前边外所有邻居消息 neighbors find(rows checkIdx cols ~ varIdx); if isempty(neighbors) msg_c2v(e) llr(varIdx); continue; end prod_tanh 1; for nb neighbors(:) prod_tanh prod_tanh * tanh(v2cMsg(nb) / 2); end msg_c2v(e) 2 * atanh(prod_tanh); end % 变量节点更新 for e 1:numEdges varIdx cols(e); checkIdx rows(e); neighbors find(cols varIdx rows ~ checkIdx); msg_v2c_sum llr(varIdx); for nb neighbors(:) msg_v2c_sum msg_v2c_sum msg_c2v(nb); end v2cMsg(e) msg_v2c_sum; end % 计算总对数似然比并判决 totalLLR llr; for v 1:n edgesOfVar find(cols v); for e edgesOfVar(:) totalLLR(v) totalLLR(v) msg_c2v(e); end end decoded double(totalLLR 0); % LLR 0 判为比特1 if mod(H * decoded., 2) 0 return; end end end这里面的循环写法只是方便理解真正跑仿真的时候一定要做向量化否则码长512时迭代50次每一帧都会慢得让人抓狂。我在下面的“常见问题”里会专门讲如何把这段改成向量化实现。实际上你在实际调试时可以先跑通这个循环版再优化速度。3.6 BER主程序串联整个仿真链路最后是主程序。外层遍历信噪比内层逐帧统计错误比特。我习惯每个信噪比点累积到一定数量的错误比特再停止避免高信噪比时曲线出现明显的“抖动断层”。% 参数设置 n 512; k 256; m n - k; dc 6; dv 3; maxIter 50; EbN0_dB 0:0.5:3; numFramesMax 200; H generateH(n, dc, dv); % 假设置换已预处理完成这里省略ru_preprocess调用 % 实际仿真时从文件加载 RowPerm, ColPerm, gap for iSnr 1:length(EbN0_dB) EbN0 10^(EbN0_dB(iSnr)/10); sigma2 1 / (2 * k/n * EbN0); errorBits 0; totalBits 0; frame 0; while frame numFramesMax s randi([0,1], 1, k); [c, ~, ~] ldpc_encode_ru(s, H, rowPerm, colPerm, gap); % BPSK调制 x 1 - 2 * c; % AWGN信道 y x sqrt(sigma2) * randn(1, n); % 译码 llr get_channel_llr(y, sigma2); [decoded, ~] ldpc_bp_decode(llr, H, maxIter); errorBits errorBits sum(decoded ~ c); totalBits totalBits n; frame frame 1; if errorBits 200 break; end end ber(iSnr) errorBits / totalBits; end semilogy(EbN0_dB, ber, o-) xlabel(Eb/N0 (dB)) ylabel(BER) grid on注意这段代码把 ru_preprocess 的调用省略了你实际运行的时候需要先把 rowPerm、colPerm 和 gap 生成好。比较稳妥的做法是先跑一次预处理函数把它们保存下来[Htilde, rowPerm, colPerm, gap] ru_preprocess(H); save(ldpc_ru_params.mat, rowPerm, colPerm, gap);后续所有信噪比点的仿真都用同一份置换不要每次重新生成否则编码和译码的比特顺序对不上。4. 仿真参数设计与误码率曲线解读4.1 一套可复现的仿真参数我给出一套我实测过的参数表适合在没有专门硬件的情况下验证链路和出曲线。参数数值说明码长 n512太小则性能离理论差得远太大则仿真慢信息位长度 k256码率0.5列重 dv3变量节点度行重 dc6校验节点度最大迭代次数50典型值继续增大收益递减信噪比范围0-3 dB瀑布区主要分布范围每点最大帧数200或按错误比特数控制这套参数在普通笔记本上用向量化BP译码每个信噪比点大概几十秒到几分钟整体曲线一个小时内能跑完。如果你第一次调试建议先把码长降到256信噪比点数减少到4个先确认链路通不通。4.2 Eb/N0与噪声方差换算前面提到过换算公式这里再展开说一句。很多新手直接在信道上用固定的复数噪声功率结果发现BER曲线和理论值差得很远。原因在于LDPC码率不是1每个信息比特对应的能量经过了编码冗余的摊薄。BPSK下每个调制符号携带1个码字比特但最终有意义的是信息比特。所以信噪比先用码率R折算再得到噪声方差。如果码率是1/2同样的Eb/N0下噪声功率要比未编码系统大一倍。这个换算出错的表现是曲线整体向右偏移看起来像性能损失了3dB。我排查过好几次这类问题最后发现都是公式里忘了乘码率。4.3 误码率曲线怎么解读跑出来的BER曲线横轴是Eb/N0纵轴是误码率用对数坐标。典型的LDPC曲线有三个区域低信噪比下的平坦区中间信噪比下的瀑布区以及高信噪比下的错误平层区。瀑布区曲线陡降是编码增益发挥作用的地方。错误平层则是译码器在某些特定错误图样上无法收敛导致BER下降趋势减缓。你判断一个LDPC码性能好不好通常看两个东西瀑布区出现的信噪比位置以及错误平层的高度。Codeword 512码长、列重3/行重6的规则码瀑布区在1dB附近和香农限大约差1.5-2dB这是正常水平。如果瀑布区出现在0.5dB以下那说明你的码构造或者译码实现很不错。如果你看到曲线在1.5dB之后突然变成垂直下降别高兴太早多半是帧数太少导致的统计假象。4.4 影响曲线好坏的几个因素迭代次数是最容易影响结果的参数。迭代次数太少译码器没有充分收敛曲线在瀑布区就会变得平缓像“拉伸”了一样。我曾经为了赶时间把迭代次数降到10结果瀑布区右移了0.5dB后来确认是迭代不够不是码的问题。反过来迭代次数从50加到100曲线改善不明显仿真时间却翻倍所以50是一个性价比很高的折中点。H矩阵的环结构也会影响性能。长度为4的环会让BP译码的消息出现明显的相关性导致性能下降。所以用generateH构造矩阵之后我建议加一个简单的4环检测把存在4环的矩阵重新生成一次。这个检测的优化目标就是让矩阵的围长至少为6。5. 实战中的坑与排查记录5.1 列置换导致信息位乱序这个坑我敢说90%第一次实现近似下三角编码的人都会踩。H矩阵经过列置换后码字的比特位置已经变了。如果你在编码端直接拿着原始信息位去套公式算出来的码字放进译码器译码结果当然对不上。解决思路是始终记得编码端输出的码字要经过逆置换恢复顺序。更稳的做法是在预处理完成后立刻记录三个东西rowPerm、colPerm、gap。编码时先对信息位做置换得到编码码字后再逆置换。译码端始终使用原始H不参与置换。5.2 phi矩阵奇异怎么办φ矩阵奇异是编码实现里比较棘手的问题。φ D − E T⁻¹B如果它不可逆整个编码过程就没法继续。我的排查顺序是确认gap是否过小。gap太小时D矩阵维度小φ更容易奇异。换一组列置换重新做RU预处理。预处理本身有一定随机性不同的置换结果可能影响φ是否奇异。换一个H矩阵重新生成。如果尝试多组变换后φ仍然奇异说明当前H矩阵本身就不适合RU编码。我用的策略是在ru_preprocess里多跑几轮记录第一个使φ非奇异的置换。如果跑了20轮还是找不到就去调整H的构造参数。这个检查最好放在预处理阶段自动完成不要在编码阶段才报错。5.3 译码不收敛或者性能极差如果你发现译码器输出的误码率居高不下先不要怀疑编码先检查LLR符号。BPSK映射有两个习惯0映射1、1映射-1和反过来的版本。只要调制、LLR计算、硬判决三处的符号保持一致系统就能工作。但如果你在某个地方搞混了最常见的现象是BER在0.5附近徘徊仿佛完全没有编码增益。另一个常见问题是校验矩阵用了非稀疏存储。H如果被Matlab转换成了full矩阵BP译码的每次迭代都要遍历大量零元速度慢还能忍更严重的是消息更新时那些零元位置也参与了计算导致算法行为异常。务必保持H为sparse类型。5.4 仿真效率太低怎么办我在初版仿真里用纯循环写BP译码码长512每帧迭代50次跑一个信噪比点要一个多小时。后来优化到向量化同样参数几分钟就跑完。向量化的核心是把所有校验节点的更新同时计算用矩阵运算替代逐边循环。标准BP的tanh规则在Matlab里可以用矩阵操作先把变量节点消息按H的行列结构展开用prod和sum聚合。如果标准BP太慢可以先换Min-Sum译码。Min-Sum只是把校验节点更新里的tanh简化成取最小绝对值性能差距不到0.2dB但速度能快好几倍。调试阶段强烈建议先用Min-Sum最后再换标准BP出最终曲线。5.5 问题速查表现象可能原因解决方向编码输出不满足校验方程列置换顺序没恢复检查colPerm使用方式编码时报phi奇异H矩阵不适合当前gap增大gap、重做预处理BER始终0.5LLR符号不一致统一调制/软信息/判决符号仿真的太慢循环未向量化改用Min-Sum或向量化BP高信噪比曲线突然断崖式下降错误比特累积太少增加每点帧数或错误门限6. 个人心得与扩展方向6.1 仿真工程组织的心得LDPC仿真做到后面真正消耗时间的不是算法本身而是工程组织。我个人的习惯是把预处理、编码、译码三个部分拆成独立模块分别用单元测试验证。编码完成后先做一步强制性检查对随机信息位编码得到码字c验证mod(H*c, 2)是否全零。这个检查每次改动代码后都要跑一遍能挡住大部分低级错误。预处理结果是离线数据我强烈建议保存成文件而不是每次仿真重新生成。我试过在帧循环里调用预处理结果大部分时间都耗在矩阵变换上编码反而没占多少时间。记住预处理是设计阶段的工作编码器是运行阶段的工作两者要分开。6.2 可以继续深入的方向如果你要把这套东西继续做下去有几个方向值得试。第一个是结构化LDPC码例如准循环LDPC它的校验矩阵由循环移位单位阵构成编码和译码都可以高度并行是5G标准采用的码型。第二个是分层译码它把校验节点分组依次更新收敛速度大约能提升一倍。第三个是归一化Min-Sum业界硬件实现基本都用它的变体仿真里加上归一化因子以后性能损失很小但复杂度大幅下降。还有一件事当我做完近似下三角编码再回头看传统G矩阵编码最大的感受是信道编码的很多“经典做法”在理论上是优雅的实际工程里却不能直接用。近似下三角编码的核心不是发明了新的数学而是用工程思维重新组织已有的线性代数让复杂度落在可控的范围内。这个思路本身比任何一个具体公式都值钱。最后给个小技巧在Matlab里调试LDPC仿真时把命令行窗口的tic/toc打开每个信噪比点跑完都看一眼时间。如果某个点突然变得特别慢先不要急着加并行先检查是不是误码率太低导致while循环里的错误比特门限一直达不到而无限跑下去。给每点加一个最大帧数上限是保护仿真脚本最简单的一招。 
长河君 · 2026/10/14 2:56:42 · 阅读 1.2k LDPC码近似下三角编码原理与Matlab实现全流程解析
LDPC码的编码实现是很多刚开始做信道编码仿真的同学绕不过去的一个坎。LDPC码本身性能逼近香农极限在5G新空口、Wi-Fi 6、DVB-S2这些标准里都是主力信道编码方案。但真正动手写Matlab仿真的时候你会发现译码端调通还算容易编码端反而最卡人。直接用生成矩阵G去编码H矩阵稀疏的优势全部丢掉一个码长上千的矩阵乘法在普通笔记本上能跑到你怀疑人生。这篇文章就把我用Matlab实现LDPC码近似下三角编码的完整过程讲清楚包括原理推导、代码实现、误码率曲线绘制还有我在仿真里踩过的那些坑。适合正在做LDPC仿真、或者准备把LDPC写进自己通信系统里的朋友参考。1. LDPC编码到底难在哪里1.1 LDPC为什么在5G时代“翻红”LDPC码是上世纪60年代提出的之后沉寂了很久直到90年代末才被重新发掘。它的核心优势在于校验矩阵H是稀疏的配合置信传播BP迭代译码性能可以非常接近香农限。相比Turbo码LDPC的译码并行度更高不需要交织器错误平层也更低所以在5G增强移动宽带场景里被选为数据信道的编码方案。不过大家平时接触LDPC基本都是从译码开始的。译码器拿到对数似然比LLR在Tanner图上做消息传递迭代几次以后硬判决得到码字。这个过程只需要H矩阵参与运算H矩阵的稀疏性让每次迭代的消息更新都只涉及少量边复杂度可控。但编码端不一样。编码的本质是要找到满足Hc^T 0的码字c这意味着你得先有办法从信息位s构造出校验位p。这就是LDPC编码比译码难的地方。1.2 从H到G是一条“死路”最直观的编码方式是把校验矩阵H通过高斯消元化成系统形式 [I P] 或者 [P I]然后得到生成矩阵G编码时直接算 c sG。这在理论上是完全没有问题的但实际做仿真的时候这一套基本跑不起来。原因很简单你的LDPC码为了保证性能H矩阵设计成稀疏的但H经高斯消元化成系统形式以后得到的 [I P] 中P部分往往是稠密的。稠密矩阵乘法的复杂度是O(n²)n是码长。算一笔账码长n 2048码率1/2一次编码需要做大约1000×1000量级的乘加运算。仿真时要跑成千上万帧还要在多个信噪比点下重复这个计算量是灾难性的。要是码长再上到4096或者8192普通的Matlab脚本根本没法在合理时间内完成。很多人第一次做LDPC仿真就卡在这里译码器已经写好了BER曲线却跑不出来因为编码端太慢。我当时也是这样后来才意识到必须绕开生成矩阵G直接从稀疏的H矩阵出发做编码。1.3 从H矩阵直接编码的思路既然不能显式构造G那就想办法直接在H上做文章。校验方程 Hc^T 0 等价于让码字分成信息位和校验位两组然后通过H的分块结构解出校验位。关键在于我们希望这个过程借助H的稀疏结构以接近线性的复杂度完成。如果H能写成一个近似下三角的形式例如H [ A B T ; C D E ]其中T是一个下三角矩阵那么编码过程就可以先从信息位推算出第一段校验位再用回代方式解出第二段校验位。这个过程的复杂度主要取决于T的大小和矩阵运算方式整体可以控制在O(n)到O(n²)之间远好于直接乘稠密G。这就是近似下三角编码方法的基本出发点。2. 近似下三角编码的原理与推导2.1 什么是近似下三角形式近似下三角编码也常被称作RU编码方法核心思路不是把H完全化成系统形式而是只做有限的行列置换让H变成一个带“缺口”的下三角结构。具体来说经过行列置换H被分成六块H [ A B T ; C D E ]这里矩阵的总列数是n总行数是m n - k。T是一个(m - g) × (m - g)的下三角矩阵g称为gap缺口大小表示右下角没有完全下三角化的列数。A的大小是(m - g) × kB是(m - g) × gC是g × kD是g × gE是g × (m - g)。如果g 0T本身是满秩的下三角矩阵编码会非常干净。但实际中随机生成的稀疏LDPC校验矩阵经过简单的行列置换不一定能让g 0通常需要允许一个小的gap比如g 16或者32。这个gap越大编码复杂度越高所以预处理的目标就是尽量让gap变小。2.2 关键公式推导从分块方程解出校验位设码字向量c分为三段c [s, p1, p2]其中s是信息位长度为kp1是第一部分校验位长度为gp2是第二部分校验位长度为m - g。把H和c代到校验方程Hc^T 0里按分块展开得到两个方程A s^T B p1^T T p2^T 0 C s^T D p1^T E p2^T 0先看第一个方程。因为T是下三角矩阵所以可以从下往上做回代用s和p1表达p2p2^T −T^{-1} (A s^T B p1^T)注意这里并不会真的计算T的逆矩阵而是利用T的下三角结构做前向代入复杂度只有O((m - g)²)量级。把这个p2代回第二个方程整理后得到(C − E T^{-1} A) s^T (D − E T^{-1} B) p1^T 0令 φ D − E T^{-1} B则有p1^T −φ^{-1} (C − E T^{-1} A) s^T也就是说只要φ矩阵可逆p1就能直接由信息位s算出来p1出来以后p2用一次回代就能解出。整个编码过程就完成了。2.3 gap值的选择与计算量影响从上面的公式可以看出gap值g直接决定复杂度。φ是一个g×g矩阵求φ的逆的复杂度是O(g³)。回代求解p2的复杂度是O((m - g)²)。整体编码复杂度大约是O(n g³)。如果g很小比如g 16那涉及φ的那些矩阵运算几乎可以忽略但如果g变得很大比如g 128φ的逆运算就会成为瓶颈。所以预处理阶段的核心目标就是在保持H稀疏性的前提下把g压缩到尽量小。我在实际仿真中一般接受g ≤ 32这个量级。对于码长256到2048的仿真这个规模下的编码时间完全可接受。2.4 为什么不用完全下三角你可能想问既然近似下三角这么绕直接想办法把H化成完全下三角不就行了吗理论上可以实际很难。完全下三角要求g 0这意味着校验矩阵本身就具备一种非常规律的结构。常规的规则LDPC码随机生成出来经过行列置换很难同时满足下三角和稀疏性两个约束。为了强行让g 0往往要做大量的行变换这个过程中H矩阵的稀疏性会被严重破坏到头来编码快了但矩阵变稠密存储和计算又变成新的瓶颈。近似下三角编码的意义就在这里稍微容忍一点点gap换来的代价其实很小但保留了H的稀疏性整体编码过程依然高效。用生活里的话说这是“花小钱办大事”的思路。3. Matlab代码实现把公式变成能跑的仿真3.1 生成稀疏校验矩阵H仿真第一步得先有一副“牌”可用。我常用的是规则码列重dv 3行重dc 6码率1/2。生成方式很多这里给一种最容易理解的下标随机交织构造法它生成的是一个行重列重都固定的正则H矩阵。function H generateH(n, dc, dv) % n: 码长 % dc: 行重, dv: 列重 % 返回 m x n 的 GF(2) 稀疏校验矩阵 m n * dv / dc; % 校验方程数 if mod(m, 1) ~ 0 error(参数必须满足 n*dv 能被 dc 整除); end H zeros(m, n); % 逐列放置1的位置保证每列有dv个1每行有dc个1 for j 1:n colIdx zeros(1, dv); for t 1:dv % 随机选择当前行重小于dc的行 candidates find(sum(H, 2) dc); colIdx(t) candidates(randi(length(candidates))); end H(colIdx, j) 1; end % 以防万一清理偶环简单策略接受现有矩阵后续译码验证 H sparse(H); end这段代码的思路是逐列确定1的位置每列随机放到行重还没满的行上。生成结果可能带有少量长度为4的环对于教学仿真影响不大。你如果追求更优的性能可以在这个基础上加环检测和重排逻辑后面再聊。3.2 预处理把H变成近似下三角结构接下来是核心的预处理环节。我直接给一个简化版的高斯消元找下三角块的方法。它的思想是选一部分列作为“主元列”对这些列做GF(2)高斯消元尽量把它们约化成下三角结构记录行置换和列置换gap就是消元过程中没能消成主元的那几列。function [Htilde, rowPerm, colPerm, gap] ru_preprocess(H) % RU预处理将H变换为近似下三角形式 % 返回 Htilde H(rowPerm, colPerm) % 使得 Htilde [A B T; C D E]T为下三角 H mod(H, 2); [m, n] size(H); rowPerm 1:m; colPerm 1:n; Htmp H; % 选主元列从右往左构造下三角 % 这里使用贪心策略优先选择当前可用列 pivCols []; % 主元列位置在当前Htmp中的列号 for pivotRow m:-1:1 % 在当前剩余矩阵中找一列在pivotRow行有1且在这行之前尽量简单 colCand find(Htmp(pivotRow, :) 1); if isempty(colCand) continue; end % 简单贪心取第一个候选列 selCol colCand(1); pivCols [pivCols, selCol]; % 消去其他行在该列上的1 rowsHaveOne find(Htmp(:, selCol) 1); rowsHaveOne(rowsHaveOne pivotRow) []; for r rowsHaveOne(:) Htmp(r, :) mod(Htmp(r, :) Htmp(pivotRow, :), 2); end end % 将主元列移到矩阵后端非主元列放到前端作为信息位和gap部分 nonPiv setdiff(1:n, pivCols); newColOrder [nonPiv, fliplr(pivCols)]; % pivCols从右下到右上fliplr保持顺序 colPerm newColOrder; % 行也重新排序将pivotRow序列恢复到自然顺序 [~, rowOrder] sort(pivCols); % 这个细节需要结合循环内记录优化 rowPerm rowOrder; Htilde H(rowPerm, colPerm); gap n - length(pivCols); % 非主元列数量这里作为近似gap end这段代码在我实际使用中做了一个关键简化它没有记录每步消元对应的行置换细节所以如果你直接复制运行行置换部分可能需要再处理一下。我更推荐的做法是把预处理视为“离线操作”在Matlab里跑一次把rowPerm、colPerm和gap保存成.mat文件后续每次仿真直接加载不必重复计算。我在仿真中就是这么组织的省了很多时间。真的要写一个严格正确的RU预处理核心是维护一个“未消元行”的集合每次选择一个合适的列把它和某一行交换位置然后对矩阵的其余部分做消元。因为涉及稀疏矩阵的频繁动态操作在Matlab里用全矩阵会更简单但码长超过2000以后会变慢。我觉得对于教学仿真直接预处理一次然后保存结果就够了没必要把预处理本身优化成生产级代码。3.3 编码函数信息位到完整码字有了近似下三角矩阵后编码函数可以直接按公式实现。这里要注意一个特别容易出错的地方预处理时对列做了置换等价于改变了码字比特的排列顺序。所以编码时信息位要先按照colPerm对应的顺序摆放编码完成后再把码字恢复到原始顺序。function [c, p1, p2] ldpc_encode_ru(s, H, rowPerm, colPerm, gap) % s: 原始信息位向量 (1 x k) % 返回完整的码字c原始比特顺序 [m, n] size(H); k n - m; % 根据colPerm把信息位放到置换后的对应位置 % 在列置换后前面部分是信息位 gap校验位后面部分是T对应的校验位 % 这里假设colPerm前k列为信息位后m列中gap个为p1其余为p2 Htilde H(rowPerm, colPerm); % 分割矩阵 A Htilde(1:end-gap, 1:k); B Htilde(1:end-gap, k1:kgap); T Htilde(1:end-gap, kgap1:end); C Htilde(end-gap1:end, 1:k); D Htilde(end-gap1:end, k1:kgap); E Htilde(end-gap1:end, kgap1:end); % 计算phi D - E * inv(T) * B实际用回代求解 X solve_tri(T, B); % X T^-1 B phi mod(D - E * X, 2); if rcond(double(phi)) 1e-10 error(phi矩阵接近奇异请重新生成H或者增大gap); end Y solve_tri(T, A); % Y T^-1 A % p1 -inv(phi) * (C * s - E * Y * s) s_perm s; % 简化处理信息位按原顺序实际应遵循colPerm p1 mod(-phi \ (mod(C * s_perm., 2) - E * mod(Y * s_perm., 2)), 2).; p1 mod(p1, 2); % p2 -inv(T) * (A * s B * p1) rhs mod(A * s_perm. B * p1., 2); p2 mod(-solve_tri(T, rhs), 2).; % 组装置换后的码字 c_perm [s_perm, p1, p2]; % 恢复原始顺序 c zeros(1, n); c(colPerm) c_perm; c mod(c, 2); end function X solve_tri(T, B) % 利用T的下三角结构解 X T^-1 B % 这里直接写满矩阵回代教学用 [r, c] size(T); X zeros(r, size(B, 2)); % 简单方式直接使用inv但代码可读性好 % 实际工程中应用向前代入 X T \ B; X mod(X, 2); end在上面的代码中关于信息位顺序的问题我做了转换在预处理时前三列已经假定是信息位集合。实际应用中你需要在预处理导出colPerm时就确定哪些列是信息位哪些列是校验位。一个比较方便的做法是选择一个初始列顺序让前k列是信息位后m列是校验位做完RU预处理后把前k列的原始列号存下来。这样编码时代码里s_perm s(:, colPerm(1:k))逻辑更清晰。为了避免混淆我这里省略了严格置换恢复但这部分值得你花时间仔细设计一次因为它是整个编码实现的“暗坑”之一。3.4 调制、AWGN信道与软信息计算编码完成后下一步是BPSK调制。我规定码字比特0映射到11映射到-1。这样接收端y x noise。噪声方差按照Eb/N0换算公式是sigma² 1 / (2 * R * 10^(EbN0/10))其中R k/n是码率。收到y后进入译码器的信道LLR初始值为Lch 2 * y / sigma²。这个公式的来源是BPSK下高斯信道LLR推导记得符号要对齐。function llr get_channel_llr(y, sigma2) llr 2 * y / sigma2; end3.5 LLR-BP译码器实现译码器我用置信传播的全信息版本。变量节点更新、校验节点更新轮流迭代。校验节点更新用tanh规则公式简洁但计算量偏大实际仿真时我喜欢先用Min-Sum快速验证链路然后再切回标准BP出性能曲线。下面给一个标准的对数域BP实现。function [decoded, iter] ldpc_bp_decode(llr, H, maxIter) % llr: 信道初始对数似然比 (1 x n) % H: 校验矩阵稀疏 [m, n] size(H); [rows, cols] find(H); % 校验矩阵非零位置 numEdges length(rows); % 初始化校验节点到变量节点的消息为0 msg_c2v zeros(1, numEdges); % 变量节点到校验节点消息初始为信道信息 msg_v2c llr; % 为方便消息传递建立索引映射表 % 这里用cell数组存每条边对应的变量节点和校验节点 v2cMsg zeros(numEdges, 1); for iter 1:maxIter % 校验节点更新 for e 1:numEdges % 实际Matlab应向量化教学演示用循环 varIdx cols(e); % 该边连接的变量节点 checkIdx rows(e); % 该边连接的校验节点 % 收集除当前边外所有邻居消息 neighbors find(rows checkIdx cols ~ varIdx); if isempty(neighbors) msg_c2v(e) llr(varIdx); continue; end prod_tanh 1; for nb neighbors(:) prod_tanh prod_tanh * tanh(v2cMsg(nb) / 2); end msg_c2v(e) 2 * atanh(prod_tanh); end % 变量节点更新 for e 1:numEdges varIdx cols(e); checkIdx rows(e); neighbors find(cols varIdx rows ~ checkIdx); msg_v2c_sum llr(varIdx); for nb neighbors(:) msg_v2c_sum msg_v2c_sum msg_c2v(nb); end v2cMsg(e) msg_v2c_sum; end % 计算总对数似然比并判决 totalLLR llr; for v 1:n edgesOfVar find(cols v); for e edgesOfVar(:) totalLLR(v) totalLLR(v) msg_c2v(e); end end decoded double(totalLLR 0); % LLR 0 判为比特1 if mod(H * decoded., 2) 0 return; end end end这里面的循环写法只是方便理解真正跑仿真的时候一定要做向量化否则码长512时迭代50次每一帧都会慢得让人抓狂。我在下面的“常见问题”里会专门讲如何把这段改成向量化实现。实际上你在实际调试时可以先跑通这个循环版再优化速度。3.6 BER主程序串联整个仿真链路最后是主程序。外层遍历信噪比内层逐帧统计错误比特。我习惯每个信噪比点累积到一定数量的错误比特再停止避免高信噪比时曲线出现明显的“抖动断层”。% 参数设置 n 512; k 256; m n - k; dc 6; dv 3; maxIter 50; EbN0_dB 0:0.5:3; numFramesMax 200; H generateH(n, dc, dv); % 假设置换已预处理完成这里省略ru_preprocess调用 % 实际仿真时从文件加载 RowPerm, ColPerm, gap for iSnr 1:length(EbN0_dB) EbN0 10^(EbN0_dB(iSnr)/10); sigma2 1 / (2 * k/n * EbN0); errorBits 0; totalBits 0; frame 0; while frame numFramesMax s randi([0,1], 1, k); [c, ~, ~] ldpc_encode_ru(s, H, rowPerm, colPerm, gap); % BPSK调制 x 1 - 2 * c; % AWGN信道 y x sqrt(sigma2) * randn(1, n); % 译码 llr get_channel_llr(y, sigma2); [decoded, ~] ldpc_bp_decode(llr, H, maxIter); errorBits errorBits sum(decoded ~ c); totalBits totalBits n; frame frame 1; if errorBits 200 break; end end ber(iSnr) errorBits / totalBits; end semilogy(EbN0_dB, ber, o-) xlabel(Eb/N0 (dB)) ylabel(BER) grid on注意这段代码把 ru_preprocess 的调用省略了你实际运行的时候需要先把 rowPerm、colPerm 和 gap 生成好。比较稳妥的做法是先跑一次预处理函数把它们保存下来[Htilde, rowPerm, colPerm, gap] ru_preprocess(H); save(ldpc_ru_params.mat, rowPerm, colPerm, gap);后续所有信噪比点的仿真都用同一份置换不要每次重新生成否则编码和译码的比特顺序对不上。4. 仿真参数设计与误码率曲线解读4.1 一套可复现的仿真参数我给出一套我实测过的参数表适合在没有专门硬件的情况下验证链路和出曲线。参数数值说明码长 n512太小则性能离理论差得远太大则仿真慢信息位长度 k256码率0.5列重 dv3变量节点度行重 dc6校验节点度最大迭代次数50典型值继续增大收益递减信噪比范围0-3 dB瀑布区主要分布范围每点最大帧数200或按错误比特数控制这套参数在普通笔记本上用向量化BP译码每个信噪比点大概几十秒到几分钟整体曲线一个小时内能跑完。如果你第一次调试建议先把码长降到256信噪比点数减少到4个先确认链路通不通。4.2 Eb/N0与噪声方差换算前面提到过换算公式这里再展开说一句。很多新手直接在信道上用固定的复数噪声功率结果发现BER曲线和理论值差得很远。原因在于LDPC码率不是1每个信息比特对应的能量经过了编码冗余的摊薄。BPSK下每个调制符号携带1个码字比特但最终有意义的是信息比特。所以信噪比先用码率R折算再得到噪声方差。如果码率是1/2同样的Eb/N0下噪声功率要比未编码系统大一倍。这个换算出错的表现是曲线整体向右偏移看起来像性能损失了3dB。我排查过好几次这类问题最后发现都是公式里忘了乘码率。4.3 误码率曲线怎么解读跑出来的BER曲线横轴是Eb/N0纵轴是误码率用对数坐标。典型的LDPC曲线有三个区域低信噪比下的平坦区中间信噪比下的瀑布区以及高信噪比下的错误平层区。瀑布区曲线陡降是编码增益发挥作用的地方。错误平层则是译码器在某些特定错误图样上无法收敛导致BER下降趋势减缓。你判断一个LDPC码性能好不好通常看两个东西瀑布区出现的信噪比位置以及错误平层的高度。Codeword 512码长、列重3/行重6的规则码瀑布区在1dB附近和香农限大约差1.5-2dB这是正常水平。如果瀑布区出现在0.5dB以下那说明你的码构造或者译码实现很不错。如果你看到曲线在1.5dB之后突然变成垂直下降别高兴太早多半是帧数太少导致的统计假象。4.4 影响曲线好坏的几个因素迭代次数是最容易影响结果的参数。迭代次数太少译码器没有充分收敛曲线在瀑布区就会变得平缓像“拉伸”了一样。我曾经为了赶时间把迭代次数降到10结果瀑布区右移了0.5dB后来确认是迭代不够不是码的问题。反过来迭代次数从50加到100曲线改善不明显仿真时间却翻倍所以50是一个性价比很高的折中点。H矩阵的环结构也会影响性能。长度为4的环会让BP译码的消息出现明显的相关性导致性能下降。所以用generateH构造矩阵之后我建议加一个简单的4环检测把存在4环的矩阵重新生成一次。这个检测的优化目标就是让矩阵的围长至少为6。5. 实战中的坑与排查记录5.1 列置换导致信息位乱序这个坑我敢说90%第一次实现近似下三角编码的人都会踩。H矩阵经过列置换后码字的比特位置已经变了。如果你在编码端直接拿着原始信息位去套公式算出来的码字放进译码器译码结果当然对不上。解决思路是始终记得编码端输出的码字要经过逆置换恢复顺序。更稳的做法是在预处理完成后立刻记录三个东西rowPerm、colPerm、gap。编码时先对信息位做置换得到编码码字后再逆置换。译码端始终使用原始H不参与置换。5.2 phi矩阵奇异怎么办φ矩阵奇异是编码实现里比较棘手的问题。φ D − E T⁻¹B如果它不可逆整个编码过程就没法继续。我的排查顺序是确认gap是否过小。gap太小时D矩阵维度小φ更容易奇异。换一组列置换重新做RU预处理。预处理本身有一定随机性不同的置换结果可能影响φ是否奇异。换一个H矩阵重新生成。如果尝试多组变换后φ仍然奇异说明当前H矩阵本身就不适合RU编码。我用的策略是在ru_preprocess里多跑几轮记录第一个使φ非奇异的置换。如果跑了20轮还是找不到就去调整H的构造参数。这个检查最好放在预处理阶段自动完成不要在编码阶段才报错。5.3 译码不收敛或者性能极差如果你发现译码器输出的误码率居高不下先不要怀疑编码先检查LLR符号。BPSK映射有两个习惯0映射1、1映射-1和反过来的版本。只要调制、LLR计算、硬判决三处的符号保持一致系统就能工作。但如果你在某个地方搞混了最常见的现象是BER在0.5附近徘徊仿佛完全没有编码增益。另一个常见问题是校验矩阵用了非稀疏存储。H如果被Matlab转换成了full矩阵BP译码的每次迭代都要遍历大量零元速度慢还能忍更严重的是消息更新时那些零元位置也参与了计算导致算法行为异常。务必保持H为sparse类型。5.4 仿真效率太低怎么办我在初版仿真里用纯循环写BP译码码长512每帧迭代50次跑一个信噪比点要一个多小时。后来优化到向量化同样参数几分钟就跑完。向量化的核心是把所有校验节点的更新同时计算用矩阵运算替代逐边循环。标准BP的tanh规则在Matlab里可以用矩阵操作先把变量节点消息按H的行列结构展开用prod和sum聚合。如果标准BP太慢可以先换Min-Sum译码。Min-Sum只是把校验节点更新里的tanh简化成取最小绝对值性能差距不到0.2dB但速度能快好几倍。调试阶段强烈建议先用Min-Sum最后再换标准BP出最终曲线。5.5 问题速查表现象可能原因解决方向编码输出不满足校验方程列置换顺序没恢复检查colPerm使用方式编码时报phi奇异H矩阵不适合当前gap增大gap、重做预处理BER始终0.5LLR符号不一致统一调制/软信息/判决符号仿真的太慢循环未向量化改用Min-Sum或向量化BP高信噪比曲线突然断崖式下降错误比特累积太少增加每点帧数或错误门限6. 个人心得与扩展方向6.1 仿真工程组织的心得LDPC仿真做到后面真正消耗时间的不是算法本身而是工程组织。我个人的习惯是把预处理、编码、译码三个部分拆成独立模块分别用单元测试验证。编码完成后先做一步强制性检查对随机信息位编码得到码字c验证mod(H*c, 2)是否全零。这个检查每次改动代码后都要跑一遍能挡住大部分低级错误。预处理结果是离线数据我强烈建议保存成文件而不是每次仿真重新生成。我试过在帧循环里调用预处理结果大部分时间都耗在矩阵变换上编码反而没占多少时间。记住预处理是设计阶段的工作编码器是运行阶段的工作两者要分开。6.2 可以继续深入的方向如果你要把这套东西继续做下去有几个方向值得试。第一个是结构化LDPC码例如准循环LDPC它的校验矩阵由循环移位单位阵构成编码和译码都可以高度并行是5G标准采用的码型。第二个是分层译码它把校验节点分组依次更新收敛速度大约能提升一倍。第三个是归一化Min-Sum业界硬件实现基本都用它的变体仿真里加上归一化因子以后性能损失很小但复杂度大幅下降。还有一件事当我做完近似下三角编码再回头看传统G矩阵编码最大的感受是信道编码的很多“经典做法”在理论上是优雅的实际工程里却不能直接用。近似下三角编码的核心不是发明了新的数学而是用工程思维重新组织已有的线性代数让复杂度落在可控的范围内。这个思路本身比任何一个具体公式都值钱。最后给个小技巧在Matlab里调试LDPC仿真时把命令行窗口的tic/toc打开每个信噪比点跑完都看一眼时间。如果某个点突然变得特别慢先不要急着加并行先检查是不是误码率太低导致while循环里的错误比特门限一直达不到而无限跑下去。给每点加一个最大帧数上限是保护仿真脚本最简单的一招。