简介:面向卫星导航与信道编码研究者的 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的“低密度”体现在哪里
LDPC(Low-Density Parity-Check)码由一个稀疏校验矩阵H定义。所谓“低密度”,是指H中1的数量相对于矩阵规模非常少。考虑一个典型的(64, 32)LDPC码,H是32×64的矩阵,每行可能有6个1,每列也可能只有3个1,1的密度低于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译码有两种等价形式:概率域(SPA,Sum-Product Algorithm)和对数域(LLR域)。概率域里传递的是后验概率,涉及大量乘法,数值容易下溢;实际代码几乎都用对数域。对数域BP中,变量节点向校验节点传递的消息是LLR值,定义为:
L = log(P(x=0)/P(x=1))
信道初始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 L1C | 1200 | 1/2 | BOC(1,1) | 外层LDPC,内层BCH |
| Galileo E1 | 1280 | 1/2 | CBOC | 采用LDPC内码 |
| 北斗B-Cnav | 1944 | 1/3 | BPSK-R | 5G 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, g+1) = -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.3dB | alpha取0.75或0.8 |
| 硬判决阈值错误 | 译出码字伴随式不全零 | 检查LLR符号约定,明确LLR>0表示0还是1 |
| H矩阵含有短环 | 高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,这样曲线会略偏保守,但用于验证译码器逻辑已经足够。
本文还有配套的精品资源,点击获取