极化码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 映射为 +1，1 映射为 -1，经过加噪后接收为 y，LLR 定义为 log(P(x=0|y)/P(x=1|y))，在 AWGN 下近似为 2y/σ²。

function llr = awgn_llr(y, sigma) % 由接收符号计算 LLR，BPSK 映射：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 译码的递归结构为什么容易写错

SC（Successive 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; % LLR<0 判 1 end return; end mid = start_idx + N/2 - 1; % 左子树：f 函数 llr_left = f_llr(llr(start_idx:mid), llr(mid+1:end_idx)); u = sc_decode_recursive(llr_left, frozen_mask, u, start_idx, mid); % 右子树：g 函数，需要左子树的部分和 llr_right = g_llr(llr(start_idx:mid), llr(mid+1:end_idx), u(start_idx:mid)); u = sc_decode_recursive(llr_right, frozen_mask, u, mid+1, 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（或码率 R=K/N）、冻结位设计信噪比。N 决定极化是否充分，N 太小（比如 64 以下）极化效果不明显，BER 曲线下降很慢；N 大到 4096 以上时 SC 译码的串行特性会让 MATLAB 跑得很慢。K 决定码率，码率越高对译码算法要求越高，SC 在高码率下和 SCL 差距明显。设计信噪比决定冻结位选得好不好，我一般取目标 BER 对应的 Eb/N0 再减 1 到 2 dB。

3.3 部分和的回传：SC 译码里最隐蔽的 bug

SC 译码递归返回时，右子树判决完需要把左右子树的结果合并回传，用于上一层 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 译码每个比特只保留一条路径，一旦前面判错后面无法纠正。SCL（Successive Cancellation List）译码保留 L 条候选路径，在每个信息比特处扩展成 2L 条，按路径度量排序后保留最优 L 条，最后选度量最小的路径作为输出。L=8 时，SCL 在中高信噪比下比 SC 有 1 dB 以上的增益，L=32 时接近最大似然。代价是内存和计算量随 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) 是当前比特的 LLR，max(0,-llr) 和 max(0,llr) 分别对应判 0 和判 1 的度量增量。排序后只保留前 L 条。注意路径矩阵是 L×N，每行一条完整路径，这样最后直接取第一行即可。

4.3 SC 与 SCL 的 BER 对比：什么时候值得上 SCL

在 N=1024、码率 0.5、AWGN 信道下，SC 译码在 Eb/N0=2 dB 时 BER 大约在 1e-2 量级，SCL L=8 可以降到 1e-3 以下，L=32 再降半个数量级。但 L 增大带来的仿真时间增长是线性的，L=32 时跑完一条 BER 曲线可能要几十分钟。我的经验是：验证算法正确性用 SC，出论文级曲线用 SCL L=8，需要逼近理论极限再上 L=32。如果只是做链路级验证，SC 加 CRC 辅助也能拿到不错的结果。

5. 极化码 MATLAB 仿真的避坑与排查：五条血泪记录

5.1 BER 曲线不下降反而平了

现象：Eb/N0 从 0 扫到 5 dB，BER 始终在 0.5 附近不动。原因：编码和译码的比特顺序不一致，最常见的是生成矩阵的 Kronecker 顺序和译码递归的索引方向相反。解决：用一个 N=2 的最小例子手工验证，输入 [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 还差

现象：L=8 的 SCL 译码 BER 反而高于 SC。原因：路径度量更新时符号搞反了，把 max(0, llr) 和 max(0, -llr) 写反，导致排序保留的是最差路径。解决：检查度量更新公式，确保 LLR 为正时判 0 的度量增量更小。可以用一个已知正确的短码测试向量验证。

5.4 仿真跑得越来越慢最后卡死

现象：N=1024、L=32 时 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-5，trial 要加到几千，这时候建议用 parfor 并行。

6.2 验证译码正确性的一个技巧

不要一上来就跑 BER 曲线，先用全零码字加噪声验证。全零码字编码后还是全零，接收 LLR 在正信噪比下应该大部分为正，pdecode 输出应该接近全零。如果全零码字都译不对，说明递归结构或部分和有问题。这个技巧能省掉大量调试时间，我每次改完 pdecode 都先跑这个。

6.3 从仿真到落地的判断

极化码 MATLAB 仿真的价值不在于跑出一条好看的曲线，而在于让你能快速回答“换码长会怎样”“换译码算法差多少”“冻结位怎么选”这三个问题。如果你要做硬件实现，MATLAB 仿真的定点化版本可以直接作为 RTL 的参考模型；如果你要做协议级验证，pdecode 的接口可以封装成模块嵌入更大的链路仿真。我的习惯是每改一个参数就存一份结果，文件名带上 N、K、L、设计信噪比，过两周回头看不会认不出哪条曲线是哪次跑的。希望帮到你。

本文还有配套的精品资源，点击获取