企业官网建设流程全解析

最近在折腾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下的乘法逆元，简称模逆。

比如p=7，a=3，因为3×5=15≡1 (mod 7)，所以5就是3在模7下的逆。大衍求一术本质上是用辗转相除的思路去解这个线性同余式，而这个思路和今天扩展欧几里得算法做的事情一模一样。换句话说，大衍求一术就是扩展欧几里得算法在古代的雏形，不是什么神秘力量，而是一套确定性的数论算法。

在MATLAB里求模逆非常方便，扩展欧几里得的系数直接由gcd函数返回：

function inva = modInv(a, p) % 基于扩展欧几里得算法求 a 在模 p 下的逆元 % 要求 gcd(a,p)=1，p通常取素数 [g, x, ~] = gcd(a, p); if g ~= 1 error('a和p不互素，无法求模逆'); end inva = mod(x, p); end

这行代码背后对应的就是大衍求一术中的计算流程，你甚至可以把这个函数直接改名叫“daYanQiuYiShu”，仿真代码里还能带着一点考据的趣味。

1.3 用模逆构造准循环稀疏校验矩阵的思路

有了模逆这个工具，怎么把它变成校验矩阵？我是按准循环LDPC（QC-LDPC）的思路来构造的。

先定三个参数：列重dv、行重dc、子矩阵大小p。设计一个dv行dc列的基矩阵B，基矩阵里的每一个元素都是一个模p下的数，代表对应位置上循环移位矩阵的位移量。然后把基矩阵的每个元素替换成一个p×p的循环移位单位阵，也就是单位矩阵循环右移若干位。这样最终得到的H矩阵就是(dv×p)行、(dc×p)列的准循环稀疏矩阵，每个子块里只有一个1，整块矩阵的行重固定为dc、列重固定为dv。

基矩阵里每个位置的“移位量”怎么定？这里就是大衍数登场的地方。一种简单且有效的方案是：对基矩阵第i行第j列的位置，取 a_{i,j} = (i+j) 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)*p+1 : i*p, (j-1)*p+1 : j*p) = block; end end H = sparse(H); % 转稀疏，后续译码提速 end

用这套构造方法，H矩阵天然稀疏，因为每个p×p子块只有一个非零元，整个矩阵的非零元密度就是1/p，p越大矩阵越稀疏。而且这种准循环结构对硬件实现特别友好，消息传递译码时可以用循环移位寄存器完成节点信息更新，这也是5G标准里LDPC采用准循环结构的原因。

2. MATLAB仿真框架搭建：从编码到BP译码

2.1 整体仿真链路怎么设计

整个仿真链路说白了就是通信系统中最经典的一段：信源产生二进制随机消息 → LDPC编码 → BPSK调制 → 加高斯白噪声AWGN → 软解调得到对数似然比LLR → BP迭代译码 → 统计误码率和误帧率。每个环节都不算难，但串联起来之后，很多细节会决定你仿真结果靠不靠谱。

重点说两个容易踩坑的地方。

第一个是信噪比定义。LDPC性能曲线一般用Eb/N0（每比特能量与噪声功率谱密度之比）来画，而不是直接用SNR。BPSK调制的符号能量Es=1时，由码率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: 信道初始LLR，1×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 仿真参数设计

为了让对比有意义，我做实验时坚持“只动一个变量”的原则。三组实验设计如下：

参数对应关系：码率通过基矩阵列重dc调节，dv固定为3，当dc=4时设计码率R=1-3/4=1/4，dc=6时R=1/2，dc=12时R=3/4。码长通过子矩阵大小p调节，对于dc=6的1/2码率，p=128对应码长n=768？这里注意，我实际仿真时为了码长精确等于2的幂次，会让dc=6，然后选n=1024，p=1024/6不是整数。所以更稳妥的办法是让dc=4或dc=8这样好整除的取值，或者允许码长为dc×p的乘积形式。我这里的实验三码长系列2064/1024/512/256，采用dc=8、dv=4的组合来保证R=1/2且整除，这也是一般QC-LDPC常用的(a,b)规则。用纯大衍基矩阵构造时，大家可以根据实际需要灵活调整p。

Eb/N0扫描范围视码率而定：1/2码率从0dB扫到6dB，步长0.5dB；每个点最少统计200帧错误。

3.2 迭代次数：从欠收敛到饱和

先看迭代次数的影响。在1/2码率、码长1024条件下，Eb/N0=3dB附近，5次迭代的误码率大约在1e-2量级，10次迭代降到2e-3左右，20次迭代到4e-4，而50次迭代相比20次提升非常有限，大概只到3e-4。这个现象很容易理解：消息传递译码在前几轮迭代时信息快速扩散，比特置信度快速提高，但到了一定轮数之后，环的存在让外信息趋于相关，再迭代也很难带来额外增益，就进入了饱和平台。

做这个对比实验的实际意义在于：工程上迭代次数直接决定译码延迟和功耗。如果系统要求1e-3的误码率，20次迭代已经够用，没必要上50次；但如果目标是1e-5以下，20次可能不够，需要结合更长的码长来换性能。所以我在报告里会把“迭代-误码率-复杂度”三者的关系一起分析，只谈误码率不谈复杂度说服力不足。

3.3 码率：以频谱效率换功率增益

码率的影响非常直观。在码长1024、迭代20次、BER=1e-4目标下，1/4码率需要的Eb/N0大约是1.2dB左右，1/2码率大约2.6dB，3/4码率大约4.1dB。也就是说，码率从1/2降到1/4，能换来大约1.4dB的编码增益，但代价是同样信息比特要占用两倍的带宽资源——这就是通信系统里经典的“带宽换功率”折中。

还有一个现象值得注意：3/4码率在低信噪比区间误码率下降得没有低码率那么陡，瀑布区的斜率更缓。这是因为高码率码字中的冗余校验比特少，纠错能力天然弱，在噪声较大时更容易出现不可纠正的错误。如果目标工作区在低信噪比，选低码率会省心很多。

3.4 码长：逼近香农限的筹码

码长的对比结果说明LDPC“码长越长性能越好”确实是成立的。在1/2码率、20次迭代下，BER=1e-3处，码长256的码字需要约4.5dB，512码长约3.8dB，1024码长约3.2dB，2048码长降到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不互素”。我在调试时发现这是因为基矩阵索引a=mod(i+j,p)取到了0，而0在模p下没有逆。解决办法是遇到0时映射到p-1，或者重新设计规则，比如 a=mod(i+j, p)+1，保证a在1到p-1范围内。

另一个大坑是四环。大衍构造虽然让多数位移量差异很大，但并不能绝对保证没有四环。判断方法很简单：随机抽两列，检查它们在任意两行上的交叠是否超过1个1。一种快速检测办法是用 H' * H，如果非对角元出现大于1的值，就说明存在四环。我在一个p=31、dc=6的配置里确实检测到过四环。解决办法是给基矩阵的位移量加一个随机扰动后重新生成并检测，循环几次总能找到无四环的配置。

还有一个容易被忽略的问题：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算法构造的矩阵放在同一套译码器下对比环长分布和瀑布区位置，那会是一篇很完整的实验报告。