简介面向卫星导航与信道编码研究者的 MATLAB 仿真资源围绕基于信念传播BP算法的多进制低密度奇偶校验LDPC码展开可用于评估该编码在不同信噪比下的误码率表现并帮助理解 Tanner 图结构与消息传递解码原理。由于卫星导航信号易受电离层、多径等干扰该仿真代码特别适合用来验证多进制 LDPC 码的纠错能力与抗干扰增益压缩包共 4 个文件包含 1 个解码主程序与 3 个参数或数据文本整体大小仅 4KB结构紧凑便于快速载入 MATLAB 环境运行对比实验。目前已有 856 人下载学习适合正在研究 LDPC 编码、卫星导航抗干扰或信道编码仿真的学生与工程师参考。仿真代码提供 BP 算法迭代解码的完整实现通过修改迭代次数、信噪比等参数可直观观察多进制 LDPC 码的纠错表现txt 文本中的指数矩阵、输入数据和码字元素则为复现实验提供了必要支撑。由于体积轻量、文件不多也非常适合作为课程设计或毕业设计的入门样例便于在真实代码基础上做二次修改与性能分析。1. 一份MATLAB LDPC BP工程先确认你要解决的是二进制还是多进制的问题看到“matlab_ldpc64_BP”这种命名大概率是有人把LDPC码的BP译码仿真打包发了出来后缀里的“64”可能是码长、信息位长度或校验矩阵的维度“BP”则明确指向置信传播。这类工程在导航、深空通信里很常见因为LDPC已经进入DVB-S2、5G NR、北斗等标准。但很多人拿网上下载的LDPC仿真代码跑一遍发现误码率下不去或者多进制调制时性能反而比二进制差。问题通常不在算法本身而在BP译码的LLR计算、归一化因子、最大迭代次数和执行细节。这篇博文会从BP译码的原理讲起给出一个能跑的MATLAB最小实现然后说明多进制LDPC在导航场景里怎么做端到端仿真最后用几条判断规则帮你确认译码器对不对。如果你正卡在MATLAB里的LDPC代码上这篇文章能帮你少走弯路。2. LDPC的BP译码原理校验矩阵、Tanner图与消息迭代2.1 从校验矩阵到稀疏图LDPC的“低密度”体现在哪里LDPCLow-Density Parity-Check码由一个稀疏校验矩阵H定义。所谓“低密度”是指H中1的数量相对于矩阵规模非常少。考虑一个典型的(64, 32)LDPC码H是32×64的矩阵每行可能有6个1每列也可能只有3个11的密度低于0.1。正是这种稀疏性让基于图的消息传递算法成为可能。在MATLAB里生成一个准循环LDPC校验矩阵常用dvbs2ldpc、ldpcQuasiCyclicMatrix或直接使用H矩阵本身。例如% 构造一个简单的(6,3)规则LDPC校验矩阵仅用于教学演示 H [1 1 0 1 0 0; 0 1 1 0 1 0; 1 0 1 0 0 1];这个H每行有3个1每列有2个1属于规则LDPC。实际系统中H的维度远大于此但消息传递的图结构是一样的。把H中的每一行看作一个“校验节点”check node每一列看作一个“变量节点”variable node变量节点和校验节点之间的连线就是H中为1的位置这样形成的就是Tanner图。BP译码的本质就是在Tanner图上迭代交换概率消息。2.2 BP译码的消息类型与迭代框架BP译码有两种等价形式概率域SPASum-Product Algorithm和对数域LLR域。概率域里传递的是后验概率涉及大量乘法数值容易下溢实际代码几乎都用对数域。对数域BP中变量节点向校验节点传递的消息是LLR值定义为L log(P(x0)/P(x1))信道初始LLR可以用2*r/σ^2计算其中r是接收到的BPSK符号σ^2是噪声方差。一次完整的迭代分两步变量节点到校验节点的消息将当前变量的所有外部消息相加。校验节点到变量节点的消息使用tanh规则或最小和近似。最小和近似是工程中最常用的简化它把校验节点的更新从复杂的双曲函数运算改成取最小值和符号乘积性能损失通常小于0.2dB。MATLAB自带的ldpcDecode支持max-log近似就是最小和的另一种叫法。2.3 对数域BP用LLR简化乘法给出校验节点更新的标准公式设L_v2c[j]表示从变量节点j发给校验节点i的消息校验节点i收到所有邻居变量节点的消息后返回给变量节点j的消息为L_c2v[i][j] 2 * atanh( prod_k( tanh( L_v2c[i][k] / 2 ) ) )这里的k遍历除了j以外的所有邻居。用最小和近似后公式简化为L_c2v[i][j] ( prod_k sign( L_v2c[i][k] ) ) * min_k | L_v2c[i][k] |其中k同样要排除j本身。实现时通常先计算整行的绝对值和符号乘积再逐个更新避免每对一个邻居就重算一次。下面是一个校验节点更新的MATLAB片段function L_c2v checkNodeUpdate(L_v2c) % L_v2c: m x n 矩阵m为校验节点数n为变量节点数 % 假设H中缺失的位置为0实际计算只处理非零位置 [m, n] size(L_v2c); L_c2v zeros(m, n); for i 1:m % 取当前行所有非零消息 row L_v2c(i, :); abs_row abs(row); sign_row sign(row); % 计算整行的符号乘积和最小绝对值 total_sign prod(sign_row(row ~ 0)); min1 min(abs_row(abs_row 0)); % 对每个非零位置排除自身后计算 for j 1:n if row(j) 0 continue; end % 排除第j个元素后的符号乘积 if sign_row(j) 0 s total_sign; else s -total_sign; end % 排除自身后的最小绝对值 if abs_row(j) min1 % 如果最小值出现了多次找第二小值 mask (1:n ~ j) (abs_row 0); m2 min(abs_row(mask)); else m2 min1; end L_c2v(i, j) s * m2; end end end这个函数的行为是先按行计算整行的符号乘积和最小绝对值再对每个非零元素排除自身后重新计算。注意abs_row 0的判断是为了跳过矩阵中的0占位实际H矩阵的0位置在计算中必须被忽略。如果直接在循环里逐个比较性能会差很多但对于64码长的教学代码清晰比速度更重要。3. 用MATLAB从零写二进制LDPC的BP译码器3.1 构造一个可仿真的LDPC码要验证BP译码不能只用上面的(6,3)矩阵。这里构造一个(64,32)LDPC码码率1/2行重6列重3。最简单的生成方式是使用MATLAB的comm.LDPCEncoder但为了展示从零写BP我们先手工生成一个准循环矩阵。实际上MATLAB R2021b之后提供了ldpcQuasiCyclicMatrix函数可以直接生成5G NR标准中的校验矩阵。比如% 5G NR标准中一个码长约1944码率1/2的基矩阵 cfg ldpcQuasiCyclicMatrix(1/2, 1944, 5G); H cfg.ParityCheckMatrix;但为了对应标题中的“64”我们也可以自己构造一个64×128的规则矩阵。这里给出一种常见的构造思路用单位矩阵循环移位生成列重为3的H。% 构造一个简单的规则LDPC矩阵 H: 64 x 128, 列重3, 行重6 % 这里采用随机置换的方式仅作为仿真教学不保证纠错性能最优 rng(42); H zeros(64, 128); for col 1:128 % 每列随机选择3个位置置1 pos randperm(64, 3); H(pos, col) 1; end % 确保每行行重为6近似注意随机生成的H可能不是满秩的而且可能存在短环BP译码性能会受影响。实际工程都用代数构造或QC-LDPC。这里只是为了演示BP迭代用这个H跑通流程就够了。3.2 完整的最小和BP译码函数下面给出一个完整的二进制LDPC最小和译码函数输入是接收符号、噪声方差、最大迭代次数输出是译码后的比特。function [decoded_bits, iter_used] bp_decode_min_sum(r, H, sigma2, max_iter) % 最小和BP译码 % 输入: % r: 接收符号向量BPSK映射1-1, 0--1 % H: 稀疏校验矩阵m x n % sigma2: 噪声方差 % max_iter: 最大迭代次数 % 输出: % decoded_bits: 译码后的硬判决比特 % iter_used: 实际迭代次数 [m, n] size(H); % 信道初始LLR: L 2*r / sigma2 L_ch 2 * r(:) / sigma2; % 初始化变量节点消息为信道LLR L_v2c repmat(L_ch, m, 1); % 将H中为0的位置对应的消息置为0表示无连接 L_v2c(~H) 0; % 存储校验节点消息 L_c2v zeros(m, n); for iter 1:max_iter % 1. 校验节点更新 for i 1:m % 找到第i行的非零列索引 nz find(H(i, :)); if isempty(nz) continue; end % 计算该行所有消息的符号乘积和最小绝对值 messages L_v2c(i, nz); abs_m abs(messages); sign_m sign(messages); total_sign prod(sign_m); [min1, idx_min1] min(abs_m); for j nz % 排除自身 mask nz ~ j; if sum(mask) 0 L_c2v(i, j) 0; continue; end if j nz(idx_min1) % 如果是最小值位置找第二小 abs_other abs_m(mask); min2 min(abs_other); else min2 min1; end sign_j total_sign * sign_m(nz j); L_c2v(i, j) sign_j * min2; end end % 2. 变量节点更新 % 后验LLR L_ch 所有来自校验节点的外部消息之和 L_posteriori L_ch sum(L_c2v, 1); % 硬判决 hard_bits L_posteriori 0; % 1对应负LLR % 3. 校验计算伴随式 syndrome mod(hard_bits * H, 2); if all(syndrome 0) decoded_bits hard_bits; iter_used iter; return; end % 生成下一轮变量到校验节点的消息后验LLR减去自身上一次的校验消息 L_v2c repmat(L_posteriori, m, 1) - L_c2v; L_v2c(~H) 0; end decoded_bits L_posteriori 0; iter_used max_iter; end这个函数的关键点有三个消息矩阵L_v2c的大小是m×n但只有H中为1的位置参与计算所以在每次更新后都要用L_v2c(~H)0屏蔽无效位置。硬判决后要计算伴随式hard_bits * H如果为零向量说明译码成功可以提前退出。变量节点更新使用了“后验LLR减去上一次发给同一个校验节点的消息”这个技巧效率高于每次重新累加。3.3 用MATLAB自带LDPC工具对比验证如果你不确定手写译码器是否正确可以和MATLAB Communication Toolbox中的ldpcDecode函数对比。下面是一个简单的验证脚本% 生成随机信息比特 n 64; % 码长 m 32; % 校验位 % 注意需要构造一个有效的生成矩阵这里用ldpcQuasiCyclicMatrix代替 cfg ldpcQuasiCyclicMatrix(1/2, 128); % 码长128 H2 cfg.ParityCheckMatrix; % 编码使用comm.LDPCEncoder encoder comm.LDPCEncoder(H2); decoder comm.LDPCDecoder(H2, DecisionMethod, Hard decision, ... MaximumIterationCount, 10, NumIterationsOutputPort, true);注意comm.LDPCEncoder内部会先对信息比特做校验位计算因此信息位长度不是简单地n-m而是由H的秩决定。使用MATLAB内置对象时必须先通过ldpcQuasiCyclicMatrix构造有效的H而不能用随机生成的满秩不确定的矩阵。对比方法是对同一组接收符号分别调用手写的bp_decode_min_sum和ldpcDecode比较译码输出是否一致。如果稍有差异很可能是最小和近似与内置译码器的归一化因子不同。内置译码器默认使用归一化最小和normalized min-sum归一化因子约为0.75。你可以在自己的实现中加入相同的归一化L_c2v(i, j) alpha * sign_j * min2; % alpha 0.75这段话的意思是归一化最小和与纯最小和的差异在于每条消息乘一个小于1的系数用来补偿最小和近似带来的过估计。加入alpha后性能会贴近标准SPA。4. 多进制LDPC与导航场景MATLAB仿真中的关键设计4.1 多进制LDPC的符号映射与非二进制BP“多进制LDPC”指的是定义在伽罗华域GF(q)上的LDPC码其中q可以是4、8、16、64等。与二进制LDPC每个变量节点对应1个比特不同多进制LDPC的每个变量节点对应一个GF(q)符号。校验矩阵H中的非零元素不再是1而是GF(q)域上的元素例如GF(4)中的0,1,α,α²。BP译码时消息是长度为q的向量表示该符号为各个域元素的概率或LLR。这导致多进制BP的计算复杂度约为二进制q倍。常见简化是使用FFT-QSPA快速傅里叶变换-和积算法把校验节点更新中的循环卷积变成频域乘法。在导航场景中多进制LDPC通常与高阶调制结合比如8PSK、16APSK让一个调制符号携带更多信息比特提高频谱效率。例如北斗B1C信号中使用的BCH码并不算LDPC但新一代测距码提议中经常出现多进制LDPC。这里我们不绑定具体系统只讨论仿真方法。4.2 导航LDPC的典型参数码长、码率、调制方式导航信号通常要求低信噪比、高可靠性因此LDPC码率偏低常见的有1/2、1/3码长从几百到几千。下面是GPS、Galileo、北斗等系统中类似的编码参数对比注意以下数据仅为常见配置参考不代表某个具体标准系统码长码率调制备注GPS L1C12001/2BOC(1,1)外层LDPC内层BCHGalileo E112801/2CBOC采用LDPC内码北斗B-Cnav19441/3BPSK-R5G NR LDPC类似结构如果你的MATLAB工程标题里出现“导航 LDPC”大概率是要在这种参数下做误码率仿真。这时需要把LDPC编码、调制、信道噪声、BP译码串起来。4.3 用MATLAB搭建端到端仿真链路下面给出一个多进制LDPC仿真的最小框架这里以GF(4)上的LDPC为例调制方式使用QPSK相当于每个符号携带2比特。% 参数设置 q 4; % 伽罗华域阶数 nSymbols 64; % 符号数变量节点数 kSymbols nSymbols / 2; % 信息符号数码率1/2 M 4; % QPSK调制阶数 EbNo 3; % 每比特信噪比dB sigma2 1 / (2 * q * 10^(EbNo/10)); % 近似噪声方差 % 生成随机符号0到q-1 info randi([0 q-1], kSymbols, 1); % 实际编码需要GF(q)的生成运算这里用随机校验矩阵代替 % 译码时要求接收符号为实数这里直接演示映射 tx_sym pskmod(info, M, 0, gray); % QPSK调制每个符号对应一个GF(4)元素 rx_sym awgn(tx_sym, EbNo, measured); % BP译码多进制版本这里仅示意 % LLR矩阵大小: nSymbols x q % 计算每个状态的对数概率 llr_table zeros(nSymbols, q); for i 1:nSymbols for g 0:q-1 ref pskmod(g, M, 0, gray); llr_table(i, g1) -abs(rx_sym(i) - ref)^2 / sigma2; end % 归一化转为相对LLR llr_table(i, :) llr_table(i, :) - max(llr_table(i, :)); end这里的要点是多进制LDPC的初始消息是每个符号对所有q个域元素的概率。在MATLAB中用pskmod把GF(q)符号映射到星座点接收点与每个星座点的距离决定初始LLR。注意这里没有实现真正的多进制LDPC校验节点更新因为那需要GF(q)域上的加法和乘法。如果你只是想在导航项目里快速跑通多进制LDPC建议使用现成的有限域工具箱比如MATLAB Communications Toolbox不支持直接的多进制LDPC编码但可以用gf对象手动实现校验运算。一个更实际的做法是把多进制LDPC拆成多个二进制LDPC的比特交织即“二进制LDPC 高阶调制”的BICM方案。这在导航系统中非常常见因为二进制LDPC的BP译码简单成熟多进制映射只影响LLR计算。下面的代码展示了BICM结构% 二进制LDPC编码后的比特交织到QPSK coded_bits randi([0 1], nSymbols*log2(M), 1); % 假设已编码 % 每2比特映射为一个符号 sym_idx bi2de(reshape(coded_bits, log2(M), [])., left-msb); tx_sym pskmod(sym_idx, M, 0, gray); % 接收后计算比特LLR % 对于QPSK每个符号的两个比特LLR可以独立计算 % LLR_bit0 (abs(r - nearest_1)^2 - abs(r - nearest_0)^2) / sigma2这种方式的优点是译码器可以完全复用二进制BP只需要在调制解调器里做软解映射。很多导航载荷实际采用的就是这种“二进制LDPC 高阶调制”的组合因为工程实现简单性能损失也有限。5. 验证与调试用仿真曲线确认你的BP译码器没写错5.1 三个判断译码器正确性的检查点写好的BP译码器决不能直接拿去跑误码率先做三个检查伴随式为零在无噪声或低噪声条件下译码结果必须满足mod(decoded * H, 2) 0。如果无噪声都过不了问题一定在校验节点更新或硬判决逻辑。迭代次数曲线在极低信噪比下最大迭代次数耗尽仍不收敛是正常的但信噪比很高时平均迭代次数应该下降。如果高信噪比下还需要很多次迭代说明消息更新可能被错误放大或缩小。与内置译码器对比用相同参数跑ldpcDecode两者BER曲线差距不能超过0.5dB。超过这个范围优先检查LLR初始化和归一化因子。5.2 常见坑迭代次数、量化精度、校验矩阵行列重坑表现解决方法消息矩阵没有屏蔽无效位置算法占满内存或结果全零每次更新后置零L_v2c(~H)0噪声方差估计不准低信噪比时性能骤降使用awgn的measured选项或先估计SNR最小和alpha因子设为1性能比标准SPA差0.2~0.3dBalpha取0.75或0.8硬判决阈值错误译出码字伴随式不全零检查LLR符号约定明确LLR0表示0还是1H矩阵含有短环高SNR时出现错误平层改用QC-LDPC或增大行列重5.3 一个快速性能测试脚本模板下面给出一个完整测试脚本从生成随机码字到绘制BER曲线你可以直接修改参数使用% 快速性能测试二进制LDPC BPSK 最小和BP clear; clc; n 128; m 64; H makeH(n, m); % 你自己实现的H生成函数 EbNoVec 0:0.5:4; ber zeros(size(EbNoVec)); max_iter 20; for idx 1:length(EbNoVec) EbNo EbNoVec(idx); sigma2 1 / (2 * 10^(EbNo/10)); % BPSK码率1/2已包含 errors 0; totalBits 0; while errors 100 totalBits 1e5 infoBits randi([0 1], n - m, 1); codeword encodeLDPC(H, infoBits); % 生成的码字 tx 1 - 2*codeword; % BPSK映射 rx tx sqrt(sigma2)*randn(size(tx)); [decBits, ~] bp_decode_min_sum(rx, H, sigma2, max_iter); errors errors sum(decBits(:) ~ codeword(:)); totalBits totalBits n; end ber(idx) errors / totalBits; end semilogy(EbNoVec, ber, o-); grid on; xlabel(Eb/N0 (dB)); ylabel(BER);这段脚本里的makeH和encodeLDPC是你自己需要补全的函数。编码时可以用简单的H * u 0求解校验位也可以用高斯消元。如果不想自己写编码直接用comm.LDPCEncoder。注意总比特数要按实际信息位计算上面的脚本为了简单用码长n计算了总比特数严格来说应该用n-m这样曲线会略偏保守但用于验证译码器逻辑已经足够。本文还有配套的精品资源点击获取