简介:该资源为LDPC(低密度奇偶校验码)误比特率仿真MATLAB源码包,面向通信工程、编码理论方向的学生与研究人员,可用于对比软判决与硬判决两类译码策略下的BER性能。包内含8个文件,以.m脚本为主,另附1个xlsx校验矩阵数据文件;脚本覆盖校验矩阵构造、奇偶校验、概率域译码、对数域译码及比特翻转硬判决等关键模块,整体约30KB,结构清晰,适合作为LDPC仿真入门与算法验证的基础工具。截至统计已有592人学习下载。通过运行这些代码,可直观观察不同信噪比条件下的误码率变化,理解软判决相对于硬判决的性能增益,同时为后续改进译码算法或搭建完整通信链路仿真提供可扩展的代码框架。
1. 为什么LDPC编码仿真要把软判决和硬判决分开跑
LDPC码性能仿真有个反直觉的开局:最大的坑不在译码器,而在你有没有把软判决和硬判决当成两套独立链路。标题里的 LPDC 是 LDPC 的常见笔误,LDPC 全称 Low-Density Parity-Check,低密度奇偶校验码,5G NR、Wi-Fi 6、DVB-S2 和光纤通信的前向纠错都依赖它。做物理层链路仿真或接收机芯片验证的人,经常要输出一份 BER 曲线评估译码器性能、迭代次数、量化位宽,而硬判决(比特翻转)和软判决(和积/置信传播)在性能和复杂度上差出好几个 dB,混在一个脚本里往往谁都调不好。这篇文章按一线工程师做链路的样子,从 H 矩阵构造到 BER 仿真代码,把软硬两条路拆开讲清楚,参数怎么设、曲线毛刺怎么排、最后怎么验证仿真没跑偏,全给你落地。
2. 软判决与硬判决LDPC译码的原理差异,以及BER到底差多少
2.1 LDPC编码的H矩阵与Tanner图,仿真前先确认码字合法
LDPC 码的定义方式和其他纠错码不一样。它不直接告诉你怎么把 k 个信息比特映射成 n 个码字比特,而是给一个稀疏奇偶校验矩阵 H,凡是满足H * c^T = 0的向量 c 都是合法码字。H 的稀疏性意味着每个校验方程只约束少数几个比特,这让迭代译码成为可能。
仿真里做 BER 测试时,常见做法是不用完整编码器,直接发全零码字。LDPC 是线性码,BPSK 调制下全零码字和其他码字的错误概率一样,误码统计完全等价。全零码字省掉了生成矩阵 G 的构造,也省掉了高斯消元,H 矩阵只要保证是个合法校验矩阵就行。
H 矩阵可以用 Tanner 图画出来:左边是变量节点(码字比特),右边是校验节点,H 里为 1 的位置就是一条边。仿真前务必检查两点:一是 H 有没有全零行或全零列,有的话说明校验位或变量位没被约束,译码器会出现永不收敛的比特;二是尽量避开长度为 4 的环,两个变量节点被同一对校验节点同时约束时会产生短环,低信噪比时两个环互相强化错误,即使软判决也救不回来。短环检查在构造完 H 之后做一次即可,规则码用 Gallager 构造法天然避免 4 环,随机稀疏矩阵就要自己筛。
2.2 硬判决:比特翻转译码器的工作方式与其边界
硬判决译码的第一步是拿接收符号直接判 0/1,BPSK 下就是 y 大于 0 判 0,小于 0 判 1。这个 0/1 序列送入比特翻转(Bit Flipping,BF)译码器后,每轮迭代做两件事:
- 计算伴随式
s = H * r mod 2,s 中为 1 的位置说明这些校验方程不满足。 - 对每个变量节点统计它连着多少个不满足的校验方程,找出不满足数最多的变量节点,把它翻转,进入下一轮。
当伴随式全为 0 或达到最大迭代次数时停止。BF 算法本质是贪心局部搜索,复杂度低到可以做高速硬件实现,但它有几个明显边界:第一,它只用校验约束的“是/否”信息,完全没用信道给的可信度,所以低信噪比下经常翻转错比特,然后越错越多;第二,两个比特失败数相同时,每次翻转哪个要靠策略,串行翻转还是并行翻转性能会差出一点;第三,陷阱集(trapping set)会让它陷入局部极值,永远收不到全零伴随式。实际写代码时我会加一条保护:连续几轮翻转同一个比特就强制跳出,避免振荡。
2.3 软判决:LLR域的和积译码,以及为什么软判决增益大
软判决译码不丢掉信道的概率信息。BPSK + AWGN 信道下,每个比特的对数似然比定义为:
$$LLR_i = \ln \frac{P(x_i=0|y_i)}{P(x_i=1|y_i)} = \frac{2 y_i}{\sigma^2}$$
y 是接收符号,σ² 是噪声方差。LLR 为正说明该比特大概率是 0,绝对值越大置信度越高。和积译码(Sum-Product Algorithm,SPA)在 Tanner 图上做消息传递,每轮迭代里校验节点收集相邻变量节点的外部信息,再把这些信息汇总后传回;变量节点把信道 LLR 和所有校验节点传来的信息相加,得到后验估计,重新判决。
对数域实现里最常见的写法是用 φ 函数:
$$\varphi(x) = -\ln\left(\tanh\frac{|x|}{2}\right)$$
校验节点的更新变成:外部信息 = 符号乘积 × φ(所有相邻变量信息的 φ 值之和)。这样避免了大量 tanh 乘法,数值稳定性比概率域实现好,也方便后续把消息截断成固定 bit 数做定点仿真。软判决的最大价值在低信噪比区,LLR 天然给置信度加权,一个被噪声拉得很偏的比特不会再像硬判决那样被当成确定值处理,瀑布区(waterfall region)直接往左移几 dB。
2.4 软判决与硬判决的BER对比表与适用场景
| 维度 | 硬判决(BF) | 软判决(LLR SPA) |
|---|---|---|
| 信道信息利用 | 只用符号 | 用完整 LLR 幅度 |
| 消息类型 | 0/1 单比特 | 浮点或量化 LLR |
| 每轮迭代复杂度 | O(E),单比特运算 | O(E),需查表和乘法 |
| 典型收敛迭代次数 | 20~50 | 5~20 |
| 低信噪比 BER | 差 | 明显更好,增益可达 3~6 dB |
| 高信噪比错误平底 | 容易出现 | 相对更干净 |
| 适合场景 | 高速硬件、光模块DSP | 5G NR基站、卫星、深空通信 |
| 迭代次数 | 未收敛时的表现 | 建议值 |
|---|---|---|
| 1~3 | 瀑布区完全出不来 | 太短 |
| 5~10 | 软判决基本收敛 | 浮点软判决可设 8 |
| 20~50 | 硬判决还能拆出一点增益 | 定点仿真按需求加大 |
| 100 以上 | 浮点无意义,定点才需要 | 物理实现兜底 |
3. LDPC的全零码字BER仿真:BPSK+AWGN+两种译码器全代码
3.1 用[n,k]系统形H矩阵生成LDPC校验矩阵
仿真里我不做完整编码器,直接构造一个含单位阵的系统形校验矩阵。设码长 n,校验位数量 m,信息位 k = n - m,H 写成 [A | I_m],左边 A 是 m×k 的随机稀疏矩阵,右边是 m×m 单位阵。这样的好处是任意信息向量 b 都能通过 p = A*b mod 2 得到合法码字,但不发非零码字时也可以只发全零,方便。
构造 A 时控制每一列有 3 个 1,对应规则 LDPC 的列重 3。为了让行重不至于太极端,随机选行位置然后做一次行列去重,去掉全零行和孤立列。完整函数如下:
import numpy as np def build_ldpc_hm(m=128, k=128, col_w=3, seed=7): """构造 [A | I_m] 形式的系统校验矩阵,A 每列恰有 col_w 个 1""" rng = np.random.default_rng(seed) n = m + k H = np.zeros((m, n), dtype=np.uint8) # 左边 A:每列随机选 col_w 个行位置放 1 for col in range(k): rows = rng.choice(m, size=col_w, replace=False) H[rows, col] = 1 # 右边 I_m:保证每行至少一个校验节点能直接约束变量 for r in range(m): H[r, k + r] = 1 # 去掉可能出现的全零行,防止孤立校验节点 nonzero_rows = np.where(H.sum(axis=1) > 0)[0] if len(nonzero_rows) < m: H = H[nonzero_rows, :] return H H = build_ldpc_hm(128, 128) print("H shape:", H.shape, " nnz:", H.sum())代码里我固定了行数和列数各 128,码率正好 1/2。参数里col_w=3对应列重 3,seed控制随机性,换一个 seed 就等于换了一组码。注意我这里没有做 4 环剔除,实际仿真数据量上去后如果曲线在瀑布区出现异常平台,再回来看 H 的围长(girth),用循环遍历变量节点对的最小公共邻居数筛一遍即可。
3.2 在一个Python文件里跑通BER仿真
下面这个完整脚本包含三部分:硬判决比特翻转译码器、软判决和积译码器和 BER 主循环。跑完会输出每个 SNR 点下两种译码器的误比特率。
import numpy as np def bit_flipping_decode(H, recv_bits, max_iter=30): """硬判决:比特翻转译码,返回译码结果和是否收敛""" m, n = H.shape x = recv_bits.copy() last_flip = -1 for it in range(max_iter): syndrome = (H @ x) % 2 if syndrome.sum() == 0: return x, True # 每个变量节点连了多少个不满足的校验 fail_count = H.T @ syndrome # 只翻失败数最大的那个比特,比一次翻多个稳定 target = int(np.argmax(fail_count)) if fail_count[target] == 0: break if target == last_flip: return x, False # 翻转振荡,直接退出 x[target] ^= 1 last_flip = target return x, (H @ x % 2).sum() == 0 def _phi(x): """log-domain 辅助函数 phi(x) = -ln(tanh(|x|/2))""" x = np.clip(np.abs(x), 1e-12, None) return -np.log(np.tanh(x / 2.0)) def ldpc_soft_decode(H, llr_in, max_iter=20): """软判决:LLR 域和积译码,返回硬判决结果""" m, n = H.shape # 初始化变量节点的后验 LLR 为信道 LLR q = llr_in.copy() # 迭代时更新 v2c 和 c2v 消息 for _ in range(max_iter): q_old = q.copy() # 每个校验节点 j 更新相邻变量节点的消息 for j in range(m): col_idx = np.where(H[j, :] == 1)[0] if len(col_idx) < 2: continue msg_in = q_old[col_idx] sign = np.prod(np.sign(msg_in)) abs_sum = np.sum([_phi(x) for x in msg_in]) for idx in col_idx: others = msg_in[msg_in != q_old[idx]] # 外部信息 if len(others) == 0: continue other_sign = np.prod(np.sign(others)) other_phi = np.sum([_phi(x) for x in others]) c2v = other_sign * _phi(other_phi) q[idx] = llr_in[idx] + c2v # 提前判决,收敛就停 hard = (q < 0).astype(int) if ((H @ hard) % 2).sum() == 0: break return (q < 0).astype(int)逻辑说明:bit_flipping_decode里我用H.T @ syndrome统计每个变量节点连接了多少个不满足的校验,np.argmax找到失败数最大的比特翻转。这里故意用“一次只翻一个”而不是“一次翻所有最大值”来避免并行翻转带来的局部振荡,last_flip是防振荡保护。ldpc_soft_decode是标准对数域 SPA,校验节点更新时对每个相邻变量节点计算外部信息,用q_old而不是 q 内部累加,避免同迭代消息互相污染。这个实现为了可读性牺牲了向量化,跑两三个 SNR 点没问题,绩效仿真再改成 numpy 矩阵形式。
3.3 BER主循环与Eb/N0到噪声方差的换算
BER 仿真最容易被忽略的是 Eb/N0 换算。比特能量 Eb 和符号能量 Es 之间差一个码率 r = k/n,BPSK 下 Es = 1,因此噪声方差为:
$$\sigma^2 = \frac{1}{2 \cdot r \cdot 10^{E_b/N_0 (dB)/10}}$$
def ber_simulate(H, ebno_db_list, blocks=2000, max_iter=30): m, n = H.shape k = n - m r = k / n results = {} for ebno in ebno_db_list: sigma2 = 1.0 / (2.0 * r * 10**(ebno / 10.0)) sigma = np.sqrt(sigma2) for mode in ['hard', 'soft']: errs, total_bits = 0, 0 for b in range(blocks): x = np.zeros(n, dtype=np.uint8) # 全零码字 s = 1 - 2 * x.astype(float) # BPSK: 0 -> +1 noise = sigma * np.random.randn(n) y = s + noise # AWGN 信道 llr = 2.0 * y / sigma2 # 信道 LLR if mode == 'hard': dec = bit_flipping_decode( H, (y < 0).astype(np.uint8), max_iter) else: dec = ldpc_soft_decode(H, llr, max_iter=max_iter) errs += int(np.sum(dec != x)) total_bits += n results[(mode, ebno)] = errs / total_bits return results res = ber_simulate(H, [0.5, 1.0, 1.5, 2.0, 2.5], blocks=2000, max_iter=30) for key, val in res.items(): print(key, f"{val:.3e}")参数说明:blocks=2000是每个 SNR 点仿真的码块数,低信噪比下每块错误比特多,2000 块统计足够;高信噪比区误码率到 1e-4 以下时,需要把 blocks 提到 50000 以上。mode='hard'和'soft'两条链路用同一个 H、同一批噪声种子,比较才有意义。循环里每次生成全零码字,所以dec != x里非零的位置就是错比特。
3.4 跑完得到的BER曲线怎么解读
正常跑出来的曲线有两个特征:一是硬判决曲线在低信噪比区比软判决高一个量级以上,差值随 SNR 提高先扩大再缩小;二是软判决曲线有一个明显的瀑布区,在某个 SNR 点附近快速下降。如果看到软判决曲线根本不起跳,先查sigma2换算;如果瀑布区出现平台,先怀疑 H 矩阵列重过低或者行重太不均匀,规则 (3,6) 码的围长保证后再跑一次;如果两条曲线几乎重合,检查llr是不是被整数化成了 0,信道 LLR 计算错误时软判决就等于硬判决。
4. LDPC BER仿真的参数设置、曲线毛刺与仿真license排错
4.1 Eb/N0与Es/N0换算、量化位宽、迭代次数那些最容易设错的参数
| 参数 | 推荐值 | 误设后果 |
|---|---|---|
| 码率 r | 按 k/n 实际计算 | 用 r=1 会让噪声偏小,曲线虚高 |
| 噪声方差 σ² | 1/(2·r·10^(EbN0/10)) | 搞成 Es/N0 则每个点偏左 3dB |
| LLR 截断 | 浮点±20;定点 6~8bit | 截断太小在低SNR直接饱和 |
| BF 迭代上限 | 30 | 太多在陷阱集上空转 |
| SPA 迭代上限 | 8~20 | 3次以下瀑布区出不来 |
| 每SNR块数 | 低SNR 500,高SNR 1e4以上 | 高SNR出现大段毛刺 |
4.2 单块错误不收敛或曲线毛刺:先查校验矩阵,再查随机种子
做定点 LDPC 硬件仿真的人常遇到一个现象:浮点 BP 跑得很好,换成定点 LLR 后 BER 曲线在高信噪比不出错,中间某个点突然冒出几块错误,曲线出现毛刺。常见原因不是量化噪声,而是 H 矩阵中有陷阱集(trapping set)——少数变量节点与校验节点形成的局部结构让迭代消息永远达不到一致的伴随式。浮点时消息幅度大,能冲出去;定点后消息被钳位,就卡在陷阱集里。
处理方式分三步:第一步把 H 的 4 环全筛掉,规则 LDPC 优先用 Gallager 构造法;第二步把量化位数从 4bit 加到 8bit,看毛刺是否消失,如果消失说明是量化问题,不是码结构问题;第三步是对同一 SNR 点换随机种子重跑,毛刺跟着种子跑就说明统计量太低,把 blocks 提高十倍再验证。
4.3 17.1版本verilog仿真license拿不到时,BER曲线还能不能用
有的团队习惯在综合工具里做 RTL 级协同仿真,把 LDPC 译码器写进 FPGA 后跑 BER。常见报错是17.1 error: failure to obtain a verilog simulation license.,这个报错来自仿真工具授权而不是设计本身,跟 LDPC 代码正确性无关。
我一般的做法是把验证拆成两层:算法层用 Python 浮点仿真出 BER 基准曲线,RTL 层只在 2~3 个 SNR 点做定点对照,采样点落在瀑布区和错误平底区。这样 RTL 仿真量小两个量级,即使 license 不可用也完全可以靠算法层曲线完成性能评估。软件代码无论如何都先跑通浮点,license 问题不应该卡住 LDPC BER 仿真本身。
5. 用密度演进校验LDPC仿真结果,抓出“看起来对其实错”的BER曲线
5.1 为什么仿真跑通了还要密度演进
BER 仿真最大的危险不是跑不出来,而是跑出错的曲线还看不出问题。全零码字 + AWGN + 标准 BP 的组合在码长不足或 H 矩阵版本不对时,瀑布区位置可能偏了 1~2 dB,而单看曲线光滑程度看不出异常。密度演进(Density Evolution)是译码领域验证理论阈值的方法:在无限码长假设下,把每个变量节点和校验节点的消息分布当作概率密度函数迭代,跟踪每条消息的 PDF,最后算出该 LDPC 码能收敛的噪声阈值。做仿真验证时只要把实际 BER 曲线跟理论阈值对齐,就能确认识码器和 H 矩阵没选错。
5.2 用高斯近似快速算出阈值,对齐瀑布区位置
规则 (3,6) LDPC 码在 AWGN 信道下可用高斯近似(Gaussian Approximation,GA)快速估阈值:假设变量节点发出的消息服从高斯分布,迭代时只更新均值和方差。这个近似下的经典结果是,BPSK 输入下 (3,6) 码的 BP 阈值大约在 σ*≈0.88 附近,对应 Eb/N0 ≈ 1.1 dB。也就是说,码长足够长的 (3,6) 码,软判决 BER 曲线应该在 1 dB 附近进入瀑布区。如果跑出来的瀑布区在 3 dB 才出现,先查 H 矩阵是不是真的规则 LDPC;如果硬判决曲线在 0 dB 和软判决重合,先查 LLR 初始化。
高斯近似给的是理论极限,有限码长实际仿真会比阈值右移 0.2~0.5 dB,这是正常的,但如果偏离超过 1 dB,一定是 H 矩阵或译码器实现有问题。把这条规则写进仿真脚本的断言里,跑完自动报瀑布区位置,能省掉大量人工看曲线的时间。
5.3 一个可落地的量化扫描脚本套路
最后给一个验证量化的套路:把 LLR 量化位宽作为扫描参数,从 4 bit 到 10 bit 扫一遍,画在同一张 BER 图上,观察 6 bit 和 8 bit 的差距。差距在 0.1 dB 以内就说明量化没有成为瓶颈;差距超过 0.3 dB 则要回到 §4.2 检查陷阱集。这套方法在浮点验证阶段就做完,比 RTL 里灰头土脸改位宽效率高得多。把浮点 BP 作为基准、密度演进阈值作为理论上界,LDPC 的 BER 仿真链路才算闭环。
本文还有配套的精品资源,点击获取