极化码MATLAB仿真：编码、SC译码与误码率统计实战

简介PCode.zip 是一份面向通信工程学习者与科研人员的极化码 MATLAB 仿真源码包聚焦 Erdal Arikan 提出的极化码编解码原理与性能验证适合具备一定 MATLAB 与信道编码基础、希望深入理解 5G 极化码实现细节的读者。压缩包共 32 个文件全部为 .m 脚本整体约 35KB涵盖码构造、编码、信道合成与解码等核心模块如 pdecode、pcencode、initPC、SC_Decoder_ver2 等函数并配有 MonteCarlo 仿真与误码率绘图脚本便于观察码率、码长、迭代次数对系统性能的影响。源码注释较为详细有助于理解比特翻转、LLR 更新与 BEC 信道处理等关键步骤。目前已有 905 人学习下载可作为极化码理论学习与仿真实验的参考工具。1. 从一份 PCode.zip 说起极化码 MATLAB 仿真到底能跑出什么如果你正在做极化码Polar Coding相关的课程设计、论文复现或者算法预研大概率会遇到一个很现实的问题理论推导看懂了SC 译码的递归公式也背下来了但真要在 MATLAB 里从零搭一套能跑通编码、信道传输、译码、误码率统计的完整链路光是索引对齐和递归边界就能耗掉两三天。PCode.zip 就是冲着这个场景来的——它把极化码的编码构造、SC 译码pdecode以及配套的仿真脚本打包在一起解压后直接是 MATLAB 工程文件不需要额外编译也不依赖通信工具箱里的冷门函数。适合谁适合已经学过极化码基本原理、需要快速拿到可运行基线代码做对比实验的人也适合刚接触极化码、想通过读一份能跑通的代码来反向理解 SC 译码递归逻辑的新手。但要注意这份资源不是“一键出论文图”的黑匣子参数配置和信道模型需要你自己根据实验目标调整下面我会把每个环节拆开讲。2. 极化码编码链路拆解从信道极化到码字生成2.1 极化码的核心机制与 MATLAB 实现映射极化码的理论基础是信道极化通过递归的信道合并与分裂操作让一部分比特信道趋于无噪声、另一部分趋于全噪声然后把信息比特放在“好”信道上。在 MATLAB 里这个过程的实现通常围绕三个核心变量展开码长 N、信息比特数 K、以及冻结比特位置集合。PCode.zip 里的编码脚本一般会先根据高斯近似或巴氏参数法计算各子信道的可靠性再排序选出 K 个最可靠的索引放信息比特其余置零作为冻结比特。这里有个容易被忽略的点极化码的编码本质上是线性分组码生成矩阵是核矩阵 F 的克罗内克积。MATLAB 里如果直接用kron函数逐级构造当 N 到 1024 以上时内存占用会明显上升。常见做法是用递归函数代替显式矩阵构造或者利用极化码编码的蝶形结构做原地运算。PCode.zip 里的实现我翻了一下用的是递归编码方式对 N1024 这种量级比较友好。% 极化码递归编码核心逻辑示意 function x polar_encode(u, N) % u: 输入比特序列长度 N % N: 码长必须是 2 的幂 if N 1 x u; return; end % 将输入分为前后两半 u1 u(1:N/2); u2 u(N/21:N); % 递归编码前半部分 x1 polar_encode(u1, N/2); x2 polar_encode(u2, N/2); % 蝶形合并前半为模二加后半直接传递 x [mod(x1 x2, 2), x2]; end上面这段代码展示的是极化码编码的递归骨架。u是经过信息位映射后的完整比特序列长度 N 必须是 2 的幂。递归终止条件是 N1此时直接返回。合并步骤里mod(x1 x2, 2)对应极化码核矩阵 F 的上支路异或操作x2对应下支路直通。参数 N 决定了递归深度N1024 时递归深度为 10 层每层做一次长度为 N 的向量操作整体复杂度是 O(N log N)。如果你把 N 设成非 2 的幂递归会在 N/2 处产生小数索引MATLAB 会直接报维度错误——这是新手最容易踩的坑之一。2.2 冻结比特位置选择三种常用构造法对比冻结比特集合的选取直接决定极化码在有限码长下的性能。PCode.zip 里默认用的可能是巴氏参数法Bhattacharyya parameter但实际做仿真时高斯近似GA和极化权重PW也很常见。三者的区别在于构造方法适用信道计算复杂度有限码长性能巴氏参数法BEC 信道低在 BEC 下最优AWGN 下偏差大高斯近似AWGN 信道中中短码长表现稳定极化权重AWGN 信道极低无需信道参数工程常用如果你做的是 AWGN 信道下的仿真建议把构造方法切到高斯近似或极化权重。PCode.zip 的脚本里通常有一个construct_method变量或者类似的配置项改成对应字符串即可。注意切换构造方法后冻结比特集合会变译码端的冻结位置必须同步更新否则译码器会把信息比特当冻结比特处理误码率直接飙到 0.5 左右——这不是代码 bug是配置不一致。% 高斯近似构造冻结比特集合简化版 function frozen_bits ga_construct(N, K, snr_dB) % N: 码长, K: 信息比特数, snr_dB: 设计信噪比 sigma 10^(-snr_dB/20); % 初始化各子信道对数似然比均值 L 2 / sigma^2; % 递归计算各分裂信道的可靠性此处省略递归细节 % 最终得到可靠性排序 reliability_rank reliability_rank compute_reliability(N, L); % 选取可靠性最低的 N-K 个位置作为冻结比特 [~, sorted_idx] sort(reliability_rank, ascend); frozen_bits sorted_idx(1:N-K); end这段代码的关键参数是snr_dB它决定了构造时假设的信道质量。设计信噪比和仿真信噪比不需要完全一致但偏差太大会导致构造失配性能损失在 0.5 dB 以内可以接受超过 1 dB 就说明构造方法需要换。compute_reliability是递归计算各子信道可靠性的函数具体实现可以用高斯近似的递归公式也可以用查表法加速。PCode.zip 里如果自带这个函数直接调用即可如果没有需要自己补上。3. SC 译码的 pdecode 实现递归、判决与回溯3.1 SC 译码的递归结构在 MATLAB 里怎么落地SCSuccessive Cancellation译码是极化码最基础的译码算法核心思想是逐比特判决每个比特的判决依赖前面所有已判决比特。在 MATLAB 里实现 SC 译码最直观的方式是递归函数但递归深度随码长对数增长N1024 时递归深度为 10MATLAB 默认递归限制通常够用但如果 N 继续增大到 4096 以上可能需要调整递归深度或者改成迭代实现。PCode.zip 里的pdecode函数大概率是递归版本。SC 译码的递归过程分两个阶段先递归计算左子节点的 LLR对数似然比判决后再递归计算右子节点。左子节点的 LLR 计算用f函数右子节点用g函数。这两个函数的数值稳定性很关键尤其是f函数在 LLR 绝对值很大时容易出现数值溢出。% SC 译码中的 f 和 g 函数数值稳定版 function out f_func(a, b) % f(a,b) 2*atanh(tanh(a/2)*tanh(b/2)) % 数值稳定实现避免直接计算 tanh 溢出 out sign(a).*sign(b).*min(abs(a), abs(b)) ... log(1 exp(-abs(ab))) - log(1 exp(-abs(a-b))); end function out g_func(a, b, u) % g(a,b,u) b (1-2u)*a out b (1 - 2*u) * a; endf_func里用的是 min-sum 近似加修正项比直接套公式稳定得多。sign(a).*sign(b).*min(abs(a),abs(b))是 min-sum 的主项后面两个 log 项是修正。如果你追求速度、对精度要求不高可以只保留 min-sum 主项性能损失大约 0.2~0.3 dB。g_func里的u是左子节点的判决结果取值 0 或 1(1-2*u)把 0 映射为 1、1 映射为 -1这是极化码译码里的标准操作。3.2 译码流程中的比特判决与冻结位处理SC 译码的判决逻辑是如果当前比特索引属于冻结比特集合直接判为 0否则根据 LLR 的符号判决LLR ≥ 0 判 0LLR 0 判 1。这个逻辑本身很简单但在递归回溯时判决顺序必须严格从左到右不能并行化——这是 SC 译码的固有串行特性。PCode.zip 的pdecode函数里冻结比特集合通常作为输入参数传入。这里有个实操细节冻结比特集合的索引方式要和编码端一致。有的实现用 0-based 索引有的用 1-basedMATLAB 默认是 1-based但如果你从 Python 或 C 代码移植过来很容易在这里翻车。我一般会在译码函数入口加一行断言检查冻结比特索引的最大值是否等于 N如果不等就直接报错省得后面误码率异常时到处找原因。% SC 译码主循环简化示意 function u_hat sc_decode(llr, frozen_bits, N) % llr: 接收到的对数似然比序列 % frozen_bits: 冻结比特索引集合1-based % N: 码长 u_hat zeros(1, N); % 递归译码 u_hat recursive_decode(llr, frozen_bits, N, u_hat); end function u recursive_decode(llr, frozen_bits, N, u) if N 1 idx length(u) 1; % 当前判决位置 if ismember(idx, frozen_bits) u(idx) 0; else u(idx) (llr 0); % LLR0 判 1 end return; end % 计算左子节点 LLR a llr(1:N/2); b llr(N/21:N); llr_left f_func(a, b); % 递归译码左子节点 u recursive_decode(llr_left, frozen_bits, N/2, u); % 计算右子节点 LLR u_left u(end-N/21:end); llr_right g_func(a, b, u_left); % 递归译码右子节点 u recursive_decode(llr_right, frozen_bits, N/2, u); end这段代码里ismember(idx, frozen_bits)是冻结比特判断idx的计算依赖length(u)这要求递归过程中u的更新顺序严格从左到右。llr_left和llr_right的计算分别对应 SC 译码的两个递归分支。注意u_left的提取位置u(end-N/21:end)取的是最近判决的左半部分比特这个索引如果写错译码结果会完全乱掉。我见过有人把end-N/21写成end-N/2结果误码率一直卡在 0.3 下不去排查了半天才发现是索引差了一位。4. 仿真脚本配置与误码率统计参数怎么设才不翻车4.1 信噪比扫描与蒙特卡洛仿真次数极化码仿真的核心输出是误码率BER或误块率BLER随信噪比变化的曲线。PCode.zip 里通常有一个主脚本负责设置信噪比范围、每个信噪比下的仿真帧数、以及调用编码和译码函数。这里有两个参数直接决定仿真结果的平滑度和可信度信噪比步长和每点帧数。信噪比步长一般取 0.5 dB 或 1 dB步长太大曲线拐点看不清楚步长太小仿真时间成倍增加。每点帧数取决于目标误码率量级如果目标 BER 在 10^-3 左右每点至少 10^4 帧才能保证统计置信度如果目标 BER 到 10^-5帧数要加到 10^6 量级。PCode.zip 默认的帧数可能偏少跑出来的曲线毛刺比较多建议根据实际需求调大。% 仿真主循环信噪比扫描 蒙特卡洛 snr_range 0:0.5:4; % 信噪比范围步长 0.5 dB num_frames 1e4; % 每个信噪比下的帧数 N 1024; % 码长 K 512; % 信息比特数 ber zeros(size(snr_range)); for snr_idx 1:length(snr_range) snr snr_range(snr_idx); error_count 0; total_bits 0; for frame 1:num_frames % 生成随机信息比特 info_bits randi([0 1], 1, K); % 编码含冻结比特插入 tx_bits polar_encode_with_frozen(info_bits, N, K, snr); % BPSK 调制 AWGN sigma 10^(-snr/20); rx_signal (1 - 2*tx_bits) sigma * randn(1, N); % 计算 LLR llr 2 * rx_signal / sigma^2; % SC 译码 u_hat sc_decode(llr, frozen_bits, N); % 提取信息比特并统计误码 info_hat u_hat(info_indices); error_count error_count sum(info_hat ~ info_bits); total_bits total_bits K; end ber(snr_idx) error_count / total_bits; fprintf(SNR %.1f dB, BER %.4e\n, snr, ber(snr_idx)); end这段代码里sigma 10^(-snr/20)是 BPSK 调制下噪声标准差的换算公式llr 2 * rx_signal / sigma^2是 AWGN 信道下 LLR 的标准计算方式。info_indices是信息比特在完整码字中的索引集合必须和编码端的映射一致。num_frames设成 1e4 时SNR4 dB 附近如果 BER 已经降到 10^-5 以下统计到的错误数可能只有个位数曲线尾部会抖动这时候要么加大帧数要么在图上标注置信区间。4.2 冻结比特集合与信息位索引的一致性校验前面提过冻结比特集合不一致会导致译码崩溃这里展开说具体怎么校验。编码端和信息位索引的映射关系通常有两种实现方式一种是先构造完整比特序列信息位 冻结位再整体编码另一种是先编码信息比特再在输出端插入冻结位对应的零。PCode.zip 用的是哪种需要看polar_encode_with_frozen函数的实现。不管哪种方式校验方法是一样的编码前记录信息比特在完整序列中的位置索引info_indices译码后从u_hat里按同样的索引提取信息比特。如果info_indices和冻结比特集合frozen_bits满足union(info_indices, frozen_bits)恰好覆盖 1 到 N 且无交集说明配置正确。我一般会在脚本开头加一段校验代码% 校验信息位与冻结位配置一致性 assert(length(info_indices) length(frozen_bits) N, ... 信息位与冻结位总数不等于码长 N); assert(isempty(intersect(info_indices, frozen_bits)), ... 信息位与冻结位存在重叠索引);这两行断言花不了几毫秒但能省掉后面几小时的排查时间。血泪经验有一次跑仿真BER 曲线在低信噪比下正常高信噪比下反而上升最后发现是冻结比特集合里混入了一个信息位索引导致译码端把信息比特强制判零高信噪比下错误反而更明显。5. 避坑与排查极化码 MATLAB 仿真里最容易翻车的五个点5.1 现象误码率曲线在高信噪比下不降反升原因冻结比特集合与信息位索引不一致译码端把信息比特当冻结比特处理信噪比越高LLR 越可靠但强制判零导致的错误越明显。解决在编码和译码之间加一致性断言检查info_indices和frozen_bits的并集是否等于 1:N 且无交集。5.2 现象N1024 时仿真跑不动内存爆掉原因显式构造生成矩阵kron多次后矩阵维度爆炸N1024 时生成矩阵是 1024×1024再往上做克罗内克积内存占用呈指数增长。解决改用递归编码函数或者利用蝶形结构做原地运算避免显式存储生成矩阵。5.3 现象SC 译码结果全为 0 或全为 1原因LLR 符号判决逻辑写反或者f_func的数值稳定实现里符号处理有误。解决检查u(idx) (llr 0)这一行确认 LLR ≥ 0 判 0、LLR 0 判 1同时用简单测试用例N2, K1手动验证f_func输出。5.4 现象仿真曲线毛刺严重无法判断性能趋势原因每点帧数不足高信噪比下统计到的错误数太少BER 估计方差大。解决将num_frames提高到 1e5 或 1e6或者在图上标注每个点的错误数错误数少于 100 的点视为不可信。5.5 现象切换构造方法后性能反而变差原因构造时假设的信道模型与仿真信道不匹配比如用巴氏参数法构造却在 AWGN 信道下仿真。解决AWGN 信道下优先用高斯近似或极化权重构造设计信噪比与仿真信噪比偏差控制在 1 dB 以内。6. 进阶技巧把 SC 译码改成 SCL 并验证增益SC 译码的局限在于逐比特硬判决没有利用路径度量做候选保留。SCLSuccessive Cancellation List译码保留 L 条候选路径在译码结束时选路径度量最小的作为输出在相同码长下能带来 0.5~1 dB 的增益。PCode.zip 里的pdecode是 SC 版本但改成 SCL 并不需要重写整个递归结构核心改动是在每层递归返回时维护一个路径列表而不是单条路径。具体做法把u从一维向量改成 L×N 的矩阵每行是一条候选路径在冻结比特位置所有路径都判 0路径度量不变在信息比特位置每条路径分裂成两条判 0 和判 1计算路径度量增量然后按度量排序保留前 L 条。路径度量增量用 LLR 的绝对值加权判 0 时增量为log(1exp(-llr))判 1 时增量为log(1exp(llr))实际实现里常用近似max(0, -llr)和max(0, llr)来降低计算量。% SCL 译码路径度量更新简化示意 function [paths, metrics] scl_update(paths, metrics, llr, L, is_frozen) if is_frozen % 冻结比特所有路径判 0度量不变 return; end % 信息比特每条路径分裂为两条 num_paths size(paths, 1); new_paths zeros(2*num_paths, size(paths, 2)); new_metrics zeros(2*num_paths, 1); for i 1:num_paths % 判 0 分支 new_paths(2*i-1, :) [paths(i, :), 0]; new_metrics(2*i-1) metrics(i) max(0, -llr); % 判 1 分支 new_paths(2*i, :) [paths(i, :), 1]; new_metrics(2*i) metrics(i) max(0, llr); end % 按度量排序保留前 L 条 [sorted_metrics, sort_idx] sort(new_metrics, ascend); keep sort_idx(1:min(L, length(sort_idx))); paths new_paths(keep, :); metrics sorted_metrics(1:min(L, length(sort_idx))); end这段代码里L是列表大小典型取值 4、8、16L 越大性能越好但计算量线性增长。max(0, -llr)和max(0, llr)是路径度量增量的近似计算比精确的 log-sum-exp 快很多性能损失可以忽略。注意paths的列数在每次更新后增加 1最终长度等于 N。SCL 译码结束后选度量最小的路径按info_indices提取信息比特。验证 SCL 增益的方法固定 N、K、构造方法和信噪比范围分别跑 SC 和 SCLL8对比 BER 曲线。如果 SCL 在 BER10^-3 处比 SC 有 0.5 dB 以上的增益说明实现正确如果增益不明显甚至更差检查路径度量更新是否在冻结比特位置错误地分裂了路径或者排序保留逻辑写反了。从那以后我每次改译码算法都会先用 N8、K4 的小规模用例手动验证每一步的路径和度量确认无误后再上大码长跑仿真。希望帮到你。本文还有配套的精品资源点击获取