Raptor码与LDPC预编码的MATLAB仿真实现与调优

简介Raptor码是一种无速率喷泉码的增强版本通过LDPC预编码与超图结构实现高鲁棒前向纠错能在数据丢包严重的信道中高效恢复原始消息常用于无线通信、数据存储与网络传输。这份20KB的小巧压缩包内共6个文件3个M脚本分别对应LDPC预编码实现、Raptor编码仿真以及AWGN信道下的性能测试3个MAT文件则保存了不同码长与信噪比条件下的实验数据便于直接绘图分析。已有536人学习下载适合具备MATLAB基础、正在研究喷泉码或信道编码的通信专业学生和工程师。运行仿真程序可观察误码率随信息长度、编码率的变化规律理解接收端随编码符号增多逐步解码的自适应过程修改信道模型或预编码参数还能进一步评估系统在衰落环境中的可靠性为实际工程选型提供参考。1. Raptor码不是一种码而是一套带LDPC预编码的fountain code框架Raptor 码在 RFC 5053、RFC 6330 和卫星组播里出现频率很高但多数人第一次仿真它是被“Raptor”这个名字误导成某种具体纠错码。实际上 Raptor 是喷泉码fountain code的一种结构外层 LT 编码把数据无限往外喷内层 LDPC 预编码把 LT 层没喷全的洞兜住。叠加后接收端只要收到略多于源符号数量的任意一组包就能以逼近 1 的概率恢复全部数据通信仿真里不再需要 ACK 重传。这篇文档把 Raptor 的 LDPC 预编码链路放到 MATLAB 里跑通从两层分工、参数调优、仿真发散排查到结果验证都能直接改 k、码率、开销和随机种子复现。适合刚接触喷泉码仿真的读者也适合把 Raptor 落到网络编码方案的工程师。2. 从fountain code到Raptor码LT喷射层与LDPC预编码层的分工逻辑2.1 纯LT码为什么在中等码长下撑不起常数开销LT 码是最早实用的喷泉码从 k 个源符号里随机取 d 个做异或得到一个输出符号d 按度分布采样。译码端维护一个“剥离”过程——找到一个当前度为 1 的输出符号它的值就是对应的源符号恢复后把这个源符号从其它输出符号里异或掉于是周围输出符号的度集体减 1循环往复。这个过程的终点只有两种全部源符号恢复或者所有剩余输出符号的度都大于 1于是卡死。卡死是纯 LT 的固有毛病。要保证以不低于 1-δ 的概率不卡死鲁棒孤子分布需要把平均度抬到 O(ln k)、把度分布的支撑延伸到 k/R 量级约 √k同时接收端还要额外多收一部分符号来维持剥皮动力。码长 k 从 100 涨到 10000这部分代价在中等码长下非常直观失败率对开销的曲线拖着很长的尾巴而不是一个干净的悬崖。Raptor 的出现就是把“每个输出符号都要能触发剥皮”这个苛刻要求拆掉一半——LT 层不需要恢复全部符号恢复 95% 就收工剩下的交给另一层。2.2 预编码层在Raptor码里的位置先垫一层再往外喷Raptor 的编码流程分两步。预编码把 k 个源符号通过一个码率为 R_p k/k 的 LDPC 码编成 k 个中间符号intermediate symbols前 k 个就是源符号本身后 k-k 个是 LDPC 校验符号。LT 编码度分布改为以 k 为参数对中间符号做异或产生输出符号发往信道。接收端收到 m 个输出符号后先按 LT 剥离恢复中间符号剥不动时把仍为 NaN 的中间符号看成擦除切到 LDPC 擦除译码某一行校验方程如果只剩一个未知变量直接解出来解出新符号后LT 层的新一轮剥离又可以继续。两层交替直到没有新符号被解出。关键在开销口径。LT 层的目标从“恢复 k 个符号”变成“恢复 k 个中间符号”度分布压力小了一大截工程实现里RFC 5053 的分布表平均度可以固定在 10 以下不再随 k 增长而 LDPC 层只管填补残余。总的符号开销是 (1ε_lt)/R_p - 1其中 ε_lt 是 LT 层相对中间符号数的额外接收比例。这里有个常见的仿真口径错误直接拿 m/k 当开销忽略预编码码率惩罚导致不同 R_p 的结果无法横向比较。下表是 k500、ε_lt0.05 时不同预编码码率下的实际总开销。预编码码率 R_p中间符号数 kLDPC 校验行数接收符号数 m相对源符号总开销0.905565658416.8%0.955272755410.8%0.98511115377.4%两级译码的主循环在结构上可以写成下面这段骨架这也是后续仿真里raptor_bp_decode的雏形。% Raptor 两级译码主循环交替处理两类约束直到没有新符号被解出 % recv_syms: 收到的输出符号值邻居表H: 预编码校验矩阵 x nan(kprime, 1); % 中间符号初始全擦除NaN 表示未知 while any(isnan(x)) x_old x; x lt_peel(x, recv_syms); % LT 层度1输出符号触发的剥离 x ldpc_peel(x, H); % 预编码层单未知校验行触发的剥离 if isequaln(x, x_old), break; end % 不动点剥皮彻底卡死 end参数说明kprime是中间符号数recv_syms里只放实际收到的输出符号擦除的符号根本不进这个表H是 m×k 的 LDPC 校验矩阵。isequaln会把 NaN 视为相等所以它能正确比较两轮迭代之间所有未知符号是否原地不动。2.3 为什么预编码层选LDPC稀疏图剥皮天然适配擦除信道预编码不一定要用 LDPC但工程和仿真里几乎都选它原因有三个。第一LDPC 的擦除译码是稀疏图剥皮和 LT 的剥皮是同一个计算范式两层可以直接无缝接力不需要把数据搬到别的域也不需要软信息迭代。第二列重 3 左右的规则 LDPC 在中等码长下就能提供稳定的擦除纠正能力校验行数只有中间符号数的百分之几开销代价可控。第三RS 码和 Turbo 码在这里都不顺手RS 在 GF(2^w) 上做矩阵求逆纠擦除能力固定为 n-k 个符号和喷泉层“收到多少用多少”的接续方式不匹配Turbo 依赖交织器和软判决迭代在 BEC 上要按擦除位置重做度量实现复杂度高一个量级。一句话概括分工LT 层负责“以接近 1 的概率把绝大多数符号剥出来”LDPC 预编码层负责“把剥不出来的残余符号用校验约束补出来”。仿真里验证某一层的好坏要分别统计 LT 层卡住时的残余擦除率以及 LDPC 层补完后是否清零这两层指标在参数调优部分会分开讨论。3. 在MATLAB里搭通Raptor码的LDPC预编码仿真链路最小可复现代码3.1 鲁棒孤子度分布采样c与delta怎么传先写度分布生成和采样这是 LT 层的核心也是最容易写错的一处。function p rsd_pmf(k, c, delta) % k: LT 层处理的符号数Raptor 里是中间符号数 kprime % c, delta: 鲁棒孤子分布参数控制平均度与失败率上界 R c * sqrt(k) * log(k / delta); special floor(k / R); % tau 的尖峰位置 p zeros(1, k); for i 1:k if i 1 rho 1 / k; else rho 1 / (i * (i - 1)); end if i special tau R / (i * k); elseif i special tau R * log(R / delta) / k; else tau 0; end p(i) rho tau; end p p / sum(p); % 归一化必须放在最后一步 endrho是理想孤子分布单独使用时期望度等于 ln k剥皮到后期动力不足tau是在大度位置补出来的“安全垫”让剥皮进行到最后十几轮还有度 1 符号冒出来。c控制 R 的大小delta是允许的译码失败率上界。归一化放在最后是因为 rhotau 的总和大于 1直接拿它当概率会让采样偏向大度。一个隐蔽的坑k 较小时special floor(k/R)可能小于 1整段循环走不到 tau 分支度分布退化成理想孤子分布译码必然卡死。k16 起步做小规模验证时先打印一下 special 是否大于等于 1。采样用逆变换法累积分布算一次即可复用function d sample_rsd(cum, k) % cum: 累积分布调用前算好k 仅用于保证 d 不超过符号数 d find(cum rand(), 1); if isempty(d), d k; end % 浮点边界保护 end3.2 构造可系统编码的LDPC预编码校验矩阵H[A|B]预编码层要满足两个条件校验矩阵能用来做系统编码擦除剥皮时列重不能太低。这里用 H[A|B] 的结构A 是 m×k 的随机部分B 是 m×m 的下双对角矩阵保证可逆。function H make_precoder_H(k, m, colw) % 构造 Raptor 预编码用 LDPC 校验矩阵 H [A | B] % A: m x k 随机列每列 colw 个 1 % B: m x m 下双对角行列式恒为 1保证系统编码可做 A zeros(m, k); for j 1:k A(randperm(m, colw), j) 1; end B eye(m); for j 2:m B(j, j - 1) 1; end H [A, B]; end function x precode_encoder(H, src, k, m) % 系统编码x [src; xp]满足 H*x 0GF(2) A H(:, 1:k); B H(:, k1:end); rhs mod(A * src, 2); % A*src 得校验方程右侧 xp zeros(m, 1); for i 1:m xp(i) rhs(i); if i 1 xp(i) mod(xp(i) xp(i-1), 2); % 回代B 是下双对角 end end x [src; xp]; end系统编码的意思是中间符号的前 k 位直接等于源符号后 m 位是 LDPC 校验位。B 是下双对角回代是 O(m)对 k500、m27 的量级可以忽略。colw取 3列重 2 的 LDPC 在擦除信道下存在大量停止集stopping set剥皮到一半就没有校验行只剩一个未知变量了列重 4 以上没有明显收益且短码长下两列容易高度重合没必要。3.3 主仿真循环擦除信道、两级剥皮译码与失败率统计LT 编码和两级剥皮的完整实现如下。function [y, nbrs] lt_encode(x, num_syms, kprime, c, delta) % x: kprime x 1 中间符号返回 num_syms 个输出符号和邻居表 p rsd_pmf(kprime, c, delta); cum cumsum(p); y zeros(num_syms, 1); nbrs cell(num_syms, 1); for s 1:num_syms d find(cum rand(), 1); idx randperm(kprime, d); % 无放回抽样保证度 d 不被稀释 nbrs{s} idx; y(s) mod(sum(x(idx)), 2); end end function x raptor_bp_decode(y_recv, nbrs_recv, H, kprime) % 两级剥皮LT 输出符号约束 LDPC 校验约束交替执行到不动点 x nan(kprime, 1); % NaN 表示中间符号未恢复 m_recv numel(y_recv); active true(m_recv, 1); % 输出符号是否还剩未知邻居 while true progress false; for s 1:m_recv % LT 层一个输出只剩1个未知邻居 if ~active(s), continue; end nb nbrs_recv{s}; unk nb(isnan(x(nb))); if isempty(unk) active(s) false; elseif numel(unk) 1 known setdiff(nb, unk); x(unk) mod(y_recv(s) sum(x(known)), 2); active(s) false; progress true; end end for i 1:size(H, 1) % LDPC 层一行只剩1个未知变量 vars find(H(i, :)); unk vars(isnan(x(vars))); if numel(unk) 1 x(unk) mod(sum(x(setdiff(vars, unk))), 2); progress true; end end if ~progress || ~any(isnan(x)) break; end end end主仿真循环统计给定开销下的失败率k 500; pre_rate 0.95; m ceil(k / pre_rate) - k; kprime k m; c 0.03; delta 0.5; colw 3; epsilon 0.15; % 相对中间符号数的接收开销 m_recv ceil((1 epsilon) * kprime); trials 200; fails 0; rng(7); H make_precoder_H(k, m, colw); for t 1:trials src randi([0 1], k, 1); x precode_encoder(H, src, k, m); pool 3 * m_recv; % 生成足够大的输出符号池 [y_all, nb_all] lt_encode(x, pool, kprime, c, delta); pick randperm(pool, m_recv); % 等价于擦除信道下收到任意 m_recv 个 x_hat raptor_bp_decode(y_all(pick), nb_all(pick), H, kprime); fails fails any(isnan(x_hat)); end fprintf(k%d, pre_rate%.2f, eps%.2f - 失败率 %.3f\n, ... k, pre_rate, epsilon, fails / trials);pool 3 * m_recv是为了让randperm抽样近似于“从无穷符号流中收到任意 m_recv 个”——输出符号独立同分布这样抽样不会引入位置偏差。剥皮循环直到无新进展才退出复杂度是 O(m_recv·平均度 H 非零元数)k500 量级在单核 MATLAB 里跑 200 次试验只要几十秒。还有一个口径问题接收符号数必须相对 k 而不是 k 去定否则预编码码率越低的方案越“占便宜”把预编码层的工作量错误地归到 LT 层头上。主要参数的作用如下表。参数演示取值作用与调整方向k500源符号数加大则失败率曲线更陡蒙特卡洛更慢pre_rate / m0.95 / 27预编码校验行数调 LDPC 兜底能力colw3LDPC 列重取 2 会引入停止集c / delta0.03 / 0.5LT 度分布安全垫epsilon0.15相对中间符号的接收开销扫描的主变量trials200失败率统计精度目标失败率 1e-2 时建议 ≥5004. Raptor码仿真的参数扫描与仿真发散排查c、预编码码率与剥皮不动点仿真跑通只是开始。工程复现 Raptor 时度分布参数、预编码码率和列重这几个旋钮互相耦合而“仿真发散”在喷泉码仿真里通常不是数值溢出而是剥皮循环提前收敛到一个仍有 NaN 的不动点。这一章把三个旋钮逐一拆开再给排查顺序。4.1 度分布参数c与delta的扫描平均度不是越大越好c 直接决定 R 和尖峰位置 special。前面例子 k527、c0.03、δ0.5 时 R ≈ 0.03×√527×ln(527/0.5) ≈ 4.8special ≈ 110平均度在 67。c 太小0.0050.01时 tau 权重过小度分布接近理想孤子分布剥皮进行到后半程动力枯竭这时 LDPC 层要补的残余超过 m失败率骤升。c 太大0.1 以上时平均度涨到 910每个输出符号的计算量变大但失败率悬崖的位置并不会因此向左移动。δ 的作用更微妙。δ 出现在 R 的对数项里对 R 的影响是对数级的k 从 500 到 5000ln(k/δ) 的变化远小于 √k。仿真里把 δ 设成和目标失败率同量级即可如果蒙特卡洛试验次数只有 200 次却把 δ 设成 1e-4分布上的意义和统计分辨率完全对不上失败率的波动会被误读成参数敏感。固定 k500、pre_rate0.95、ε0.2 时扫描 c趋势如下表。cδ失败率趋势平均度趋势0.010.5偏高剥皮后期卡死约 6.30.030.5明显下降约 6.90.100.5不再改善约 9.50.030.05与 0.5 档接近约 7.1数值会随随机种子浮动但趋势稳定c 存在一个平台区平台左侧失败率快速劣化右侧只有复杂度上升。扫描时保持其它参数和随机种子不动具体做法在第 5 章。4.2 预编码码率与LDPC列重的取舍兜底能力和总开销的交换预编码码率 R_p 决定 m 和总开销的上限列重 colw 决定停止集密度。总开销公式 (1ε_lt)/R_p - 1 说明R_p 从 0.95 降到 0.90多付出约 6 个百分点的总开销多拿到一倍左右的 LDPC 校验行27 行变 56 行。k500、colw3 时各码率的对照如下。预编码码率 R_p校验行数 m校验行平均行重理论可补残余上限ε_lt0.05 总开销0.9056约 305616.8%0.9527约 592710.8%0.9811约 139117.4%行重随 R_p 上升是随机构造的代价校验行越重一个未知变量同时被其它变量约束的概率越大剥皮越依赖“该行只剩 1 个未知”这个巧合。这也是工程上RFC 5053 的 LDPC1/LDPC2 结构用稀疏双对角而不是纯随机矩阵的原因。演示代码保留随机 A是为了直观看到这个行为R_p0.98 时LT 层残余一旦超过 11 个LDPC 层必失败同样残余在 R_p0.90 时完全可补。仿真中建议固定 colw3先扫 R_p 看失败率悬崖移动方向再回头调 c。colw 的另一个坑是列重 2。列重 2 的 LDPC 在 BEC 剥皮下会形成大量停止集剥到最后剩余未恢复的变量和校验行构成环状二部图没有任何行只剩 1 个未知。验证方法很简单同一组参数分别跑 colw2 和 colw3各 200 次试验colw2 的失败率会高出一到两个数量级而且提高 ε 也压不下来——这是停止集在起作用不是参数没调好。4.3 仿真发散排查剥皮提前卡死的四个顺序检查喷泉码仿真里的“不收敛”几乎都是剥皮循环在还剩 NaN 时退出。按下面的顺序排查。第一擦除标记。擦除必须用 NaN不能用 0。真实信道里收到的符号只有 0/10 是合法值把擦除写成 0 会让译码器把大量未知当已知失败率假性下降且恢复内容错误。检查方式是统计x_hat中 NaN 数量的同时核对已恢复符号与源符号是否逐位一致。第二随机流互相污染。randperm、randi、rand共享全局随机流编码器里新增一个抽样调用信道擦除模式就全变了。固定做法是把三个随机流拆开s_deg RandStream(mt19937ar, Seed, 11); % 度分布采样 s_nbr RandStream(mt19937ar, Seed, 22); % 邻居选择 s_chn RandStream(mt19937ar, Seed, 33); % 信道/接收符号选择 % 度采样用 rand(s_deg, 1)邻居用 randperm(s_nbr, kprime, d) % 接收符号选择用 tmp randperm(s_chn, pool); pick tmp(1:m_recv);这样改参数做 A/B 对比时只有你想变的随机性在变。第三邻居抽样必须是无放回。用randi放回抽样时两个邻居可能相同实际度小于名义度等效于度分布整体向小度方向偏移剥皮动力不足。统一用randperm(kprime, d)。第四校验矩阵自身缺陷。检查 H 是否有全零行、全零列以及 B 第一校验列权重退化为 1 的情形。全零行意味着一个校验约束不存在全零列意味着某个中间符号永远不被预编码约束其可靠性被降回纯 LT。提示把预编码层“关掉”做对照是最快的排查动作。raptor_bp_decode里 H 传空矩阵、kprime 传 k走同一条剥皮路径就是纯 LT。如果纯 LT 能恢复而加了预编码反而失败问题一定出在 H 构造或编码一致性上如果两者都失败问题在度分布或擦除标记。最后补一个统计层问题失败率仿真要跑足试验次数。目标失败率 1e-2 时200 次试验平均只能看到 2 次失败置信区间很宽曲线抖动得像“随参数发散”。习惯上每档参数至少跑 500 次或固定跑到 20 次失败再停对比参数时用同一批试验编号。5. 验证Raptor码仿真结果的三个自查技巧消元对照、随机流与曲线口径5.1 小码长下用GF(2)消元对照剥皮译码剥皮算法实现容易出隐蔽错误先别急着跑 500 次蒙特卡洛。取 k16、m4把每个输出符号的异或关系写成 GF(2) 系数矩阵 G_out行对应输出符号列对应中间符号整个链路变成 G_out·x y。对收到的方程做高斯消元解出的 x 应该和raptor_bp_decode的结果完全一致。只要小 k 下两者吻合译码器逻辑基本可信放大 k 时才不会把实现 bug 误判成码的性能问题。function x gf2_gauss(A, y) % A: m x n 系数矩阵(GF(2))y: m x 1未确定的位置返回 NaN [m, n] size(A); M [double(mod(A,2)), double(mod(y,2))]; r 1; for c 1:n p find(M(r:m, c) 1, 1); if isempty(p), continue; end p p r - 1; M([r p], :) M([p r], :); for i 1:m if i ~ r M(i, c) M(i, :) mod(M(i, :) M(r, :), 2); end end r r 1; if r m, break; end end x nan(n, 1); for c 1:n rr find(M(1:m, c) 1); if numel(rr) 1 x(c) M(rr, end); end end end对照时让接收符号数取 3k此时 G_out 几乎必然满秩消元能恢复全部中间符号。另外单独加一句assert(~mod(H * x, 2))验证预编码编码器本身满足校验约束把“编码写错”和“译码写错”这一层先分开。5.2 独立随机流做A/B对比把调参效果和随机波动分开把整条链路包成一个run_raptor_case(k, pre_rate, c, delta, eps, seeds)函数内部用三个独立RandStream分别管度、邻居和信道。跑参数对比时保持另外两组 seed 不变只改目标参数同一参数重复 5 组 seed取失败率的区间而不是单点值。这样“把 c 从 0.03 改成 0.02 后失败率变高”到底是参数影响还是随机波动一眼能判断。5.3 开销-失败率曲线的两个易错点画 pfail 对 ε 的曲线时先检查两个口径。横轴必须换算成相对源符号的总开销 ε_total (1ε_lt)/R_p - 1否则预编码码率不同时曲线不能放在同一张图里对比。纵轴用“源符号恢复出错的比例”而不是“至少一个中间符号未恢复的比例”因为 LT 层和 LDPC 层的失败在预编码码率不同时会互相掩盖。正确做法是最终把src_hat与src逐位比对失败判定一律以源符号为准。曲线形态上Raptor 应该在某个 ε 阈值后呈悬崖式下降然后进入平层。看到线性慢降时先查 c 是否太小再查 H 是否有全零列最后用 5.1 的小码长对照筛掉实现层 bug。平层高度就是该参数组合下的错误底error floor真实系统里由 LDPC 停止集决定。仿真里区分“底来自 LDPC 还是 LT 残余”有个技巧把预编码码率提到 0.98如果平层不降反升说明问题出在 LT 层残余过多如果平层下降说明问题出在预编码层的停止集。把这三步写进仿真脚本当固定检查项之后换信道模型、加码长扫描时才有底气把曲线上的每个拐点解释成编码结构本身的行为而不是脚本或随机性带进来的假象。本文还有配套的精品资源点击获取