极化码MATLAB仿真：从pdecode到SC/SCL译码实战

简介PCode.zip是一份面向通信工程学习者与科研人员的极化码MATLAB仿真程序包聚焦Polar Coding的编解码原理与性能验证。包内共32个m文件压缩包约35KB涵盖码构造、编码、信道合成与解码等核心环节如pdecode、pcencode、initPC、SC_Decoder_ver2等函数并配有systematic编解码与BEC信道相关脚本便于对比不同实现方式。代码注释较为详细可帮助初学者理解生成矩阵、比特翻转与LLR更新等关键步骤的数学含义。已有905人学习下载适合作为课程设计、论文复现或5G极化码入门实践的参考素材通过调整码率、码长与迭代次数读者可观察误码率变化并评估系统性能。1. 极化码 MATLAB 仿真到底在算什么从 PCode.zip 里的 pdecode 说起如果你手上正好有一个叫 PCode.zip 的压缩包里面躺着 Polar Coding 相关的 MATLAB 脚本其中一个文件叫 pdecode那你大概率已经站在了极化码仿真的门口。极化码Polar Code是第一种被严格证明能达到二进制离散无记忆信道容量的构造性编码方案核心思想是信道极化把 N 个独立信道通过递归组合分裂成一部分接近无噪的“好信道”和一部分接近全噪的“坏信道”再把信息比特只放在好信道上。MATLAB 在这个方向上的价值很直接——它让你不用碰硬件就能把编码、调制、信道、译码、误码率统计整条链路跑通改一个参数就能看到 BER 曲线的变化。PCode.zip 这种命名方式通常意味着一个可独立运行的仿真包pdecode 则是其中的译码入口函数。这篇文章面向的是想用 MATLAB 把极化码从公式变成可跑代码的工程师和学生不讲空泛的编码理论史只讲怎么把 pdecode 这类函数用对、参数怎么设、SC 和 SCL 译码在代码层面差在哪、哪些坑会让你的 BER 曲线看起来像玄学。2. 极化码 MATLAB 仿真的最小链路编码、信道与 pdecode 的接口约定2.1 为什么先定接口再写算法极化码仿真翻车最常见的原因不是译码算法写错了而是编码端和译码端对“比特怎么排”的理解不一致。MATLAB 里没有强类型约束一个行向量传进去函数内部按列处理结果全错但代码不报错。所以第一步不是打开 pdecode 看算法而是把整条链路的接口约定写死。我一般会固定四个约定第一所有比特序列统一用行向量长度 N 为 2 的幂第二信息比特位置用逻辑索引向量 frozen_mask 表示true 表示冻结位false 表示信息位第三编码输出是行向量 x长度 N第四pdecode 的输入是接收软信息 LLR 行向量和 frozen_mask输出是估计比特行向量。这四个约定写进每个函数的注释头后面换任何译码算法都不会乱。2.2 编码矩阵的生成不用递归也能构造极化码的生成矩阵是 Kronecker 幂结构MATLAB 里可以用递归或直接矩阵运算两种方式。递归写法直观但 N 大时栈深直接写法更稳。function G polar_gen_matrix(N) % 生成 N×N 极化码生成矩阵N 必须是 2 的幂 % 采用 Kronecker 幂的直接构造避免深递归 n log2(N); if mod(n,1) ~ 0 error(N 必须是 2 的幂); end F [1 0; 1 1]; % 极化核 G 1; for i 1:n G kron(G, F); % Kronecker 积逐层扩展 end % 注意MATLAB 的 kron 顺序与部分教材相反 % 若与你的编码结果对不上先检查这里 end这段代码的逻辑说明F 是 2×2 极化核kron(G,F) 逐层做 Kronecker 积n 层之后得到 N×N 矩阵。参数说明N 必须是 2 的幂否则 log2 不是整数直接报错比后面 BER 异常更容易排查。这里有一个血泪经验——不同教材对 Kronecker 积的左右顺序定义不同如果你的编码结果和参考实现差一个比特反转排列先别怀疑译码回头检查 kron 的参数顺序和是否需要在最后做 bit-reversal 置换。2.3 信道模型与 LLR 计算仿真里最常用的是 AWGN 信道加 BPSK 调制。编码比特 0 映射为 11 映射为 -1经过加噪后接收为 yLLR 定义为 log(P(x0|y)/P(x1|y))在 AWGN 下近似为 2y/σ²。function llr awgn_llr(y, sigma) % 由接收符号计算 LLRBPSK 映射0→1, 1→-1 % sigma 是噪声标准差不是方差 llr 2 * y / (sigma^2); end参数说明sigma 是标准差方差是 sigma^2这两个搞混是 BER 曲线整体平移的经典原因。如果你用 Eb/N0 作为横轴换算关系是 sigma 1/sqrt(2R10^(EbN0/10))其中 R 是码率。我一般会在脚本开头把 Eb/N0 转 sigma 单独写一个函数避免每个脚本里重复推导。2.4 冻结位选择巴氏参数法的最小实现冻结位选哪些位置放信息直接决定性能。工程上最常用的是巴氏参数Bhattacharyya parameter递归计算选巴氏参数最小的 K 个位置作为信息位。function frozen_mask select_frozen(N, K, design_snr_db) % 用巴氏参数递归选择冻结位 % N: 码长, K: 信息比特数, design_snr_db: 设计信噪比 sigma 1/sqrt(2*10^(design_snr_db/10)); z zeros(1, N); z(1) exp(-1/(2*sigma^2)); % 初始巴氏参数 % 递归极化 len 1; while len N for i 1:len z_new_1 2*z(i) - z(i)^2; z_new_2 z(i)^2; z(len i) z_new_2; z(i) z_new_1; end len len * 2; end [~, idx] sort(z, ascend); frozen_mask true(1, N); frozen_mask(idx(1:K)) false; % 巴氏参数最小的 K 个位置放信息 end逻辑说明初始巴氏参数由设计信噪比决定递归过程中每个位置分裂成两个一个变好一个变坏。参数说明design_snr_db 是设计信噪比通常取仿真目标 BER 附近的值取太低会导致高信噪比下性能损失取太高则低信噪比下不划算。注意这里返回的 frozen_mask 中 false 表示信息位和 2.1 的约定一致。3. pdecode 里的 SC 译码递归、LLR 更新与三个必调参数3.1 SC 译码的递归结构为什么容易写错SCSuccessive Cancellation译码是极化码最基础的译码算法按比特顺序逐个判决每个比特的判决依赖前面所有已判决比特。它的核心是 f 和 g 两个 LLR 更新函数以及一个深度为 log2(N) 的递归树。写错的地方通常集中在三处递归的左右子树索引范围、g 函数里部分和partial sum的符号、以及叶子节点判决后部分和的回传。function u_hat pdecode(llr, frozen_mask) % SC 译码主入口 % llr: 接收 LLR 行向量, frozen_mask: 冻结位掩码 N length(llr); u_hat zeros(1, N); u_hat sc_decode_recursive(llr, frozen_mask, u_hat, 1, N); end function u sc_decode_recursive(llr, frozen_mask, u, start_idx, end_idx) % 递归 SC 译码 N end_idx - start_idx 1; if N 1 if frozen_mask(start_idx) u(start_idx) 0; % 冻结位直接判 0 else u(start_idx) llr(start_idx) 0; % LLR0 判 1 end return; end mid start_idx N/2 - 1; % 左子树f 函数 llr_left f_llr(llr(start_idx:mid), llr(mid1:end_idx)); u sc_decode_recursive(llr_left, frozen_mask, u, start_idx, mid); % 右子树g 函数需要左子树的部分和 llr_right g_llr(llr(start_idx:mid), llr(mid1:end_idx), u(start_idx:mid)); u sc_decode_recursive(llr_right, frozen_mask, u, mid1, end_idx); end function out f_llr(a, b) % f 函数符号稳定实现 out sign(a).*sign(b).*min(abs(a), abs(b)); end function out g_llr(a, b, u_left) % g 函数依赖左子树判决结果 out b (1 - 2*u_left) .* a; end逻辑说明f 函数取两个 LLR 的“符号一致取最小绝对值”g 函数根据左子树判决结果决定是加还是减。参数说明f_llr 里用 sign 和 min 组合而不是直接公式是为了避免数值溢出这是工程实现和教科书公式的差别。g_llr 里 u_left 是左子树已判决比特1-2*u_left 把 0/1 映射为 1/-1。3.2 三个必调参数N、K、设计信噪比极化码仿真里真正影响结果的参数就三个码长 N、信息比特数 K或码率 RK/N、冻结位设计信噪比。N 决定极化是否充分N 太小比如 64 以下极化效果不明显BER 曲线下降很慢N 大到 4096 以上时 SC 译码的串行特性会让 MATLAB 跑得很慢。K 决定码率码率越高对译码算法要求越高SC 在高码率下和 SCL 差距明显。设计信噪比决定冻结位选得好不好我一般取目标 BER 对应的 Eb/N0 再减 1 到 2 dB。参数典型取值影响调整建议N256 / 512 / 1024极化充分度与仿真速度先 256 验证再上 1024KN/2 或 N/4码率与性能从 0.5 码率起步设计信噪比目标 Eb/N0 - 1.5 dB冻结位选择质量扫描 0~4 dB 找最优3.3 部分和的回传SC 译码里最隐蔽的 bugSC 译码递归返回时右子树判决完需要把左右子树的结果合并回传用于上一层 g 函数。如果只在叶子节点写 u不在递归返回时做部分和更新译码结果会在 N 大于 2 时全错。常见做法是在递归函数里维护一个局部部分和向量或者在每层返回后重新计算。我一般用后者代码简单但多一次循环N 不大时完全可接受。function ps compute_partial_sum(u, start_idx, end_idx) % 计算部分和用于上层 g 函数 N end_idx - start_idx 1; ps u(start_idx:end_idx); len N; while len 1 len len / 2; for i 1:len ps(i) mod(ps(i) ps(len i), 2); end end ps ps(1); end逻辑说明部分和是极化码递归结构里的模 2 加逐层折半合并。参数说明u 是当前已判决比特start_idx 和 end_idx 界定当前子树范围。注意 mod 2 不能省否则数值会累积成非 0/1。4. 从 SC 到 SCL列表译码在 MATLAB 里的实现要点与性能差距4.1 SCL 为什么值得多写一百行代码SC 译码每个比特只保留一条路径一旦前面判错后面无法纠正。SCLSuccessive Cancellation List译码保留 L 条候选路径在每个信息比特处扩展成 2L 条按路径度量排序后保留最优 L 条最后选度量最小的路径作为输出。L8 时SCL 在中高信噪比下比 SC 有 1 dB 以上的增益L32 时接近最大似然。代价是内存和计算量随 L 线性增长MATLAB 里 L 超过 32 后仿真时间会明显变长。4.2 路径度量的更新与排序SCL 的核心是路径度量path metric的更新。每条路径维护一个度量值扩展时根据判决比特和 LLR 符号是否一致决定度量增加量。function [metrics, paths] scl_extend(metrics, paths, llr, frozen_mask, idx) % SCL 路径扩展与剪枝 % metrics: 当前 L 条路径度量, paths: L×N 路径矩阵 L size(paths, 1); new_metrics zeros(2*L, 1); new_paths zeros(2*L, size(paths, 2)); for l 1:L if frozen_mask(idx) % 冻结位两条扩展路径都判 0度量不变 new_metrics(2*l-1) metrics(l); new_metrics(2*l) metrics(l); new_paths(2*l-1, :) paths(l, :); new_paths(2*l, :) paths(l, :); else % 信息位分别尝试 0 和 1 pm0 metrics(l) max(0, -llr(idx)); pm1 metrics(l) max(0, llr(idx)); new_metrics(2*l-1) pm0; new_metrics(2*l) pm1; new_paths(2*l-1, :) paths(l, :); new_paths(2*l, :) paths(l, :); new_paths(2*l, idx) 1; end end [metrics, sort_idx] sort(new_metrics, ascend); metrics metrics(1:L); paths new_paths(sort_idx(1:L), :); end逻辑说明冻结位不分裂信息位分裂成两条并更新度量。参数说明llr(idx) 是当前比特的 LLRmax(0,-llr) 和 max(0,llr) 分别对应判 0 和判 1 的度量增量。排序后只保留前 L 条。注意路径矩阵是 L×N每行一条完整路径这样最后直接取第一行即可。4.3 SC 与 SCL 的 BER 对比什么时候值得上 SCL在 N1024、码率 0.5、AWGN 信道下SC 译码在 Eb/N02 dB 时 BER 大约在 1e-2 量级SCL L8 可以降到 1e-3 以下L32 再降半个数量级。但 L 增大带来的仿真时间增长是线性的L32 时跑完一条 BER 曲线可能要几十分钟。我的经验是验证算法正确性用 SC出论文级曲线用 SCL L8需要逼近理论极限再上 L32。如果只是做链路级验证SC 加 CRC 辅助也能拿到不错的结果。5. 极化码 MATLAB 仿真的避坑与排查五条血泪记录5.1 BER 曲线不下降反而平了现象Eb/N0 从 0 扫到 5 dBBER 始终在 0.5 附近不动。原因编码和译码的比特顺序不一致最常见的是生成矩阵的 Kronecker 顺序和译码递归的索引方向相反。解决用一个 N2 的最小例子手工验证输入 [0 0]、[0 1]、[1 0]、[1 1] 四种情况看编码输出和译码输出是否一致不一致就查 kron 顺序和 bit-reversal。5.2 高信噪比下 BER 卡在 1e-3 不再下降现象低信噪比正常Eb/N0 超过 3 dB 后 BER 曲线出现错误地板。原因冻结位选择用的设计信噪比太低好信道和坏信道的区分不够部分信息位被放到了实际较差的信道上。解决把设计信噪比提高到目标 BER 对应值再加 1 dB重新生成 frozen_mask错误地板会下移。5.3 SCL 译码结果比 SC 还差现象L8 的 SCL 译码 BER 反而高于 SC。原因路径度量更新时符号搞反了把 max(0, llr) 和 max(0, -llr) 写反导致排序保留的是最差路径。解决检查度量更新公式确保 LLR 为正时判 0 的度量增量更小。可以用一个已知正确的短码测试向量验证。5.4 仿真跑得越来越慢最后卡死现象N1024、L32 时 MATLAB 内存占用飙升最后无响应。原因SCL 路径矩阵每步都做全量复制L×N 矩阵在循环里反复分配。解决预分配路径矩阵用索引操作代替整行复制或者把路径存储改成只存判决比特的稀疏结构。MATLAB 的 profiler 可以定位到具体哪一行在吃内存。5.5 换一个 Eb/N0 定义结果全变现象同样的代码别人论文里 BER 在 2 dB 到 1e-4你的在 2 dB 还是 1e-2。原因Eb/N0 和 Es/N0 混用或者码率没有计入换算。解决统一用 Eb/N0换算 sigma 时显式乘上码率 R并在脚本注释里写清楚定义。这个坑没有技术含量但极其常见我一般会在脚本第一行就写明 Eb/N0 定义。6. 把 pdecode 用出工程价值从 BER 曲线到可复现的仿真脚本6.1 一个可复现的仿真脚本骨架把前面所有函数串起来一个最小可复现脚本长这样% polar_sim_main.m clear; clc; N 1024; K 512; L 8; ebn0_list 0:0.5:3; ber zeros(size(ebn0_list)); for i 1:length(ebn0_list) ebn0 ebn0_list(i); R K / N; sigma 1 / sqrt(2 * R * 10^(ebn0/10)); frozen_mask select_frozen(N, K, ebn0); G polar_gen_matrix(N); err 0; total 0; for trial 1:200 info randi([0 1], 1, K); u zeros(1, N); u(~frozen_mask) info; x mod(u * G, 2); y (1 - 2*x) sigma * randn(1, N); llr awgn_llr(y, sigma); u_hat pdecode(llr, frozen_mask); err err sum(u_hat(~frozen_mask) ~ info); total total K; end ber(i) err / total; fprintf(Eb/N0%.1f dB, BER%.4e\n, ebn0, ber(i)); end semilogy(ebn0_list, ber, o-); grid on; xlabel(Eb/N0 (dB)); ylabel(BER);逻辑说明外层扫 Eb/N0内层蒙特卡洛跑 200 帧统计信息位错误率。参数说明trial 数决定 BER 的统计置信度1e-3 量级至少需要 1e5 个信息比特200 帧×512 比特刚好够。如果要看 1e-5trial 要加到几千这时候建议用 parfor 并行。6.2 验证译码正确性的一个技巧不要一上来就跑 BER 曲线先用全零码字加噪声验证。全零码字编码后还是全零接收 LLR 在正信噪比下应该大部分为正pdecode 输出应该接近全零。如果全零码字都译不对说明递归结构或部分和有问题。这个技巧能省掉大量调试时间我每次改完 pdecode 都先跑这个。6.3 从仿真到落地的判断极化码 MATLAB 仿真的价值不在于跑出一条好看的曲线而在于让你能快速回答“换码长会怎样”“换译码算法差多少”“冻结位怎么选”这三个问题。如果你要做硬件实现MATLAB 仿真的定点化版本可以直接作为 RTL 的参考模型如果你要做协议级验证pdecode 的接口可以封装成模块嵌入更大的链路仿真。我的习惯是每改一个参数就存一份结果文件名带上 N、K、L、设计信噪比过两周回头看不会认不出哪条曲线是哪次跑的。希望帮到你。本文还有配套的精品资源点击获取