简介面向通信工程专业学生与科研人员这份关于准循环低密度奇偶校验码的仿真源码围绕编译码误码率测试流程完整实现了从校验矩阵构造、生成矩阵转换、置信传播迭代译码到误码率统计与曲线绘制的全链路。代码按功能拆分为多个脚本包括主控程序、校验矩阵生成函数、生成矩阵转换函数、译码函数并附有数据文件便于直接运行与参数调整。整个压缩包共五个文件包含四个m脚本文件与一个mat数据文件大小仅五十一千字节轻量便携。目前已有二百七十二人学习下载。通过调整信道信噪比等参数可以直观观察不同条件下的误码率变化曲线深入理解准循环低密度奇偶校验码的构造特点与迭代译码性能适合作为课程设计、毕业设计或科研验证的参考实现。1. QC-LDPC 误码率仿真为什么值得把编译码链路完整搭一遍很多人以为 LDPC 仿真的难点在译码器但实际动手后才发现拉开差距的是校验矩阵构造和噪声方差设定。本资源是一套可直接跑的 MATLAB 误码率仿真源码覆盖 QC-LDPC 码的完整链路基矩阵扩展校验矩阵、GF(2) 高斯消元求生成矩阵、BPSK 调制加 AWGN、对数域置信传播BP译码最后汇总 BER/FER 曲线。QC-LDPC 是 5G NR 和 Wi-Fi 6 都在用的结构化 LDPC 码它的准循环特性让 MATLAB 仿真代码可以做到紧凑且可复现。适合两类人一类是做编码理论对比、需要快速出曲线的研究生另一类是正在做通信课程设计、需要完整可运行代码的从业者。文件拿到手就能跑重点是理解每个环节为什么这么写。2. 把基矩阵扩成校验矩阵QC-LDPC 构造的循环移位展开2.1 基矩阵与扩展因子 Z先看懂“图纸”QC-LDPC 的名字里“QC”就是 Quasi-Cyclic准循环。它的校验矩阵 H 不是随机稀疏矩阵而是由 mb 行、nb 列的块矩阵组成每个块是 Z×Z 的方阵。这个 mb×nb 的块矩阵叫基矩阵Base Matrix每个元素表示对应块的循环移位值取值为 -1 时代表全零块取值为 0 到 Z-1 时代表单位阵循环移位后的结果。码率与块数的关系是 R (nb - mb) / nb。比如 mb4、nb8码率就是 1/2。扩展因子 Z 决定码长n nb×Z信息位 k (nb-mb)×Z。这个设计的核心价值在于只要更换 Z 就能在相同基矩阵下得到不同码长而基矩阵本身决定了码的度分布和纠错性能硬件实现时也只需要存储基矩阵和移位值不需要存整张 H。2.2 扩展函数从基矩阵到完整校验矩阵function H qc_ldpc_parity_check(base_matrix, Z) % base_matrix: mb x nb 维基矩阵 % 元素的含义-1全零块s(0sZ)单位阵循环右移s位 % Z: 扩展因子每个块为 Z x Z 方阵 [mb, nb] size(base_matrix); H zeros(mb*Z, nb*Z); % 预分配 for r 1:mb for c 1:nb s base_matrix(r, c); if s 0 continue; % 全零块 end rows (r-1)*Z 1 : r*Z; cols (c-1)*Z 1 : c*Z; unit eye(Z); H(rows, cols) circshift(unit, [0, s]); end end end逻辑说明遍历基矩阵每个元素非负值就在对应位置放入一个循环移位单位阵。circshift(unit, [0, s])把单位阵的列循环右移 s 位行不动。这个“列循环移位”方向要统一不同论文可能用行移位代码里保持一致即可。参数说明Z 一般取 2 的幂或标准里指定的值比如 16、32、64。基矩阵里移位值的上限是 Z-1写代码的人不需要额外校验但构造基矩阵时要小心。我用 4×8 的基矩阵做演示Z16 时得到 64×128 的 H码率正好 0.5跑一圈仿真只需要几十秒非常适合调试整条链路。2.3 一个重要的构造约定后 mb 列做下三角托底随机生成基矩阵的坑在于扩展出的 H 在 GF(2) 上很可能秩亏导致后续求生成矩阵直接报错。我一般会让基矩阵的后 mb 个块列构成下三角结构对角块全为 0非对角块用随机移位值。这样后 mb×Z 列天然构成一个在 GF(2) 上可逆的方阵H 的行秩必然等于 mmb×Z从根上避免了秩亏问题。下三角托底还有个附带好处后 m 列天然对应校验位前 k 列对应信息位。这样求生成矩阵时只需要做行变换不需要做列置换记账复杂度低很多。% 4x8 基矩阵Z16后4列为下三角托底 base_matrix [ 0 1 2 -1 0 -1 -1 -1; -1 2 0 3 -1 0 -1 -1; 3 -1 1 -1 -1 -1 0 -1; -1 0 -1 2 -1 -1 -1 0 ];逻辑说明看后四列第 5 列第 1 行为 0第 6 列第 2 行为 0第 7 列第 3 行为 0第 8 列第 4 行为 0这就是下三角。前四列的移位值随意只要在合法范围内。参数说明这个基矩阵每行前四列有 2 到 3 个非负块度分布不算规则但对演示链路足够。想换成 Wi-Fi 或 5G 标准基矩阵时只需要保证“后 mb 列可逆”这一条否则就得动列置换。3. 编码器的两种落地姿势GF(2) 高斯消元与下三角化3.1 为什么编码要用生成矩阵 G 而不是直接用 H校验矩阵 H 满足 H·c^T0它描述的是码字约束不能直接拿它把信息比特映射成码字。编码需要生成矩阵 G满足 c u·G其中 u 是信息比特。LDPC 码的常规做法是先把 H 通过高斯消元化成系统形式 [P | I_m]然后生成矩阵就是 [I_k | P^T]。这里的算术全部在 GF(2) 上进行加法是异或。对 QC-LDPC 来说如果 Z 不大比如 16 或 32G 矩阵只有 k×n 大小直接乘就行。如果 Z 到了 256 以上G 是稠密矩阵存储和乘法开销暴涨就要换下三角化编码后面讨论。3.2 GF(2) 高斯消元的可跑实现function [G, H_sys] generator_from_H(H) % 在 GF(2) 上对 H 做初等行变换化成 [P | I_m] 系统形式 % 约定H 的后 m 列必须可逆基矩阵下三角托底保证 [m, n] size(H); Hw double(H); % 逐列处理目标是把右侧 m 列化成单位阵 for col 1:m pcol n - m col; % 当前行对应的校验列 pivot find(Hw(col:end, pcol), 1) col - 1; if isempty(pivot) error(右侧 m 列不可逆请检查基矩阵构造); end if pivot ~ col Hw([col pivot], :) Hw([pivot col], :); end for r 1:m if r ~ col Hw(r, pcol) 1 Hw(r, :) mod(Hw(r, :) Hw(col, :), 2); end end end if ~isequal(Hw(:, n-m1:end), eye(m)) error(消元结果异常右侧未化为单位阵); end H_sys Hw; P Hw(:, 1:n-m); G [eye(n-m), P.]; G mod(G, 2); end逻辑说明外层循环逐列找主元找到后把含主元的行换到当前行再用异或消除其他行的同列元素。与普通浮点高斯消元的差异在于两点一是没有除法主元只能是 1二是消元用mod(AB,2)等价于 XOR。消元结束后右侧 m 列应该正好变成单位阵这个isequal断言能拦住绝大多数构造错误。参数说明Hw double(H)这一步很关键。如果 H 是 uint8 存进来的MATLAB 的mod和比较运算会产生类型混乱后面避坑章会细说。G 返回的是 double 类型的 0/1 矩阵编码时mod(u*G, 2)直接得到码字。3.3 大 Z 时的后备方案近似下三角编码的思路当 Z 超过 128G 矩阵稠密得让人后悔。业界标准的处理是 Richardson-Urbanke 近似下三角编码把 H 通过行列置换化成 [A B; C D] 的块下三角形式其中 T 部分保持下三角。编码时先算校验位中间变量再回代复杂度从 O(n²) 降到 O(n g²)其中 g 是下三角化后的间隙。这套源码包里并没有把完整实现塞进来因为对 Z≤64 的仿真场景G 矩阵法的内存开销完全可以接受。如果你要跑标准码长 n1944才需要考虑这个方案。我一般处理方式是小 Z 用 G 矩阵法验证算法正确性验证完换大参数时再单独写编码器两条路互不干扰。4. 对数域 BP 译码从消息传递公式到可跑代码4.1 置信传播的直觉变量节点和校验节点互相投票LDPC 译码的本质是稀疏图上迭代。变量节点代表码字比特校验节点代表 H 的行约束。每一轮迭代中变量节点把当前软信息发给相邻校验节点校验节点根据“校验和必须为 0”的约束更新消息再传回。对数域表示下信道初始信息是 Lch 2y/σ²校验节点更新用双曲正切乘积公式Lcv 2·atanh(∏tanh(Lvc/2))。这个公式来自概率域的乘法在对数域变成 tanh 的乘法。变量节点更新更简单把某个校验节点传回的消息从自己的总后验里减掉避免消息回流放大这就是外信息。迭代若干轮后对后验 LLR 做硬判决LLR0 判为 1LLR0 判为 0。4.2 对数域 BP 译码器实现function [bits, iter_used] bp_decode_logdomain(y, H, max_iter, sigma2) % y: 信道观测值BPSK 映射0-1, 1--1 % H: 校验矩阵double % max_iter: 最大迭代次数 % sigma2: 噪声方差注意是方差不是标准差 [m, n] size(H); Lch 2 .* y ./ sigma2; Lvc repmat(Lch, m, 1); % 变量-校验消息初始为信道消息 Lcv zeros(m, n); % 校验-变量消息 Lv Lch; % 变量节点后验 LLR for it 1:max_iter % 校验节点更新 for i 1:m idx find(H(i, :)); if isempty(idx), continue; end for j idx others idx(idx ~ j); temp prod(tanh(0.5 .* Lvc(i, others))); % 数值稳定atanh 输入必须严格在 (-1,1) 内 temp max(min(temp, 1 - 1e-12), -1 1e-12); Lcv(i, j) 2 * atanh(temp); end end % 变量节点更新 for j 1:n idx find(H(:, j)); if isempty(idx), continue; end Lv(j) Lch(j) sum(Lcv(idx, j)); for i idx Lvc(i, j) Lv(j) - Lcv(i, j); end end % 硬判决 提前停止 bits double(Lv 0); if mod(H * bits(:), 2) 0 iter_used it; return; end end iter_used max_iter; end逻辑说明消息矩阵是 m×n 的稀疏结构只在 H 非零位置有意义。校验节点更新里的others表示除目标节点外的所有邻居乘积里必须排除自己这一路消息。变量节点更新里的Lv(j) - Lcv(i,j)就是外信息防止自己刚发出去的消息被自己收回来放大。提前停止条件是校验方程成立等价于找到一个有效码字。参数说明sigma2是噪声方差不是标准差。BPSK 解调软信息公式里 Lch 2y/σ² 对应的是方差。还有1 - 1e-12这个 clamp 很重要高信噪比下 tanh 乘积容易超过 0.9999999999不加这行直接atanh(1)Inf一轮迭代后整条消息矩阵全是 NaN。4.3 蒙特卡洛误码率主脚本与参数设定clear; clc; Z 16; base_matrix [ 0 1 2 -1 0 -1 -1 -1; -1 2 0 3 -1 0 -1 -1; 3 -1 1 -1 -1 -1 0 -1; -1 0 -1 2 -1 -1 -1 0 ]; H qc_ldpc_parity_check(base_matrix, Z); [m, n] size(H); k n - m; R k / n; [G, ~] generator_from_H(H); EbN0_dB 0:0.5:4; num_frames 2000; max_iter 50; ber zeros(size(EbN0_dB)); fer zeros(size(EbN0_dB)); for idx 1:length(EbN0_dB) EbN0 10^(EbN0_dB(idx)/10); sigma2 1 / (2 * R * EbN0); % BPSK 的噪声方差必须除以码率 n_bit_err 0; n_frm_err 0; frames_done 0; for f 1:num_frames u randi([0 1], 1, k); c mod(u * G, 2); x 1 - 2*c; y x sqrt(sigma2) * randn(1, n); [d, ~] bp_decode_logdomain(y, H, max_iter, sigma2); n_bit_err n_bit_err sum(d ~ c); n_frm_err n_frm_err any(d ~ c); frames_done f; if n_bit_err 100 break; % 误码数足够即停止节省仿真时间 end end ber(idx) n_bit_err / (frames_done * k); fer(idx) n_frm_err / frames_done; end semilogy(EbN0_dB, ber, o-, EbN0_dB, fer, s-); grid on; xlabel(Eb/N0 (dB)); ylabel(BER / FER); legend(BER, FER);逻辑说明主循环里信号映射是 0→1、1→-1加噪声后直接送译码器。每帧译码后对比码字c与判决结果d分别累计误比特数和误帧数。frames_done记录实际完成的帧数因为中途可能触发停止条件分母不能用num_frames。参数说明噪声方差的公式是 σ² 1/(2·R·10^(EbN0/10))。很多人直接用 σ² 10^(-EbN0/10)那条曲线会整体右移 0.5 到 1 dB看起来像译码器有问题实际上是能量定义错了。BPSK 信号能量 Es1EbEs/RN02σ²联立解出上面的式子。5. 误码率仿真避坑五个让我重跑一整天的常见问题5.1 NaN 传播导致判决全错现象高信噪比Eb/N0≥3dB下译码器输出的bits大量为 0BER 反而比低信噪比还高。原因校验节点更新里的atanh输入逼近 1MATLAB 返回 Inf下一轮变量更新变成 InfInf - Inf直接产生 NaN之后整个消息矩阵被污染硬判决全错。解决在temp prod(tanh(...))之后立刻 clamp 到 ±(1-1e-12)。这是标准 log-BP 实现里必须有的保护不是可有可无。浮点数的精度通常到 1e-16clamp 到 1e-12 既不会显著影响迭代行为又完全避开 Inf。5.2 曲线右移 1 dB噪声方差漏乘码率现象仿真的 BER 曲线与参考论文对比时整体右移约 1 dB但瀑布区斜率看起来没问题。原因能量定义搞混了。AWGN 信道下 BPSK 的 Es/N0 与 Eb/N0 之间差一个码率 Rσ²N0/2所以 σ² 1/(2·R·(Eb/N0))。漏了 R等效于把 Eb/N0 的高估了 1/R 倍1/2 码率时就是 3 dB 的偏差。解决主脚本里显式写出sigma2 1/(2*R*EbN0)并加上注释。我自己的习惯是每换一个码率就打印一次 sigma2 核对确认随码率变化趋势正确后再开跑。5.3 G 矩阵求不出来H 在 GF(2) 上秩亏现象调用generator_from_H时直接报 error提示找不到主元。原因随机生成的基矩阵扩展后的 H 在校验位部分不可逆。特别常见于从论文里抄的基矩阵——论文通常只保证码的最小距离和度分布不保证系统形式化方便。解决按本文 2.3 节的约定把基矩阵的后 mb 列设计为下三角托底结构从构造上保证后 m 列可逆。替换标准基矩阵时先跑一行assert(rank(double(generator_from_H(H))) size(H,1))做回归不满足就换另一套标准基矩阵参数别硬杠。5.4 Min-Sum 近似导致瀑布区后移现象用了 min-sum 近似后BER 曲线比 log-BP 后移约 0.3 dB 到 0.5 dB高信噪比平台更明显。原因min-sum 用最小绝对值代替 tanh 乘积本质上高估了校验消息的幅度。这种高估在低信噪比时被噪声掩盖但在瀑布区会显著恶化判决质量。解决做精确曲线对比用 log-BP需要提速时用归一化 min-sum给校验消息整体乘一个缩放因子 0.75 到 0.8。这个因子理论上跟码率和度分布有关工程上直接 0.75 起步看曲线微调。5.5 uint8 矩阵引发的类型暗雷现象H 用zeros(...,uint8)存储加速校验isequal(mod(H*c.,2), zeros(m,1))恒为 false但手动打印结果明明是 0。原因uint8 与 double 混算时MATLAB 会做隐式转换mod()的结果类型可能不是 double。isequal(uint8(0), double(0))返回 false因为它同时比较类型和值。解决统一用 double 存储 H或者判断改用all(mod(H*c.,2) 0)。前者一劳永逸后者省内存但不省心。仿真代码里 H 矩阵通常只有几万到几十万个元素double 完全扛得住没必要为省这点内存踩类型坑。6. 进阶验证与提速自检脚本、Min-Sum 与迭代次数调参6.1 自检脚本十分钟确认整条链路没白写在跑任何统计仿真之前先跑这个无噪声自检能把编码器和译码器里的逻辑错误一次性暴露出来。Z 16; base_matrix [ 0 1 2 -1 0 -1 -1 -1; -1 2 0 3 -1 0 -1 -1; 3 -1 1 -1 -1 -1 0 -1; -1 0 -1 2 -1 -1 -1 0 ]; H qc_ldpc_parity_check(base_matrix, Z); [m, n] size(H); k n - m; [G, ~] generator_from_H(H); assert(isequal(mod(H * G., 2), zeros(m, k))); % 1. 正交性 u randi([0 1], 1, k); c mod(u * G, 2); assert(isequal(mod(H * c., 2), zeros(m, 1))); % 2. 码字满足校验 y 1 - 2*c; % 无噪声 [d, it] bp_decode_logdomain(y, H, 50, 1e-6); assert(isequal(d, c)); % 3. 译码收敛 fprintf(自检通过译码迭代 %d 次\n, it);三个断言分别验证 H 与 G 正交、码字满足校验方程、译码器能收敛到原码字。这个习惯后来救了我很多次每次换了新基矩阵我都先跑一遍这个脚本再上蒙特卡洛。6.2 换 Min-Sum 提速算法对比与取舍log-BP 的校验节点更新有 tanh 和 atanh 运算循环慢。工程上更常用 min-sum校验消息的幅度取邻居消息的最小绝对值符号取邻居符号的乘积。这样整轮校验更新可以向量化同样帧数下提速 5 到 10 倍。算法每帧耗时Z16, 50次迭代BER 精度log-BP基准精确min-sum约 1/5瀑布区后移约 0.3 dB归一化 min-sum约 1/5接近 log-BP缩放 0.75~0.8我调试时的习惯是先用归一化 min-sum 加 10 次迭代抓链路 bug参数设大之后才换回 log-BP 跑正式曲线。两条路并存的价值就在这。6.3 最大迭代次数与平均迭代次数一个常被忽略的玄学参数最大迭代次数设得越大BER 曲线瀑布区越靠左但 50 次以上增益就很小了。更值得关注的是平均迭代次数你会发现它随 Eb/N0 升高而快速下降因为校验方程提前满足的概率变大。如果某个信噪比点平均迭代次数接近 max_iter说明那里还在瀑布区边缘曲线不稳。把 max_iter 从 20 提高到 50再观察平均迭代次数是否饱和能直接判断当前区域的运行状态。从那以后我每次跑误码率之前都强制走一遍自检脚本、确认 sigma2 公式、再开蒙特卡洛这套流程少走了太多弯路。希望帮到你。本文还有配套的精品资源点击获取