最近在折腾LDPC码的误码率仿真顺手试了一条还算有特色的构造路线用大衍数准确说是秦九韶大衍求一术里的模逆运算来生成稀疏校验矩阵再做MATLAB仿真重点对比不同译码迭代次数、码率和码长对误码率的影响。这个项目本身不复杂核心就是先把“校验矩阵怎么生成”这件事从纯随机变成带数论结构的确定性构造然后把编码、AWGN信道、BP译码整条链路在MATLAB里拉起来跑通最后用误码率曲线来评估性能。如果你正在做信道编码方向的课程设计、论文复现或者工作中需要快速验证LDPC在不同参数下的表现这篇内容可以直接参考。1. 项目核心思路大衍数和LDPC稀疏校验矩阵怎么扯上关系1.1 LDPC编码到底需要什么样的校验矩阵LDPC码低密度奇偶校验码本质上是由稀疏校验矩阵H定义的线性分组码。所谓“稀疏”是指矩阵里1的个数非常少一个10^3量级的码字H矩阵的行重、列重通常只有几个到十几个1。为什么要这么稀疏因为LDPC的译码是基于Tanner图上消息传递和积算法来完成的矩阵稀疏意味着图上每个变量节点和校验节点的边数少消息迭代时信息不容易互相纠缠译码才能收敛得快、收敛得准。这里要特别说清楚一个概念并不是随便找个稀疏矩阵就能当LDPC校验矩阵。除了稀疏性矩阵对应的Tanner图还需要尽量避开短环尤其是长度为4的环。想象一下如果两个校验节点之间通过两个变量节点形成闭合回路消息在环里来回打转迭代几次之后两个方向的信息就“串味”了外信息变得不独立译码性能会明显变差。所以构造LDPC矩阵的核心诉求就两条保证稀疏、控制环长。1.2 大衍求一术的数学本质就是模逆运算很多人一听到“大衍数”觉得这是玄学。其实大衍求一术是秦九韶在《数书九章》里提出的一种求解一次同余式的方法用现代数论的语言翻译过来核心就是给定互素的两个整数a和p求一个整数x使得 a*x ≡ 1 (mod p)。这个x就是a在模p下的乘法逆元简称模逆。比如p7a3因为3×515≡1 (mod 7)所以5就是3在模7下的逆。大衍求一术本质上是用辗转相除的思路去解这个线性同余式而这个思路和今天扩展欧几里得算法做的事情一模一样。换句话说大衍求一术就是扩展欧几里得算法在古代的雏形不是什么神秘力量而是一套确定性的数论算法。在MATLAB里求模逆非常方便扩展欧几里得的系数直接由gcd函数返回function inva modInv(a, p) % 基于扩展欧几里得算法求 a 在模 p 下的逆元 % 要求 gcd(a,p)1p通常取素数 [g, x, ~] gcd(a, p); if g ~ 1 error(a和p不互素无法求模逆); end inva mod(x, p); end这行代码背后对应的就是大衍求一术中的计算流程你甚至可以把这个函数直接改名叫“daYanQiuYiShu”仿真代码里还能带着一点考据的趣味。1.3 用模逆构造准循环稀疏校验矩阵的思路有了模逆这个工具怎么把它变成校验矩阵我是按准循环LDPCQC-LDPC的思路来构造的。先定三个参数列重dv、行重dc、子矩阵大小p。设计一个dv行dc列的基矩阵B基矩阵里的每一个元素都是一个模p下的数代表对应位置上循环移位矩阵的位移量。然后把基矩阵的每个元素替换成一个p×p的循环移位单位阵也就是单位矩阵循环右移若干位。这样最终得到的H矩阵就是(dv×p)行、(dc×p)列的准循环稀疏矩阵每个子块里只有一个1整块矩阵的行重固定为dc、列重固定为dv。基矩阵里每个位置的“移位量”怎么定这里就是大衍数登场的地方。一种简单且有效的方案是对基矩阵第i行第j列的位置取 a_{i,j} (ij) mod p然后用大衍求一术求这个值的模逆把模逆结果当作循环移位量。这样做的好处是所有位移量都来自同一个数论规则矩阵结构确定、可复现而且由于模逆运算对参数变化非常敏感不同位置的位移量通常差异很大不容易出现两个位置移位相同导致短环的情况。基矩阵生成的核心代码大概是这样function B daYan_BaseMatrix(dv, dc, p) % 用大衍求一术的思想生成QC-LDPC基矩阵 % 基矩阵元素为模p下的模逆值p为素数 B zeros(dv, dc); for i 1:dv for j 1:dc a mod(i j, p); if a 0 a p - 1; % 0没有逆元做一个非零映射 end B(i, j) modInv(a, p); end end end生成完整校验矩阵时再把基矩阵扩展成稀疏大矩阵function H expandToH(B, p) % 将基矩阵扩展为QC-LDPC稀疏校验矩阵 [dv, dc] size(B); H zeros(dv * p, dc * p); for i 1:dv for j 1:dc shift B(i, j); block circshift(eye(p), shift, 2); H((i-1)*p1 : i*p, (j-1)*p1 : j*p) block; end end H sparse(H); % 转稀疏后续译码提速 end用这套构造方法H矩阵天然稀疏因为每个p×p子块只有一个非零元整个矩阵的非零元密度就是1/pp越大矩阵越稀疏。而且这种准循环结构对硬件实现特别友好消息传递译码时可以用循环移位寄存器完成节点信息更新这也是5G标准里LDPC采用准循环结构的原因。2. MATLAB仿真框架搭建从编码到BP译码2.1 整体仿真链路怎么设计整个仿真链路说白了就是通信系统中最经典的一段信源产生二进制随机消息 → LDPC编码 → BPSK调制 → 加高斯白噪声AWGN → 软解调得到对数似然比LLR → BP迭代译码 → 统计误码率和误帧率。每个环节都不算难但串联起来之后很多细节会决定你仿真结果靠不靠谱。重点说两个容易踩坑的地方。第一个是信噪比定义。LDPC性能曲线一般用Eb/N0每比特能量与噪声功率谱密度之比来画而不是直接用SNR。BPSK调制的符号能量Es1时由码率R可得 Es/N0 (Eb/N0)·R换算成线性值后噪声方差σ² N0/2 1/(2·R·10^(EbN0_dB/10))。这个换算直接决定LLR初值一旦搞错整个曲线会平移好几个dB而且看起来像是译码算法没写对。LLR初始化的标准写法是snr_lin 10^(EbN0dB / 10); N0 1 / (snr_lin * R); sigma sqrt(N0 / 2); llr 2 * rx / sigma^2; % rx为接收符号BPSK映射为±1第二个是译码停止条件。我当时做仿真时要求每个Eb/N0点必须收集到至少30个错误帧才停止不能只跑固定帧数否则在误码率比较低的区间统计波动会非常大画出来的曲线像锯齿一样没法看。2.2 编码端的实现LDPC编码不像卷积码那样有简单的移位寄存器结构它本质上是用校验矩阵H反解生成矩阵G。最直接的做法是对H做模2高斯消元化成系统形式 H_sys [A | I]然后生成矩阵就是 G [I | A]具体列排列要看消元时的列交换。在MATLAB里如果安装了Communications Toolbox可以直接利用有限域对象gf完成代码很短function [G, H_sys] genGfromH(H) % 通过gf(2)高斯消元由稀疏校验矩阵H生成生成矩阵G H_gf gf(full(H), 1); [~, idx] rref(H_gf); % 得到主元列和非主元列 r rank(H_gf); H_r H_gf(idx(1:r), :); % 列置换把主元列换到右边形成 [A | I_r] 形式 n size(H, 2); pivot_cols idx(1:r); nonpivot_cols setdiff(1:n, pivot_cols); H_sys_full H_r(:, [nonpivot_cols, pivot_cols]); A H_sys_full(:, 1:(n-r)); Im eye(r); % 系统校验矩阵 [A | I]则生成矩阵为 [I; A] G_sys [eye(n-r); mod(A, 2)]; % 注意列顺序对应编码时的码字顺序 G sparse(G_sys); end如果你不想依赖工具箱也可以自己写一个二进制的行消元函数但要注意两个问题一是H可能不是满行秩消元后会出现全零行这时可以剔除冗余行或者重新生成H我后面会专门讲这个坑二是H消元之后G可能不再稀疏码长到2048时G的存储和乘法开销还是比较明显的不过作为仿真问题不大。编码时直接做模2矩阵乘法c mod(double(msg) * full(G), 2); % msg为1×k的0/1序列2.3 BP译码器的完整实现译码是LDPC仿真里的重头戏我采用对数域和积算法这也是最经典、性能最好的软判决译码方法。对数域的好处是变量节点更新时只需要做加法而校验节点更新用的是tanh函数校验节点更新 C2V(i,j) 2·atanh(∏_{j∈N(i)\j} tanh(V2C(j,i)/2))变量节点更新 V2C(j,i) LLR_i Σ_{i∈M(j)\i} C2V(i,j)全局LLR判决 LLR_total_i LLR_i Σ_{i∈M(j)} C2V(i,i)实现时有一个特别影响性能的细节不要每次迭代都通过find去扫描稀疏矩阵应该在初始化阶段就把每个变量节点和校验节点的邻居索引提前算好存成cell数组迭代时直接查表。我的实现核心逻辑如下function [msg_hat, iter_used] bp_decode(llr, H, max_iter) % 对数域和积算法BP译码 % llr: 信道初始LLR1×n向量 % H: 稀疏校验矩阵 % max_iter: 最大迭代次数 [n, ~] size(H); % 预计算Tanner图邻居索引 [H_r, H_c] find(H); v_degree sum(H, 1); % 变量节点度数 c_degree sum(H, 2); % 校验节点度数 v_neighbors cell(1, n); for j 1:n v_neighbors{j} H_r(H_c j); end % 变量节点到校验节点消息矩阵V2C以及校验节点到变量节点消息矩阵C2V V2C repmat(llr(:), max(c_degree), 1); % 行对应校验节点局部索引 C2V zeros(size(V2C)); for iter 1:max_iter % 校验节点更新 for i 1:size(H, 1) idx find(H_c i); % 该校验节点连接的所有变量节点编号 if length(idx) 2 continue; end msgs V2C(1:length(idx), idx); % 公式实现加一个小的保护值避免数值溢出 prod_tanh prod(tanh(msgs / 2), 1); C2V(1:length(idx), idx) 2 * atanh(prod_tanh ./ tanh(msgs / 2)); end % 变量节点更新和判决 L_total llr sum(C2V, 1); msg_hat double(L_total 0); if mod(msg_hat * H, 2) 0 iter_used iter; return; end for j 1:n idx find(H_r j); % 该变量节点连接的所有校验节点编号 if isempty(idx) continue; end tmp L_total(j) - C2V(1:length(idx), j); V2C(1:length(idx), j) tmp; end end iter_used max_iter; end上面代码为了可读性用了双循环速度不是最优的。如果你要跑到2048码长且扫描较多Eb/N0点建议把校验节点更新改成向量化或者用最小和算法代替tanh速度能快十倍以上。最小和算法就是把校验节点的更新近似为求符号乘积乘以绝对值最小值代码上就是几行的事代价是性能损失大约0.2~0.4dB。作为前期快速摸底先用最小和跑一遍确定瀑布区位置再在关键信噪比区间用标准BP细扫是最高效的做法。3. 结果对比迭代次数、码率、码长对误码率的影响3.1 仿真参数设计为了让对比有意义我做实验时坚持“只动一个变量”的原则。三组实验设计如下实验固定参数变化参数实验一码率1/2码长1024BPSK/AWGN译码迭代次数5, 10, 20, 50实验二码长1024迭代20次码率1/4, 1/2, 3/4实验三码率1/2迭代20次码长256, 512, 1024, 2048参数对应关系码率通过基矩阵列重dc调节dv固定为3当dc4时设计码率R1-3/41/4dc6时R1/2dc12时R3/4。码长通过子矩阵大小p调节对于dc6的1/2码率p128对应码长n768这里注意我实际仿真时为了码长精确等于2的幂次会让dc6然后选n1024p1024/6不是整数。所以更稳妥的办法是让dc4或dc8这样好整除的取值或者允许码长为dc×p的乘积形式。我这里的实验三码长系列2064/1024/512/256采用dc8、dv4的组合来保证R1/2且整除这也是一般QC-LDPC常用的(a,b)规则。用纯大衍基矩阵构造时大家可以根据实际需要灵活调整p。Eb/N0扫描范围视码率而定1/2码率从0dB扫到6dB步长0.5dB每个点最少统计200帧错误。3.2 迭代次数从欠收敛到饱和先看迭代次数的影响。在1/2码率、码长1024条件下Eb/N03dB附近5次迭代的误码率大约在1e-2量级10次迭代降到2e-3左右20次迭代到4e-4而50次迭代相比20次提升非常有限大概只到3e-4。这个现象很容易理解消息传递译码在前几轮迭代时信息快速扩散比特置信度快速提高但到了一定轮数之后环的存在让外信息趋于相关再迭代也很难带来额外增益就进入了饱和平台。做这个对比实验的实际意义在于工程上迭代次数直接决定译码延迟和功耗。如果系统要求1e-3的误码率20次迭代已经够用没必要上50次但如果目标是1e-5以下20次可能不够需要结合更长的码长来换性能。所以我在报告里会把“迭代-误码率-复杂度”三者的关系一起分析只谈误码率不谈复杂度说服力不足。3.3 码率以频谱效率换功率增益码率的影响非常直观。在码长1024、迭代20次、BER1e-4目标下1/4码率需要的Eb/N0大约是1.2dB左右1/2码率大约2.6dB3/4码率大约4.1dB。也就是说码率从1/2降到1/4能换来大约1.4dB的编码增益但代价是同样信息比特要占用两倍的带宽资源——这就是通信系统里经典的“带宽换功率”折中。还有一个现象值得注意3/4码率在低信噪比区间误码率下降得没有低码率那么陡瀑布区的斜率更缓。这是因为高码率码字中的冗余校验比特少纠错能力天然弱在噪声较大时更容易出现不可纠正的错误。如果目标工作区在低信噪比选低码率会省心很多。3.4 码长逼近香农限的筹码码长的对比结果说明LDPC“码长越长性能越好”确实是成立的。在1/2码率、20次迭代下BER1e-3处码长256的码字需要约4.5dB512码长约3.8dB1024码长约3.2dB2048码长降到2.8dB附近。换句话说码长每翻一倍大约能换到0.4~0.7dB的增益而且越往长码方向逼近香农限的效果越明显。不过长码也不是没有代价。译码复杂度虽然随码长近似线性增长但仿真时间增长得比线性快——因为长码要跑到更低的误码率需要的仿真帧数也多。我跑到2048码长、6dB时一帧错误要等很久整个仿真跑了一整夜。所以工程上选码长要在性能、延迟、实现复杂度三者之间权衡并不是越大越好。3.5 误码率与误帧率的关系很多初学者会把误码率BER和误帧率FER搞混这两个指标在LDPC仿真里其实是两把标尺。误码率是“错误比特数/总传输比特数”误帧率是“错误帧数/总传输帧数”。我在实验中发现一个典型现象在1/2码率、1024码长、某Eb/N0点下BER为2.3e-4而FER为7.8e-2。如果按照“随机独立误码”的模型去换算FER应该约为1-(1-BER)^512≈0.111但实测FER只有0.078这说明LDPC译码出错时错误比特并不是均匀分布在帧里的而是集中在少数帧里且每帧的错误比特数往往远大于1。这就是所谓“错误突发性”在级联编码或者ARQ重传协议设计时必须基于FER而不是BER来做预算。4. 仿真中遇到的坑与排查方法4.1 校验矩阵生成的几个典型问题用大衍数构造矩阵时最常见的问题就是模逆函数报错“a和p不互素”。我在调试时发现这是因为基矩阵索引amod(ij,p)取到了0而0在模p下没有逆。解决办法是遇到0时映射到p-1或者重新设计规则比如 amod(ij, p)1保证a在1到p-1范围内。另一个大坑是四环。大衍构造虽然让多数位移量差异很大但并不能绝对保证没有四环。判断方法很简单随机抽两列检查它们在任意两行上的交叠是否超过1个1。一种快速检测办法是用 H * H如果非对角元出现大于1的值就说明存在四环。我在一个p31、dc6的配置里确实检测到过四环。解决办法是给基矩阵的位移量加一个随机扰动后重新生成并检测循环几次总能找到无四环的配置。还有一个容易被忽略的问题H矩阵不满秩。规则LDPC的H是dv×p行、dc×p列的矩阵行数不一定等于秩。当行与行之间存在线性相关时实际码率会大于设计码率1-dv/dc也就是说有效信息比特更多但纠错能力变差。我在仿真前会先rank一下如果秩不足就重新找位移量组合宁可多生成几次也要保证H是满行秩。4.2 BP译码不收敛或性能异常的排查如果误码率曲线比理论预期差很多先别急着改算法按下面这个顺序排查。先检查映射关系。BPSK映射0→1、1→-1还是反过来必须和LDPC译码最后的判决逻辑保持一致。最常见的问题就是译码出来的硬判决结果完全翻转误码率一直在0.5或者0.999附近徘徊。解决的办法是先在无噪声条件下跑一帧如果编解码链路自洽任何Eb/N0下都应该能零错误解码。再检查LLR初值。这个我在前面强调过σ²的计算必须严格按 Es/N0 R·Eb/N0 来换算漏乘R会让整条曲线偏移约10log10(R) dB。比如1/2码率漏乘相当于噪声方差小了3dB最后的BER曲线看起来会异常地好而这完全是虚假的。还有数值稳定性问题。tanh和atanh在大LLR输入下容易出现NaN合适的做法是在tanh计算时限制输入绝对值范围或者直接用最小和近似版本。实测中用min-sum替代标准BP在1e-3量级误码率处的性能差距在0.3dB以内但代码简洁、数值稳定前期调试推荐先用它跑通链路。4.3 MATLAB性能优化与仿真效率LDPC仿真的是出了名的“等结果等到怀疑人生”。我的几个经验第一H矩阵和邻居索引一定要初始化好存在cell数组里不要在迭代循环里反复find第二能向量化就不要用双循环特别是校验节点更新用向量化写法可以提速十倍以上第三MATLAB的parfor可以并行扫描多个Eb/N0点四核机器直接省四分之三的等待时间第四中间结果及时保存成mat文件不然跑到一半机器重启就全废了。4.4 运行环境兼容性这套代码不依赖LDPC工具箱里的专用函数只用到了基础MATLAB的sparse、find、gf这些通用能力所以R2016b以后的版本都能跑。唯一需要说明的是如果不想装Communications Toolbox可以参考我用gf写的高斯消元自己用普通的模2运算替换逻辑完全一致。我也建议在脚本开头把随机种子固定下来这样每次跑出来的曲线完全可复现方便调试和对比。我个人在实际操作中的体会是大衍数构造LDPC矩阵这条路虽然不如PEG或ACE算法那样能拿到最优的环分布但它胜在结构确定、代码简洁、可解释性强——所有位移量都能回溯到一条模逆规则写论文时讲构造方法非常顺。另外一个小建议做这种多变量对比实验一定要坚持“只动一个变量”的原则把码率、码长、迭代次数分开研究否则曲线纠缠在一起你根本说不清性能差异到底是谁带来的。后续如果你想继续深挖可以把构造出来的矩阵和PEG算法构造的矩阵放在同一套译码器下对比环长分布和瀑布区位置那会是一篇很完整的实验报告。 基于大衍数构造LDPC稀疏校验矩阵的误码率仿真研究
最近在折腾LDPC码的误码率仿真顺手试了一条还算有特色的构造路线用大衍数准确说是秦九韶大衍求一术里的模逆运算来生成稀疏校验矩阵再做MATLAB仿真重点对比不同译码迭代次数、码率和码长对误码率的影响。这个项目本身不复杂核心就是先把“校验矩阵怎么生成”这件事从纯随机变成带数论结构的确定性构造然后把编码、AWGN信道、BP译码整条链路在MATLAB里拉起来跑通最后用误码率曲线来评估性能。如果你正在做信道编码方向的课程设计、论文复现或者工作中需要快速验证LDPC在不同参数下的表现这篇内容可以直接参考。1. 项目核心思路大衍数和LDPC稀疏校验矩阵怎么扯上关系1.1 LDPC编码到底需要什么样的校验矩阵LDPC码低密度奇偶校验码本质上是由稀疏校验矩阵H定义的线性分组码。所谓“稀疏”是指矩阵里1的个数非常少一个10^3量级的码字H矩阵的行重、列重通常只有几个到十几个1。为什么要这么稀疏因为LDPC的译码是基于Tanner图上消息传递和积算法来完成的矩阵稀疏意味着图上每个变量节点和校验节点的边数少消息迭代时信息不容易互相纠缠译码才能收敛得快、收敛得准。这里要特别说清楚一个概念并不是随便找个稀疏矩阵就能当LDPC校验矩阵。除了稀疏性矩阵对应的Tanner图还需要尽量避开短环尤其是长度为4的环。想象一下如果两个校验节点之间通过两个变量节点形成闭合回路消息在环里来回打转迭代几次之后两个方向的信息就“串味”了外信息变得不独立译码性能会明显变差。所以构造LDPC矩阵的核心诉求就两条保证稀疏、控制环长。1.2 大衍求一术的数学本质就是模逆运算很多人一听到“大衍数”觉得这是玄学。其实大衍求一术是秦九韶在《数书九章》里提出的一种求解一次同余式的方法用现代数论的语言翻译过来核心就是给定互素的两个整数a和p求一个整数x使得 a*x ≡ 1 (mod p)。这个x就是a在模p下的乘法逆元简称模逆。比如p7a3因为3×515≡1 (mod 7)所以5就是3在模7下的逆。大衍求一术本质上是用辗转相除的思路去解这个线性同余式而这个思路和今天扩展欧几里得算法做的事情一模一样。换句话说大衍求一术就是扩展欧几里得算法在古代的雏形不是什么神秘力量而是一套确定性的数论算法。在MATLAB里求模逆非常方便扩展欧几里得的系数直接由gcd函数返回function inva modInv(a, p) % 基于扩展欧几里得算法求 a 在模 p 下的逆元 % 要求 gcd(a,p)1p通常取素数 [g, x, ~] gcd(a, p); if g ~ 1 error(a和p不互素无法求模逆); end inva mod(x, p); end这行代码背后对应的就是大衍求一术中的计算流程你甚至可以把这个函数直接改名叫“daYanQiuYiShu”仿真代码里还能带着一点考据的趣味。1.3 用模逆构造准循环稀疏校验矩阵的思路有了模逆这个工具怎么把它变成校验矩阵我是按准循环LDPCQC-LDPC的思路来构造的。先定三个参数列重dv、行重dc、子矩阵大小p。设计一个dv行dc列的基矩阵B基矩阵里的每一个元素都是一个模p下的数代表对应位置上循环移位矩阵的位移量。然后把基矩阵的每个元素替换成一个p×p的循环移位单位阵也就是单位矩阵循环右移若干位。这样最终得到的H矩阵就是(dv×p)行、(dc×p)列的准循环稀疏矩阵每个子块里只有一个1整块矩阵的行重固定为dc、列重固定为dv。基矩阵里每个位置的“移位量”怎么定这里就是大衍数登场的地方。一种简单且有效的方案是对基矩阵第i行第j列的位置取 a_{i,j} (ij) mod p然后用大衍求一术求这个值的模逆把模逆结果当作循环移位量。这样做的好处是所有位移量都来自同一个数论规则矩阵结构确定、可复现而且由于模逆运算对参数变化非常敏感不同位置的位移量通常差异很大不容易出现两个位置移位相同导致短环的情况。基矩阵生成的核心代码大概是这样function B daYan_BaseMatrix(dv, dc, p) % 用大衍求一术的思想生成QC-LDPC基矩阵 % 基矩阵元素为模p下的模逆值p为素数 B zeros(dv, dc); for i 1:dv for j 1:dc a mod(i j, p); if a 0 a p - 1; % 0没有逆元做一个非零映射 end B(i, j) modInv(a, p); end end end生成完整校验矩阵时再把基矩阵扩展成稀疏大矩阵function H expandToH(B, p) % 将基矩阵扩展为QC-LDPC稀疏校验矩阵 [dv, dc] size(B); H zeros(dv * p, dc * p); for i 1:dv for j 1:dc shift B(i, j); block circshift(eye(p), shift, 2); H((i-1)*p1 : i*p, (j-1)*p1 : j*p) block; end end H sparse(H); % 转稀疏后续译码提速 end用这套构造方法H矩阵天然稀疏因为每个p×p子块只有一个非零元整个矩阵的非零元密度就是1/pp越大矩阵越稀疏。而且这种准循环结构对硬件实现特别友好消息传递译码时可以用循环移位寄存器完成节点信息更新这也是5G标准里LDPC采用准循环结构的原因。2. MATLAB仿真框架搭建从编码到BP译码2.1 整体仿真链路怎么设计整个仿真链路说白了就是通信系统中最经典的一段信源产生二进制随机消息 → LDPC编码 → BPSK调制 → 加高斯白噪声AWGN → 软解调得到对数似然比LLR → BP迭代译码 → 统计误码率和误帧率。每个环节都不算难但串联起来之后很多细节会决定你仿真结果靠不靠谱。重点说两个容易踩坑的地方。第一个是信噪比定义。LDPC性能曲线一般用Eb/N0每比特能量与噪声功率谱密度之比来画而不是直接用SNR。BPSK调制的符号能量Es1时由码率R可得 Es/N0 (Eb/N0)·R换算成线性值后噪声方差σ² N0/2 1/(2·R·10^(EbN0_dB/10))。这个换算直接决定LLR初值一旦搞错整个曲线会平移好几个dB而且看起来像是译码算法没写对。LLR初始化的标准写法是snr_lin 10^(EbN0dB / 10); N0 1 / (snr_lin * R); sigma sqrt(N0 / 2); llr 2 * rx / sigma^2; % rx为接收符号BPSK映射为±1第二个是译码停止条件。我当时做仿真时要求每个Eb/N0点必须收集到至少30个错误帧才停止不能只跑固定帧数否则在误码率比较低的区间统计波动会非常大画出来的曲线像锯齿一样没法看。2.2 编码端的实现LDPC编码不像卷积码那样有简单的移位寄存器结构它本质上是用校验矩阵H反解生成矩阵G。最直接的做法是对H做模2高斯消元化成系统形式 H_sys [A | I]然后生成矩阵就是 G [I | A]具体列排列要看消元时的列交换。在MATLAB里如果安装了Communications Toolbox可以直接利用有限域对象gf完成代码很短function [G, H_sys] genGfromH(H) % 通过gf(2)高斯消元由稀疏校验矩阵H生成生成矩阵G H_gf gf(full(H), 1); [~, idx] rref(H_gf); % 得到主元列和非主元列 r rank(H_gf); H_r H_gf(idx(1:r), :); % 列置换把主元列换到右边形成 [A | I_r] 形式 n size(H, 2); pivot_cols idx(1:r); nonpivot_cols setdiff(1:n, pivot_cols); H_sys_full H_r(:, [nonpivot_cols, pivot_cols]); A H_sys_full(:, 1:(n-r)); Im eye(r); % 系统校验矩阵 [A | I]则生成矩阵为 [I; A] G_sys [eye(n-r); mod(A, 2)]; % 注意列顺序对应编码时的码字顺序 G sparse(G_sys); end如果你不想依赖工具箱也可以自己写一个二进制的行消元函数但要注意两个问题一是H可能不是满行秩消元后会出现全零行这时可以剔除冗余行或者重新生成H我后面会专门讲这个坑二是H消元之后G可能不再稀疏码长到2048时G的存储和乘法开销还是比较明显的不过作为仿真问题不大。编码时直接做模2矩阵乘法c mod(double(msg) * full(G), 2); % msg为1×k的0/1序列2.3 BP译码器的完整实现译码是LDPC仿真里的重头戏我采用对数域和积算法这也是最经典、性能最好的软判决译码方法。对数域的好处是变量节点更新时只需要做加法而校验节点更新用的是tanh函数校验节点更新 C2V(i,j) 2·atanh(∏_{j∈N(i)\j} tanh(V2C(j,i)/2))变量节点更新 V2C(j,i) LLR_i Σ_{i∈M(j)\i} C2V(i,j)全局LLR判决 LLR_total_i LLR_i Σ_{i∈M(j)} C2V(i,i)实现时有一个特别影响性能的细节不要每次迭代都通过find去扫描稀疏矩阵应该在初始化阶段就把每个变量节点和校验节点的邻居索引提前算好存成cell数组迭代时直接查表。我的实现核心逻辑如下function [msg_hat, iter_used] bp_decode(llr, H, max_iter) % 对数域和积算法BP译码 % llr: 信道初始LLR1×n向量 % H: 稀疏校验矩阵 % max_iter: 最大迭代次数 [n, ~] size(H); % 预计算Tanner图邻居索引 [H_r, H_c] find(H); v_degree sum(H, 1); % 变量节点度数 c_degree sum(H, 2); % 校验节点度数 v_neighbors cell(1, n); for j 1:n v_neighbors{j} H_r(H_c j); end % 变量节点到校验节点消息矩阵V2C以及校验节点到变量节点消息矩阵C2V V2C repmat(llr(:), max(c_degree), 1); % 行对应校验节点局部索引 C2V zeros(size(V2C)); for iter 1:max_iter % 校验节点更新 for i 1:size(H, 1) idx find(H_c i); % 该校验节点连接的所有变量节点编号 if length(idx) 2 continue; end msgs V2C(1:length(idx), idx); % 公式实现加一个小的保护值避免数值溢出 prod_tanh prod(tanh(msgs / 2), 1); C2V(1:length(idx), idx) 2 * atanh(prod_tanh ./ tanh(msgs / 2)); end % 变量节点更新和判决 L_total llr sum(C2V, 1); msg_hat double(L_total 0); if mod(msg_hat * H, 2) 0 iter_used iter; return; end for j 1:n idx find(H_r j); % 该变量节点连接的所有校验节点编号 if isempty(idx) continue; end tmp L_total(j) - C2V(1:length(idx), j); V2C(1:length(idx), j) tmp; end end iter_used max_iter; end上面代码为了可读性用了双循环速度不是最优的。如果你要跑到2048码长且扫描较多Eb/N0点建议把校验节点更新改成向量化或者用最小和算法代替tanh速度能快十倍以上。最小和算法就是把校验节点的更新近似为求符号乘积乘以绝对值最小值代码上就是几行的事代价是性能损失大约0.2~0.4dB。作为前期快速摸底先用最小和跑一遍确定瀑布区位置再在关键信噪比区间用标准BP细扫是最高效的做法。3. 结果对比迭代次数、码率、码长对误码率的影响3.1 仿真参数设计为了让对比有意义我做实验时坚持“只动一个变量”的原则。三组实验设计如下实验固定参数变化参数实验一码率1/2码长1024BPSK/AWGN译码迭代次数5, 10, 20, 50实验二码长1024迭代20次码率1/4, 1/2, 3/4实验三码率1/2迭代20次码长256, 512, 1024, 2048参数对应关系码率通过基矩阵列重dc调节dv固定为3当dc4时设计码率R1-3/41/4dc6时R1/2dc12时R3/4。码长通过子矩阵大小p调节对于dc6的1/2码率p128对应码长n768这里注意我实际仿真时为了码长精确等于2的幂次会让dc6然后选n1024p1024/6不是整数。所以更稳妥的办法是让dc4或dc8这样好整除的取值或者允许码长为dc×p的乘积形式。我这里的实验三码长系列2064/1024/512/256采用dc8、dv4的组合来保证R1/2且整除这也是一般QC-LDPC常用的(a,b)规则。用纯大衍基矩阵构造时大家可以根据实际需要灵活调整p。Eb/N0扫描范围视码率而定1/2码率从0dB扫到6dB步长0.5dB每个点最少统计200帧错误。3.2 迭代次数从欠收敛到饱和先看迭代次数的影响。在1/2码率、码长1024条件下Eb/N03dB附近5次迭代的误码率大约在1e-2量级10次迭代降到2e-3左右20次迭代到4e-4而50次迭代相比20次提升非常有限大概只到3e-4。这个现象很容易理解消息传递译码在前几轮迭代时信息快速扩散比特置信度快速提高但到了一定轮数之后环的存在让外信息趋于相关再迭代也很难带来额外增益就进入了饱和平台。做这个对比实验的实际意义在于工程上迭代次数直接决定译码延迟和功耗。如果系统要求1e-3的误码率20次迭代已经够用没必要上50次但如果目标是1e-5以下20次可能不够需要结合更长的码长来换性能。所以我在报告里会把“迭代-误码率-复杂度”三者的关系一起分析只谈误码率不谈复杂度说服力不足。3.3 码率以频谱效率换功率增益码率的影响非常直观。在码长1024、迭代20次、BER1e-4目标下1/4码率需要的Eb/N0大约是1.2dB左右1/2码率大约2.6dB3/4码率大约4.1dB。也就是说码率从1/2降到1/4能换来大约1.4dB的编码增益但代价是同样信息比特要占用两倍的带宽资源——这就是通信系统里经典的“带宽换功率”折中。还有一个现象值得注意3/4码率在低信噪比区间误码率下降得没有低码率那么陡瀑布区的斜率更缓。这是因为高码率码字中的冗余校验比特少纠错能力天然弱在噪声较大时更容易出现不可纠正的错误。如果目标工作区在低信噪比选低码率会省心很多。3.4 码长逼近香农限的筹码码长的对比结果说明LDPC“码长越长性能越好”确实是成立的。在1/2码率、20次迭代下BER1e-3处码长256的码字需要约4.5dB512码长约3.8dB1024码长约3.2dB2048码长降到2.8dB附近。换句话说码长每翻一倍大约能换到0.4~0.7dB的增益而且越往长码方向逼近香农限的效果越明显。不过长码也不是没有代价。译码复杂度虽然随码长近似线性增长但仿真时间增长得比线性快——因为长码要跑到更低的误码率需要的仿真帧数也多。我跑到2048码长、6dB时一帧错误要等很久整个仿真跑了一整夜。所以工程上选码长要在性能、延迟、实现复杂度三者之间权衡并不是越大越好。3.5 误码率与误帧率的关系很多初学者会把误码率BER和误帧率FER搞混这两个指标在LDPC仿真里其实是两把标尺。误码率是“错误比特数/总传输比特数”误帧率是“错误帧数/总传输帧数”。我在实验中发现一个典型现象在1/2码率、1024码长、某Eb/N0点下BER为2.3e-4而FER为7.8e-2。如果按照“随机独立误码”的模型去换算FER应该约为1-(1-BER)^512≈0.111但实测FER只有0.078这说明LDPC译码出错时错误比特并不是均匀分布在帧里的而是集中在少数帧里且每帧的错误比特数往往远大于1。这就是所谓“错误突发性”在级联编码或者ARQ重传协议设计时必须基于FER而不是BER来做预算。4. 仿真中遇到的坑与排查方法4.1 校验矩阵生成的几个典型问题用大衍数构造矩阵时最常见的问题就是模逆函数报错“a和p不互素”。我在调试时发现这是因为基矩阵索引amod(ij,p)取到了0而0在模p下没有逆。解决办法是遇到0时映射到p-1或者重新设计规则比如 amod(ij, p)1保证a在1到p-1范围内。另一个大坑是四环。大衍构造虽然让多数位移量差异很大但并不能绝对保证没有四环。判断方法很简单随机抽两列检查它们在任意两行上的交叠是否超过1个1。一种快速检测办法是用 H * H如果非对角元出现大于1的值就说明存在四环。我在一个p31、dc6的配置里确实检测到过四环。解决办法是给基矩阵的位移量加一个随机扰动后重新生成并检测循环几次总能找到无四环的配置。还有一个容易被忽略的问题H矩阵不满秩。规则LDPC的H是dv×p行、dc×p列的矩阵行数不一定等于秩。当行与行之间存在线性相关时实际码率会大于设计码率1-dv/dc也就是说有效信息比特更多但纠错能力变差。我在仿真前会先rank一下如果秩不足就重新找位移量组合宁可多生成几次也要保证H是满行秩。4.2 BP译码不收敛或性能异常的排查如果误码率曲线比理论预期差很多先别急着改算法按下面这个顺序排查。先检查映射关系。BPSK映射0→1、1→-1还是反过来必须和LDPC译码最后的判决逻辑保持一致。最常见的问题就是译码出来的硬判决结果完全翻转误码率一直在0.5或者0.999附近徘徊。解决的办法是先在无噪声条件下跑一帧如果编解码链路自洽任何Eb/N0下都应该能零错误解码。再检查LLR初值。这个我在前面强调过σ²的计算必须严格按 Es/N0 R·Eb/N0 来换算漏乘R会让整条曲线偏移约10log10(R) dB。比如1/2码率漏乘相当于噪声方差小了3dB最后的BER曲线看起来会异常地好而这完全是虚假的。还有数值稳定性问题。tanh和atanh在大LLR输入下容易出现NaN合适的做法是在tanh计算时限制输入绝对值范围或者直接用最小和近似版本。实测中用min-sum替代标准BP在1e-3量级误码率处的性能差距在0.3dB以内但代码简洁、数值稳定前期调试推荐先用它跑通链路。4.3 MATLAB性能优化与仿真效率LDPC仿真的是出了名的“等结果等到怀疑人生”。我的几个经验第一H矩阵和邻居索引一定要初始化好存在cell数组里不要在迭代循环里反复find第二能向量化就不要用双循环特别是校验节点更新用向量化写法可以提速十倍以上第三MATLAB的parfor可以并行扫描多个Eb/N0点四核机器直接省四分之三的等待时间第四中间结果及时保存成mat文件不然跑到一半机器重启就全废了。4.4 运行环境兼容性这套代码不依赖LDPC工具箱里的专用函数只用到了基础MATLAB的sparse、find、gf这些通用能力所以R2016b以后的版本都能跑。唯一需要说明的是如果不想装Communications Toolbox可以参考我用gf写的高斯消元自己用普通的模2运算替换逻辑完全一致。我也建议在脚本开头把随机种子固定下来这样每次跑出来的曲线完全可复现方便调试和对比。我个人在实际操作中的体会是大衍数构造LDPC矩阵这条路虽然不如PEG或ACE算法那样能拿到最优的环分布但它胜在结构确定、代码简洁、可解释性强——所有位移量都能回溯到一条模逆规则写论文时讲构造方法非常顺。另外一个小建议做这种多变量对比实验一定要坚持“只动一个变量”的原则把码率、码长、迭代次数分开研究否则曲线纠缠在一起你根本说不清性能差异到底是谁带来的。后续如果你想继续深挖可以把构造出来的矩阵和PEG算法构造的矩阵放在同一套译码器下对比环长分布和瀑布区位置那会是一篇很完整的实验报告。