从零实现QC-LDPC编译码：MATLAB仿真与置信传播算法详解-编程实验室

简介：本资源是一份面向通信工程专业学生与研究人员的QC-LDPC码MATLAB误码率仿真完整实现，聚焦于准循环低密度奇偶校验码的编码、AWGN信道传输、BP迭代译码及BER性能评估等核心环节，解决LDPC码理论学习与工程仿真脱节问题。压缩包共5个文件（4个.m脚本+1个.mat数据），总大小51KB，其中m文件分别承担系统主流程控制（tops.m）、QC校验矩阵构造（func_QC_H.m）、H矩阵转G矩阵（func_H2G.m）、LDPC译码核心算法（func_Ldpc_dec.m）等关键功能，结构清晰、模块解耦，便于理解算法逻辑与调试优化。已有272人学习下载，读者可直接运行复现典型SNR下的误码率曲线，掌握QC-LDPC从矩阵构造、编码实现到迭代译码的全流程MATLAB建模方法，并借助内置信道建模与误比特统计逻辑快速开展参数对比实验。

1. 项目背景与核心价值

最近在通信算法圈子里，和几个朋友聊起信道编码的仿真实现，发现一个挺有意思的现象：很多人一提到LDPC码，第一反应就是去找现成的工具箱，比如MATLAB自带的comm.LDPCEncoder和comm.LDPCDecoder，或者网上一些封装好的函数。用起来确实方便，几行代码就能跑出误码率曲线。但问题也来了，一旦遇到性能不理想，或者想针对特定结构（比如QC-LDPC）做点优化，就完全不知道从何下手，成了一个“调包侠”。这让我想起几年前自己刚接触LDPC时，也是这么过来的，直到后来硬着头皮啃协议、手写编译码器，才真正搞明白那些参数和迭代过程背后的门道。

所以，今天我想分享的，就是围绕“QC-LDPC的编译码的matlab误码率仿真-源码”这个主题，进行一次从理论到代码的彻底拆解。QC-LDPC，也就是准循环低密度奇偶校验码，现在是5G NR数据信道和Wi-Fi 6/7等主流标准的宠儿，它的核心优势在于其校验矩阵具有规则的结构，能用非常简洁的方式描述一个巨大的矩阵，从而让编码变得高效，硬件实现也更友好。但光知道这些概念没用，你得能把它“跑”起来，能看到在不同的信噪比下，它的误码率曲线到底长什么样，这才是工程师该干的事。

这个项目的价值，绝不仅仅是给你一段能运行的MATLAB代码。它的核心在于，通过从零构建一个完整的QC-LDPC仿真链路，让你亲手触摸到编码、调制、过信道、解调、译码每一个环节的细节。你会清楚地知道，校验矩阵H是怎么从基矩阵和循环移位系数扩展来的；编码器是如何利用QC结构实现线性复杂度的；更重要的是，译码器里那个核心的置信传播（BP）算法，消息到底是怎么在变量节点和校验节点之间“流动”并最终纠正错误的。当你自己把这条链路搭通，看到仿真曲线一点点趋近理论值的时候，那种对算法“掌控感”的提升，是任何现成工具箱都给不了的。这份源码，就是一个带你穿透迷雾、直达核心的脚手架。

2. QC-LDPC的核心结构：从基矩阵到校验矩阵

要仿真QC-LDPC，第一步绝不是打开MATLAB就开始写for循环，而是必须彻底理解它的数据结构。QC-LDPC的巧妙之处，在于它用一个小巧的“模板”，通过“复制”和“移位”操作，生成了一个巨大的校验矩阵。

2.1 基矩阵与扩展因子：构建巨人的蓝图

我们可以把最终的校验矩阵H想象成一栋宏伟的建筑。而建造它，我们不需要画出每一块砖的位置，只需要一张设计蓝图和一个复制模组。这张蓝图就是基矩阵Hb，而复制模组的大小就是扩展因子Z。

基矩阵Hb是一个mb x nb的矩阵。但注意，它的元素不再是传统的0或1，而是整数。这个整数代表了“循环移位”的值。通常，-1代表一个Z x Z的全零矩阵，而0, 1, 2, ..., Z-1则代表一个Z x Z的单位矩阵，并向右循环移位相应的位数。

举个例子，假设我们有一个简单的基矩阵和扩展因子Z=3：

Hb = [ 0, -1, 2; -1, 1, 0 ]

那么，根据规则，我们将生成最终的校验矩阵H：

Hb(1,1)=0-> 替换为一个3x3的单位矩阵I。
Hb(1,2)=-1-> 替换为一个3x3的全零矩阵0。
Hb(1,3)=2-> 替换为单位矩阵向右循环移位2位得到的矩阵。
以此类推...

最终，H的维度将是(mb*Z) x (nb*Z)。在这个例子里，就是6 x 9。通过这种方式，我们仅用mb*nb个整数，就描述了一个可能非常稀疏（低密度）的大矩阵。这种描述方式在标准协议（如5G NR）中极为常见，协议文档里直接给出的就是基矩阵和移位值表。

在MATLAB中构建这个矩阵，是仿真的基石。下面是一个核心的构建函数片段，它清晰地展示了这个替换过程：

function H = expandHMatrix(baseMatrix, Z) % 根据基矩阵和扩展因子Z，生成QC-LDPC校验矩阵H % baseMatrix: mb x nb的基矩阵，元素为整数，-1表示全零矩阵 % Z: 扩展因子 % H: 生成的 (mb*Z) x (nb*Z) 的二进制校验矩阵 [mb, nb] = size(baseMatrix); H = zeros(mb*Z, nb*Z); % 预分配最终矩阵 for i = 1:mb for j = 1:nb shift = baseMatrix(i, j); startRow = (i-1)*Z + 1; startCol = (j-1)*Z + 1; if shift == -1 % 放置ZxZ的全零子矩阵，什么都不用做，因为已初始化为0 continue; else % 放置循环移位后的单位矩阵 % 先生成一个ZxZ的单位矩阵 I = eye(Z); % 然后进行循环右移shift位 shiftedI = circshift(I, [0, shift]); % 将子矩阵填充到H的对应位置 H(startRow:startRow+Z-1, startCol:startCol+Z-1) = shiftedI; end end end % 确保H是逻辑矩阵或双精度矩阵，便于后续运算 H = logical(H); end

注意：在实际的5G NR等标准中，移位值可能需要进行模Z运算，并且单位矩阵的循环移位方向（左移还是右移）需要根据协议规定严格确认。上述代码是一个通用示例，实际使用时需对照标准文档调整circshift的参数。

2.2 校验矩阵的稀疏性与Tanner图表示

生成的大矩阵H是高度稀疏的，这正是“低密度”名称的由来。这种稀疏性直接关联到它的译码算法——置信传播（BP）算法。为了直观理解BP算法，我们需要引入Tanner图表示法。

Tanner图是一个二分图，包含两类节点：

变量节点：对应校验矩阵H的每一列，也就是编码后码字的一个比特。总共有N = nb*Z个。
校验节点：对应校验矩阵H的每一行，也就是一个校验方程。总共有M = mb*Z个。

如果H矩阵在第i行、第j列的元素是1，那么在Tanner图上，第i个校验节点和第j个变量节点之间就有一条边相连。

QC结构使得其对应的Tanner图也具有非常规则的结构，这带来了两大好处：

编码便利：可以利用稀疏矩阵的准循环特性，使用移位寄存器等硬件友好结构实现线性复杂度的编码。
译码并行：规则的图结构使得在硬件实现译码器时，多个节点或边的计算可以并行处理，极大提升吞吐量。

理解Tanner图至关重要，因为后续的译码算法本质上就是信息在这个图上沿着边进行迭代传递的过程。在写仿真代码时，我们虽然不直接画图，但脑子里必须时刻有这个图景：消息从变量节点发出，经过校验节点处理，再返回变量节点，如此循环。

3. QC-LDPC编码器的实现：利用结构优势

有了校验矩阵H，接下来就是编码。对于一般的LDPC码，编码需要解决线性方程组H * x^T = 0^T，其中x是待求的码字。直接高斯消元法复杂度是O(N^3)，对于长码（N可能几千）是不可接受的。幸运的是，QC-LDPC的结构允许我们进行高效编码。

3.1 校验矩阵的系统化与高效编码原理

我们的目标是将校验矩阵H通过行列变换，化为系统形式[P | I]，其中I是M x M的单位矩阵，P是M x K的矩阵（K = N - M是信息位长度）。这样，对于信息向量s（长度为K），码字就可以直接写成x = [s, p]，其中校验位p = P * s (mod 2)。

对于QC-LDPC，由于其分块循环特性，我们可以将大矩阵的运算转化为对小矩阵（Z x Z循环矩阵）的运算。一个Z x Z的循环矩阵完全由其第一行（或第一列）决定，乘以一个向量等价于该向量与生成序列的循环卷积。在二进制域上，这可以进一步简化为移位和异或操作。

在实际操作中，我们并不总是需要对H进行显式的系统化。一种更常用的方法是利用校验矩阵的双对角结构。许多标准设计的QC-LDPC基矩阵，其最后几列会构成一个双对角结构（类似于[I; I I; I I; ...]的块形式）。这种结构使得我们可以通过“前向代入”的方式，递归地计算出校验位。

假设我们将码字向量x和校验矩阵H按扩展因子Z进行分块：

x = [x1, x2, ..., x_nb] ，每个xi是长度为Z的向量。 H = [H1,1 H1,2 ... H1,nb; H2,1 H2,2 ... H2,nb; ... H_mb,1 H_mb,2 ... H_mb,nb]

其中每个H_i,j是一个Z x Z的循环矩阵或零矩阵。

编码方程H*x^T = 0就转化为一系列子向量的方程。利用双对角结构，我们可以先计算最后几个与校验位直接相关的块，再逐步向前推导。下面给出一个基于此思想的简化编码函数框架：

function codeword = qc_ldpc_encode(info_bits, baseMatrix, Z) % QC-LDPC 编码函数 % info_bits: 行向量，长度为 K = (nb - mb) * Z % baseMatrix: mb x nb 的基矩阵 % Z: 扩展因子 % codeword: 编码后的行向量，长度为 N = nb * Z [mb, nb] = size(baseMatrix); K = (nb - mb) * Z; % 信息位长度 N = nb * Z; % 码字长度 % 1. 将信息比特分块，每块Z比特 info_blocks = reshape(info_bits, Z, nb - mb).'; % 每行是一个信息块 % 2. 初始化码字块 (nb个块，每个块Z比特) cw_blocks = zeros(nb, Z); % 前(nb-mb)块是信息块，后mb块是校验块（待计算） cw_blocks(1:nb-mb, :) = info_blocks; % 3. 核心编码计算（这里需要根据具体基矩阵结构设计计算顺序） % 假设基矩阵的最后mb列具有双对角结构，便于递归计算。 % 这是一个高度简化的示意流程，实际代码需根据基矩阵非零位置具体推导。 for i = 1:mb % 计算第 (nb-mb+i) 个校验块 % 原理：对于第i个校验方程， sum_over_j( H_ij * cw_block_j^T ) = 0 % 因此， H_i, (nb-mb+i) * cw_block_(nb-mb+i)^T = sum_over_other_j( H_ij * cw_block_j^T ) % 由于H_i, (nb-mb+i) 是单位矩阵的移位，求逆很容易（就是反移位）。 % 我们需要先计算等式右边的累加和。 accum = zeros(Z, 1); % 累加器 for j = 1:nb shift = baseMatrix(i, j); if shift >= 0 && j ~= (nb-mb+i) % 忽略当前待求的校验块本身 % 获取对应的码字块（可能是信息块或已求出的校验块） block = cw_blocks(j, :).'; % 计算 H_ij * block: 即对block进行循环移位shift位 shifted_block = circshift(block, shift); accum = mod(accum + shifted_block, 2); end end % 现在， H_i, (nb-mb+i) * p_i = accum % 假设 H_i, (nb-mb+i) 是单位矩阵循环右移shift_p位，则 p_i = circshift(accum, -shift_p) shift_p = baseMatrix(i, nb-mb+i); if shift_p >= 0 p_i = circshift(accum, -shift_p); % 反移位得到校验块 cw_blocks(nb-mb+i, :) = p_i.'; else error('基矩阵校验部分不符合双对角或近似结构，编码失败。'); end end % 4. 将码字块拼接成最终码字向量 codeword = reshape(cw_blocks.', 1, N); end

实操心得：编码器的实现强烈依赖于基矩阵的具体结构。上述代码是一个通用性框架，展示了分块计算的思想。在实际项目中（例如实现5G NR LDPC），你需要严格按照协议文档附录中给出的基矩阵和编码流程来实现，那通常是一个优化过的、步骤固定的算法，会利用更多的结构特性来减少计算量。自己实现通用编码器是很好的学习过程，但产品化时务必遵循标准。

3.2 编码正确性验证

在完成编码器后，第一时间不是去跑仿真，而是要做单元测试：验证生成的码字是否满足所有校验方程。

% 验证编码正确性 H = expandHMatrix(baseMatrix, Z); % 使用2.1节的函数生成H info = randi([0,1], 1, K); cw = qc_ldpc_encode(info, baseMatrix, Z); syndrome = mod(H * cw.', 2); % 计算伴随式 if all(syndrome == 0) disp('编码正确！伴随式全为零。'); else error('编码错误！伴随式非零。'); end

这个步骤至关重要，它能确保你的编码器逻辑基础是牢固的，避免将编码错误带入后续的仿真，导致误码率结果完全失真。

4. 置信传播译码算法详解与MATLAB实现

编码只是发送端的故事，通信链路的核心挑战在接收端：如何从被噪声污染的信号中恢复出原始信息。这就是LDPC码大放异彩的地方，其译码算法——置信传播（BP），也称为和积算法（SPA），是它性能优越的关键。

4.1 从概率到对数似然比：算法的基础

接收机收到信号后，经解调（如BPSK），我们得到每个比特的初始软信息。对于加性高斯白噪声（AWGN）信道，BPSK调制（映射：0->+1， 1->-1），这个软信息通常用对数似然比来表示：

LLR = ln( P(bit=0 | received_sample) / P(bit=1 | received_sample) )

经过推导，对于AWGN信道，初始LLR可以简化为：

LLR_n = (2 / sigma^2) * y_n

其中，y_n是解调后的采样值（比如，+1或-1加上噪声），sigma^2是噪声方差。LLR的正负表示硬判决策（正为0，负为1），绝对值大小表示置信度。BP算法就是在Tanner图上，迭代地更新这些LLR值。

算法的核心是两类节点之间的消息传递：

变量节点到校验节点消息：变量节点收集来自信道和其他校验节点的信息，综合后发给相邻的校验节点。
校验节点到变量节点消息：校验节点根据与之相连的所有变量节点发来的消息，计算出一个新的、基于校验约束的信息，发回给各个变量节点。

4.2 对数域BP算法实现与优化

在概率域直接计算涉及大量乘法，容易下溢。因此工程上普遍采用对数域BP算法，将乘法变为加法，更加稳定。以下是核心更新公式的简化版（假设使用tanh运算的log域近似，即min-sum算法或其改进型，这是硬件实现的主流）：

设L(q_mn)是从变量节点n发送给校验节点m的消息，L(r_mn)是从校验节点m发送给变量节点n的消息。N(m)表示与校验节点m相连的所有变量节点集合，N(m)\n表示除去节点n。

校验节点更新（min-sum算法，计算L(r_mn)）：
```
L(r_mn) = ( Π_{n'∈N(m)\n} sign(L(q_mn')) ) * min_{n'∈N(m)\n} ( |L(q_mn')| )
```
这个公式非常直观：输出消息的符号是所有输入消息符号的乘积（奇偶校验的体现），幅度是输入消息幅度中最小的那个（最不可靠的链路决定整个校验的可靠性）。为了补偿min-sum算法的性能损失，通常会引入一个衰减因子α（如0.8），即L(r_mn) = α * L(r_mn)，这称为归一化min-sum算法。
变量节点更新（计算L(q_mn)）：
```
L(q_mn) = L_channel_n + Σ_{m'∈M(n)\m} L(r_m'n)
```
其中L_channel_n是信道的初始LLR，M(n)是与变量节点n相连的所有校验节点集合。变量节点发给某个校验节点的消息，等于初始信道信息加上来自其他所有校验节点消息的总和。
后验LLR计算与硬判决：在每次迭代后，计算每个比特的后验LLR用于硬判决：
```
L_total_n = L_channel_n + Σ_{m∈M(n)} L(r_mn)
```
硬判决：hat_x_n = 0 if L_total_n >= 0 else 1。

然后检查H * hat_x^T == 0，如果成立则译码成功，提前终止迭代。

下面是在MATLAB中实现归一化min-sum算法译码器的核心循环部分：

function [decoded_bits, iter_used] = qc_ldpc_decode(llr_received, H, max_iter, alpha) % QC-LDPC 译码函数 (归一化 min-sum 算法) % llr_received: 接收到的LLR值，行向量，长度N % H: 校验矩阵 (M x N 逻辑矩阵) % max_iter: 最大迭代次数 % alpha: 归一化因子 (通常0.6~0.9) % decoded_bits: 译码输出的硬比特 % iter_used: 实际使用的迭代次数 [M, N] = size(H); % 初始化变量节点到校验节点的消息 L(q_mn) L_q = zeros(M, N); % 初始化校验节点到变量节点的消息 L(r_mn) L_r = zeros(M, N); % 获取每个校验节点和变量节点的邻居信息，避免在循环中频繁调用find [cn_neighbors, vn_neighbors] = get_tanner_graph_neighbors(H); % 将信道LLR赋值给变量节点的初始值 L_channel = llr_received; % 迭代译码主循环 for iter = 1:max_iter % --- 步骤1: 变量节点更新，计算 L(q_mn) --- for n = 1:N % 找到与变量节点n相连的所有校验节点 connected_cns = vn_neighbors{n}; for m_idx = 1:length(connected_cns) m = connected_cns(m_idx); % L(q_mn) = L_channel_n + sum of L(r) from all OTHER check nodes sum_L_r = sum(L_r(connected_cns, n)) - L_r(m, n); L_q(m, n) = L_channel(n) + sum_L_r; end end % --- 步骤2: 校验节点更新，计算 L(r_mn) (归一化 min-sum) --- for m = 1:M % 找到与校验节点m相连的所有变量节点 connected_vns = cn_neighbors{m}; num_neighbors = length(connected_vns); % 计算所有L(q_mn)的符号和绝对值 signs = sign(L_q(m, connected_vns)); abs_vals = abs(L_q(m, connected_vns)); % 计算总符号（所有输入符号的乘积） total_sign = prod(signs); for n_idx = 1:num_neighbors n = connected_vns(n_idx); % 对于当前边(m,n)，找出除n外其他所有邻居的最小绝对值 other_abs = abs_vals; other_abs(n_idx) = []; % 移除当前节点n对应的值 min_abs = min(other_abs); % 计算当前边输出消息的符号（总符号除以当前输入符号） current_sign = total_sign * signs(n_idx); % 等价于 prod(signs(其他节点)) % 应用归一化min-sum: L(r_mn) = alpha * current_sign * min_abs L_r(m, n) = alpha * current_sign * min_abs; end end % --- 步骤3: 计算后验LLR并尝试硬判决 --- decoded_bits = zeros(1, N); for n = 1:N connected_cns = vn_neighbors{n}; L_total = L_channel(n) + sum(L_r(connected_cns, n)); decoded_bits(n) = L_total < 0; % LLR<0 判为1 end % --- 步骤4: 校验伴随式，若全零则提前终止 --- syndrome = mod(H * decoded_bits.', 2); if all(syndrome == 0) iter_used = iter; return; % 译码成功，退出 end end % 达到最大迭代次数仍未成功 iter_used = max_iter; % 返回最后一次迭代的硬判决结果 end function [cn_neighbors, vn_neighbors] = get_tanner_graph_neighbors(H) % 预处理函数：获取Tanner图中每个节点的邻居列表，大幅提升循环效率 [M, N] = size(H); cn_neighbors = cell(M, 1); vn_neighbors = cell(N, 1); for m = 1:M cn_neighbors{m} = find(H(m, :)); end for n = 1:N vn_neighbors{n} = find(H(:, n)); end end

踩坑实录与性能调优：
邻居预计算：在译码循环中，反复使用find(H(m,:))或find(H(:,n))查找非零元素是巨大的性能瓶颈。上述代码中的get_tanner_graph_neighbors函数在译码前一次性计算出所有邻居关系，是必须做的优化，对于长码，性能可能有数量级的提升。
归一化因子α：min-sum算法是BP算法的近似，会损失一些性能。引入归一化因子α（α<1）可以补偿一部分损失。这个值需要根据具体的码率和信噪比区间进行微调，通常在0.6到0.9之间。可以通过仿真，在目标误码率（如1e-5）附近寻找使性能最优的α值。
迭代终止：每次迭代后都进行伴随式校验是必要的，一旦成功立即退出，可以节省大量不必要的计算。特别是在高信噪比下，多数帧可能几次迭代就成功了。
矩阵运算 vs. 节点运算：上述实现是“节点为中心”的，逻辑清晰。对于追求极致速度的仿真，可以考虑“层调度”或“部分并行”的更新方式，甚至将消息用矩阵形式存储，利用MATLAB的向量化运算加速，但代码会复杂很多。建议初学者先从节点式实现开始，确保正确性。

5. 构建完整的误码率仿真链路

有了可靠的编译码器，我们就可以搭建一个完整的蒙特卡洛仿真系统，来评估QC-LDPC码在不同信噪比下的性能。仿真的目标是得到一条误码率（BER）和误帧率（FER）随信噪比（Eb/N0）变化的曲线。

5.1 仿真流程与参数设置

一个标准的仿真链路包含以下步骤：

生成随机信息比特。
QC-LDPC编码。
调制（如BPSK：0->+1, 1->-1）。
添加高斯白噪声。
解调并计算初始LLR。
QC-LDPC译码。
统计误比特数和误帧数。

仿真的核心参数包括：

信噪比范围：EbN0_dB，通常是一个向量，如0:0.5:3。
每个信噪比下的仿真帧数：num_frames。为了在低误码率（如1e-5）下得到统计可靠的结果，需要足够的错误事件。一个经验法则是，至少需要收集100个错误比特或50个错误帧，仿真才可停止。因此，在低信噪比区（高误码率）可以跑少一些帧数，在高信噪比区必须跑非常多的帧数。
最大迭代次数：max_iter，通常设为50或100。实际迭代次数由译码器提前终止条件控制。

下面是一个仿真主函数的框架：

function [ber, fer] = run_ldpc_simulation(baseMatrix, Z, EbN0_dB_list, max_iter, alpha, max_frames_per_snr, max_bit_errors) % 运行QC-LDPC误码率仿真 % 返回每个信噪比下的误码率(BER)和误帧率(FER) [mb, nb] = size(baseMatrix); K = (nb - mb) * Z; N = nb * Z; code_rate = K / N; % 计算码率 num_snr = length(EbN0_dB_list); ber = zeros(1, num_snr); fer = zeros(1, num_snr); % 预先展开H矩阵，并获取邻居信息（译码用） H = expandHMatrix(baseMatrix, Z); [cn_neighbors, vn_neighbors] = get_tanner_graph_neighbors(H); fprintf('开始仿真，码率=%.3f, 码长=%d, 信息位=%d\n', code_rate, N, K); for snr_idx = 1:num_snr EbN0_dB = EbN0_dB_list(snr_idx); % 将Eb/N0 (dB) 转换为信噪比 SNR per bit (线性值) EbN0 = 10^(EbN0_dB / 10); % 对于BPSK，符号能量Es = Eb * code_rate? 注意关系。 % 更直接的方式：计算噪声方差 sigma^2 % 对于BPSK，每个符号能量Es = 1。 SNR = Es / (N0) = 1 / (2 * sigma^2) % 而 Eb/N0 = (Es / code_rate) / N0 = SNR / code_rate % 所以 sigma^2 = 1 / (2 * SNR) = 1 / (2 * (EbN0 * code_rate)) sigma2 = 1 / (2 * EbN0 * code_rate); % 噪声方差 sigma = sqrt(sigma2); fprintf('处理 Eb/N0 = %.2f dB ... ', EbN0_dB); bit_errors = 0; frame_errors = 0; total_bits_simulated = 0; frames_simulated = 0; while (frames_simulated < max_frames_per_snr) && (bit_errors < max_bit_errors) % 1. 生成随机信息比特 info_bits = randi([0, 1], 1, K); % 2. 编码 tx_codeword = qc_ldpc_encode(info_bits, baseMatrix, Z); % 3. BPSK调制: 0 -> +1, 1 -> -1 tx_signal = 1 - 2 * tx_codeword; % 4. 添加AWGN噪声 noise = sigma * randn(1, N); rx_signal = tx_signal + noise; % 5. 计算初始LLR: LLR = (2 / sigma^2) * y llr_received = (2 / sigma2) * rx_signal; % 6. 译码 [decoded_bits, ~] = qc_ldpc_decode(llr_received, H, max_iter, alpha, cn_neighbors, vn_neighbors); % 7. 提取信息位进行比较（假设系统码，前K位为信息位） decoded_info = decoded_bits(1:K); % 统计错误 frame_error = any(decoded_info ~= info_bits); bit_error = sum(decoded_info ~= info_bits); frame_errors = frame_errors + frame_error; bit_errors = bit_errors + bit_error; total_bits_simulated = total_bits_simulated + K; frames_simulated = frames_simulated + 1; end ber(snr_idx) = bit_errors / total_bits_simulated; fer(snr_idx) = frame_errors / frames_simulated; fprintf('完成 %d 帧, BER=%.2e, FER=%.2e\n', frames_simulated, ber(snr_idx), fer(snr_idx)); end end

5.2 结果可视化与性能分析

仿真结束后，我们需要绘制BER/FER曲线图，这是评估编码性能最直观的方式。

% 假设仿真结果已存储在 ber, fer 变量中 figure; semilogy(EbN0_dB_list, ber, 'b-o', 'LineWidth', 1.5, 'MarkerFaceColor', 'b', 'DisplayName', 'BER'); hold on; semilogy(EbN0_dB_list, fer, 'r-s', 'LineWidth', 1.5, 'MarkerFaceColor', 'r', 'DisplayName', 'FER'); grid on; xlabel('Eb/N0 (dB)'); ylabel('误码率/误帧率'); title('QC-LDPC码性能仿真'); legend('Location', 'best'); % 可以添加未编码BPSK的性能曲线作为参考 EbN0_linear = 10.^(EbN0_dB_list/10); ber_uncoded = qfunc(sqrt(2*EbN0_linear)); % 未编码BPSK的理论BER semilogy(EbN0_dB_list, ber_uncoded, 'k--', 'DisplayName', '未编码BPSK'); hold off;

如何解读曲线：

瀑布区：在低信噪比时，曲线陡峭下降，这是LDPC码的“陡峭译码阈值”特性的体现，说明一旦信噪比超过某个门限，性能会急剧改善。
错误平层：在高信噪比时，曲线下降变得平缓。这通常是由于码字中存在的“陷阱集”或“停止集”造成的，是LDPC码结构的固有特性。通过优化校验矩阵设计（如避免短环）可以降低错误平层。
与未编码对比：可以看到，在相同Eb/N0下，编码后的BER远低于未编码，这就是编码增益。但注意，由于码率小于1，比较应在相同的信息比特能量上进行，即横坐标是Eb/N0。
BER vs FER：通常FER曲线在BER曲线上方。当FER很低时（如<1e-4），要获得可靠的FER统计需要仿真海量帧，计算代价很高。此时，BER可以作为辅助评估指标。

6. 源码工程化与高级话题探讨

完成基本仿真后，一个健壮的、可复用的源码工程还需要考虑更多细节。

6.1 代码模块化与性能优化

将代码拆分为独立、功能清晰的模块是专业工程的习惯：

generate_qc_ldpc_matrix.m: 根据协议生成或读取基矩阵，并扩展为H。
qc_ldpc_encoder.m: 编码器模块。
qc_ldpc_decoder.m: 译码器模块（支持不同的算法，如Log-BP, Min-Sum, Normalized Min-Sum）。
simulate_ber.m: 仿真主循环和链路。
plot_results.m: 绘图和结果分析。

性能优化技巧：

向量化：在可能的地方用矩阵运算代替for循环。例如，校验节点更新中的min和prod操作，可以对整个L_q矩阵的切片进行操作。
提前终止：如前所述，译码成功的帧立即停止迭代。
并行计算：如果拥有Parallel Computing Toolbox，可以使用parfor循环并行仿真不同信噪比点或不同数据帧，这是加速蒙特卡洛仿真最有效的手段。注意，每个并行 worker 需要独立的数据和函数副本。
使用更快的算法：Min-Sum比Log-BP快，但性能有损失。Offset Min-Sum是另一种更好的折中。

6.2 不同QC-LDPC码的尝试与标准实现

你可以尝试仿真不同码率、不同码长的QC-LDPC码。一个很好的起点是去实现5G NR 标准中定义的LDPC码。3GPP TS 38.212协议文档的附录中给出了完整的基矩阵定义（包括BG1和BG2，以及各种扩展因子Zc）。实现标准码的意义在于，你可以将自己的结果与公开发表的论文或报告进行对比，验证仿真链路的正确性。

实现5G NR LDPC编码器时需要注意，其编码过程有详细的步骤，包括信息比特的填充、CRC附加、基图选择、比特选择等，远比我们上面的通用示例复杂。译码器则可以相对通用。

6.3 常见问题排查与调试

如果你的仿真曲线出现异常，比如BER不随信噪比下降，或者比未编码性能还差，请按以下步骤排查：

编码验证：这是第一步，也是最重要的一步。确保你的编码器输出的每一个码字都满足H * cw^T = 0。随机测试成千上万个随机信息向量。
信道模型验证：确保你的噪声方差sigma^2计算正确。对于BPSK，关系是sigma^2 = 1/(2 * SNR) = 1/(2 * (EbN0 * code_rate))。一个简单的验证方法是，关闭编码和译码，直接比较发送的BPSK符号和接收符号的误符号率，看是否符合ber_uncoded = qfunc(sqrt(2*EbN0))的理论值。
LLR计算验证：在无噪声情况下（sigma极小），发送全零码字（BPSK映射为全+1），接收端LLR应该是一个很大的正数。加入噪声后，检查LLR的分布是否符合预期。
译码器单步调试：构造一个简单的、已知的错误图案。例如，发送全零码字，在信道中人为翻转几个比特。然后单步运行译码器，观察LLR消息的传递过程，看是否能纠正这些错误。检查校验节点和变量节点的更新公式是否实现正确。
矩阵稀疏性检查：确保生成的H矩阵是稀疏的。如果密度过高，不仅性能差，译码复杂度也会爆炸。
提前终止逻辑：检查伴随式校验H * decoded_bits'的计算是否正确，确保译码成功时能正确退出。