极化码Matlab仿真从原理到实践：编码、SC/SCL译码与BER曲线

简介：这份资源是面向通信工程、信息论学习者与科研人员的极化码基础仿真代码包，围绕Erdal Arikan提出的信道极化理论，提供从编码、调制、加噪到解码与误码率评估的完整Matlab实现，适合希望动手验证极化码原理、复现经典算法或开展课程设计的中高级读者。压缩包共23个文件，约6MB，以20个m脚本为核心，覆盖编码、SC解码、蝶形运算、冻结比特初始化、信道容量计算等模块，另含1个pdf学位论文、1个mat数据文件和1个md说明文档，便于理论对照与代码调试。目前已有225人学习下载。读者可借助这些脚本理解极化码的构造与译码流程，结合外文论文深入掌握SC及列表译码思路，并在此基础上修改参数、扩展算法或探索5G控制信道等应用场景。

1. 从一份 Polar-Codes-Software-Matlab 压缩包说起：极化码到底能跑出什么

很多人第一次接触极化码，是在 5G 控制信道标准里看到「Polar Code」这个词，然后顺手搜到一个叫Polar-Codes-Software-Matlab--master.zip的压缩包，解压一看全是.m文件，不知道从哪下手。我当年也是这样，把文件夹翻了三遍才搞明白：这套东西的核心价值，是让你在本地用 Matlab 把极化码从「信道极化」到「SC 译码」整条链路跑通，而不是只停留在公式推导。它解决的是「我懂极化码原理，但不知道怎么用代码验证误码率」这个断层。适合两类人：一是通信、电子信息方向的研究生和工程师，需要快速搭一个可改参数的仿真平台；二是准备做毕设或课程设计的人，想找一个能跑、能改、能出 BER 曲线的极化码 Matlab 实现。这一章先把「极化码 + Matlab」这件事讲清楚，后面再拆怎么跑、怎么调、怎么避坑。

极化码的理论根基是「信道极化」：把 N 个相同的二进制离散无记忆信道，通过一个线性变换，变成一组容量两极分化的比特信道——一部分趋近于无噪，一部分趋近于全噪。然后只把信息比特放在那些「好信道」上，坏信道上放冻结比特（通常全 0）。这个思路听起来简单，但真正落地时，你会遇到三个绕不开的问题：怎么算每个比特信道的可靠性、怎么构造生成矩阵、怎么译码。Matlab 在这件事上的优势是矩阵操作天然贴合极化码的递归结构，kronecker积、reshape、mod这些函数几行就能把编码器写出来，比 C 语言少很多指针和内存管理的麻烦。所以Polar-Codes-Software-Matlab这类包，本质上是一个「教学 + 验证」工具，不是工业级编码器，但足够你理解每一步在干什么。

2. 极化码在 Matlab 里的最小可跑链路：编码、信道、译码三段拆解

2.1 先搞懂极化码的三个核心参数：N、K、A

在动手写代码之前，必须把三个参数钉死，否则后面所有仿真都是玄学。N 是码长，必须是 2 的幂，常见 256、512、1024、2048；K 是信息比特数，决定码率 R = K/N；A 是信息位集合，也就是「好信道」的下标集合，大小等于 K。这三个参数一旦确定，编码器和译码器的结构就定了。很多新手翻车，是因为把 A 当成随便取的，结果译码端和编码端用的 A 不一致，BER 直接 0.5，还以为是算法错了。

A 怎么来？常见做法是用高斯近似（Gaussian Approximation, GA）或者密度进化（Density Evolution, DE）算每个比特信道的错误概率，然后从小到大排序，取前 K 个作为信息位。Matlab 里如果不想自己实现 DE，可以用简化版：在 AWGN 信道下用 GA 递推，核心是计算对数似然比的均值和方差。下面这段代码是我常用的 GA 计算 A 的最小实现，假设信噪比用 Eb/N0 表示，先转成设计信噪比。

function A = polar_ga_construct(N, K, EbN0_dB) % 高斯近似构造极化码信息位集合 % N: 码长, K: 信息比特数, EbN0_dB: 设计信噪比(dB) % 返回 A: 信息位下标(列向量), 按可靠性从高到低排列 % 将 Eb/N0 转为噪声方差, 假设 BPSK 调制, 码率 R=K/N R = K / N; EbN0 = 10^(EbN0_dB / 10); sigma2 = 1 / (2 * R * EbN0); % 每个符号的噪声方差 % 初始化: 第0层只有一个信道, LLR 均值 = 2/sigma2 L = zeros(N, 1); L(1) = 2 / sigma2; % 逐层递推, 每层信道数翻倍 n = log2(N); for level = 1:n half = 2^(level - 1); for i = 1:half % 公式: 两个子信道 LLR 均值满足递推关系 % 这里用简化的 GA 递推, 实际工程中常用查表或近似函数 L_new1 = phi_inv(1 - (1 - phi(L(i)))^2); L_new2 = 2 * L(i); L(i) = L_new1; L(half + i) = L_new2; end end % 按 LLR 均值从大到小排序, 取前 K 个作为信息位 [~, idx] = sort(L, 'descend'); A = sort(idx(1:K)); end function y = phi(x) % 高斯近似中的 phi 函数近似 if x < 0.001 y = 1; elseif x > 20 y = 0; else y = exp(-0.4527 * x^0.86 + 0.0218); end end function x = phi_inv(y) % phi 的反函数近似, 用于递推 if y >= 1 x = 0; elseif y <= 0 x = 20; else x = (-log(y) / 0.4527)^(1/0.86); end end

这段代码的逻辑是：从根信道出发，每一层把当前信道的 LLR 均值分裂成两个，一个变好一个变差，递推 n 层后得到 N 个信道的可靠性度量。phi和phi_inv是 GA 里常用的近似函数，精度够用但不是唯一写法。参数说明：EbN0_dB是设计信噪比，一般取你目标工作点附近的值，比如 2~4 dB；A返回的是排序后的下标，编码时把信息比特放在这些位置上。注意A必须同时传给编码和译码，否则必翻车。

2.2 编码器：三行核心代码和生成矩阵的构造

极化码编码公式是x = u * G，其中G = F^{⊗n}，F = [1 0; 1 1]，⊗是克罗内克积。Matlab 里构造G最直接的方式是递归做kron，但 N 大了以后G是 N×N 矩阵，1024 码长就是 100 万个数，内存吃得消但没必要。实际编码可以用「蝶形」结构，复杂度 O(N log N)。不过为了先跑通，我建议先用矩阵乘法验证正确性，再换快速算法。

function x = polar_encode(u, A, N) % 极化码编码 % u: 信息比特(长度 K), A: 信息位下标, N: 码长 % 返回 x: 编码后比特(长度 N) % 步骤1: 把信息比特放到信息位, 其余置0(冻结比特) u_full = zeros(N, 1); u_full(A) = u(:); % 步骤2: 构造生成矩阵 G = F^{⊗n}, 用递归 kron n = log2(N); F = [1 0; 1 1]; G = 1; for i = 1:n G = kron(G, F); end % 步骤3: 编码并做模2运算 x = mod(G' * u_full, 2); % 注意这里用 G' 因为 Matlab 列向量习惯 end

逻辑说明：先把信息比特映射到u_full的A位置上，冻结比特保持 0；然后构造生成矩阵，kron递归 n 次得到F^{⊗n}；最后矩阵乘加模 2。参数说明：u长度必须等于length(A)，A是列向量或行向量都行，但u_full(A)要求A是正整数下标。这里有个坑：Matlab 的kron顺序和某些教材定义的F^{⊗n}可能差一个转置，如果你发现译码结果全错，先把G转置试一下。另外mod(G' * u_full, 2)在 N=1024 时矩阵乘是 1024×1024，跑一次几毫秒，仿真几万帧就慢了，后面可以换成蝶形快速编码。

2.3 SC 译码：递归函数怎么写、LLR 怎么传

SC（Successive Cancellation）译码是极化码的经典译码算法，核心是递归的f和g函数。f函数处理两个 LLR 的软信息合并，g函数利用已译出的比特做干扰消除。Matlab 里用递归实现最直观，但要注意递归深度是 log2(N)，N=1024 时深度 10，不会爆栈。下面是一个可跑的 SC 译码核心。

function u_hat = polar_sc_decode(y, A, N) % SC 译码 % y: 接收到的软信息(LLR), 长度 N % A: 信息位下标, N: 码长 % 返回 u_hat: 译码比特(长度 N) n = log2(N); u_hat = zeros(N, 1); % 递归译码, 从第 n 层开始 u_hat = sc_recursive(y, n, A, u_hat, 1, N); end function u = sc_recursive(y, n, A, u, start_idx, N) % 递归 SC 译码 if n == 0 % 叶子节点: 根据 LLR 符号判决 if y >= 0 u(start_idx) = 0; else u(start_idx) = 1; end return; end half = 2^(n - 1); % 计算 f 函数: 左子节点 LLR y_left = f_llr(y(1:half), y(half+1:2*half)); % 递归译码左子树 u = sc_recursive(y_left, n-1, A, u, start_idx, N); % 计算 g 函数: 右子节点 LLR, 需要左子树已译出的比特 u_left = u(start_idx : start_idx + half - 1); y_right = g_llr(y(1:half), y(half+1:2*half), u_left); % 递归译码右子树 u = sc_recursive(y_right, n-1, A, u, start_idx + half, N); end function out = f_llr(a, b) % f 函数: 两个 LLR 的软合并 out = 2 * atanh(tanh(a/2) .* tanh(b/2)); % 数值稳定版本, 避免 tanh 溢出 out(isnan(out)) = 0; end function out = g_llr(a, b, u) % g 函数: 利用已译比特做干扰消除 out = b + (1 - 2*u) .* a; end

逻辑说明：sc_recursive是核心递归体，每次把 LLR 分成左右两半，先算f得到左子节点 LLR 并递归，再用左子树译出的比特算g得到右子节点 LLR 并递归。f_llr用双曲正切公式，数值上容易出 NaN，所以加了一行保护。g_llr里(1-2*u)把 0/1 映射成 +1/-1。参数说明：y是 LLR，正值倾向 0，负值倾向 1；A在递归里其实没直接用，因为冻结比特在叶子节点会被强制置 0——上面代码为了简洁没写这个判断，实际用时要在叶子节点加一句「如果当前下标不在 A 中，强制判 0」。这个坑很隐蔽，不加的话冻结比特会随机判，BER 直接崩。

3. 把 BER 曲线跑出来：仿真框架、信噪比扫描和结果验证

3.1 仿真主循环：帧数、Eb/N0 范围和停止条件

跑通编码和译码只是第一步，真正要验证性能，得写一个仿真主循环，扫 Eb/N0 从 0 到 4 dB，每个点跑足够多的帧，统计误比特率。帧数怎么定？经验值是每个信噪比点至少 1000 帧，或者错误比特数达到 100 就停。N=1024、K=512 时，Matlab 跑 1000 帧大概几秒到十几秒，取决于机器。下面是一个可复用的仿真框架。

function [EbN0_dB, BER] = polar_simulate(N, K, EbN0_dB_range, max_frames) % 极化码 BER 仿真主循环 % N: 码长, K: 信息比特数 % EbN0_dB_range: 信噪比扫描范围(dB) % max_frames: 每个信噪比点最大帧数 % 构造信息位集合, 设计信噪比取扫描范围中点 A = polar_ga_construct(N, K, mean(EbN0_dB_range)); R = K / N; BER = zeros(size(EbN0_dB_range)); for idx = 1:length(EbN0_dB_range) EbN0 = 10^(EbN0_dB_range(idx) / 10); sigma = sqrt(1 / (2 * R * EbN0)); err_bits = 0; total_bits = 0; for frame = 1:max_frames % 生成随机信息比特 u = randi([0 1], K, 1); % 编码 x = polar_encode(u, A, N); % BPSK 调制: 0 -> +1, 1 -> -1 s = 1 - 2 * x; % AWGN 信道 y = s + sigma * randn(N, 1); % LLR 计算 llr = 2 * y / sigma^2; % SC 译码 u_hat_full = polar_sc_decode(llr, A, N); u_hat = u_hat_full(A); % 统计误比特 err_bits = err_bits + sum(u ~= u_hat); total_bits = total_bits + K; % 提前停止: 错误比特超过 100 就跳出 if err_bits >= 100 break; end end BER(idx) = err_bits / total_bits; fprintf('Eb/N0 = %.1f dB, BER = %.4e, frames = %d\n', ... EbN0_dB_range(idx), BER(idx), frame); end end

逻辑说明：外层循环扫信噪比，内层循环跑帧。每帧生成随机信息比特、编码、BPSK 调制、加高斯噪声、算 LLR、SC 译码、统计错误。A只在循环外算一次，因为构造只依赖设计信噪比，不随每帧变。参数说明：max_frames建议 1000~5000，EbN0_dB_range取 0:0.5:4 比较合适；sigma的计算里R是码率，因为 Eb/N0 和符号信噪比差一个码率因子。注意llr = 2*y/sigma^2这个公式假设 BPSK 且噪声方差是sigma^2，如果你改了调制方式，LLR 公式要跟着改。

3.2 结果对不对：三个验证手段和常见异常

跑出 BER 曲线后，怎么判断结果可信？我一般用三个手段交叉验证。第一，看趋势：Eb/N0 增大时 BER 必须单调下降，如果出现平台或者上升，说明译码或 LLR 有问题。第二，和理论值对比：极化码在 SC 译码下，N=1024、K=512 时，Eb/N0=2 dB 附近 BER 大概在 1e-2 量级，3 dB 到 1e-3 左右，如果你跑出来 2 dB 就 1e-5，大概率是信息位集合构造错了或者噪声加错了。第三，换码长验证：N 翻倍时，相同码率下 BER 应该更好，如果 N=2048 比 N=1024 还差，检查A的构造是否随 N 重新计算。

常见异常有三种。一是 BER 恒为 0.5，通常是A在编码和译码端不一致，或者冻结比特没强制置 0。二是 BER 曲线在高信噪比出现「地板」，降不下去，可能是 SC 译码本身的局限——SC 译码不是最优的，码长较短时性能受限，换 SCL 译码会改善。三是仿真跑得特别慢，N=1024 时矩阵乘编码是瓶颈，换成蝶形快速编码能快一个数量级。这些现象背后都有明确原因，不要靠调参数碰运气。

4. 避坑与排查：极化码 Matlab 仿真里最容易翻车的五件事

4.1 现象：BER 一直是 0.5，译码完全无效

原因：最常见的是信息位集合A在编码和译码两端不一致。比如编码用 GA 构造的A，译码却用了默认的全 1 或者别的排序。另一个原因是冻结比特没有在译码叶子节点强制置 0，SC 译码把冻结位当信息位判了。解决：把A存下来，编码和译码传同一个变量；在sc_recursive的n==0分支加判断，如果start_idx不在A里，直接u(start_idx)=0，不看 LLR。

4.2 现象：BER 曲线在高信噪比不下降，出现平台

原因：SC 译码是次优算法，尤其在中短码长下，错误传播会导致性能地板。另外如果A的构造信噪比和仿真信噪比偏差太大，信息位选择不是最优，也会提前出现平台。解决：换 SCL 译码，保留 L 条候选路径，L 取 4 或 8 就能明显改善；或者把A的构造信噪比设到目标工作点附近，比如你关心 2.5 dB 的性能，就用 2.5 dB 构造。

4.3 现象：仿真速度极慢，跑一个点要几分钟

原因：编码用了 N×N 矩阵乘，N=1024 时每次编码 100 万次乘加，几万帧就是几百亿次操作。Matlab 虽然矩阵快，但也扛不住。解决：换成蝶形快速编码，复杂度从 O(N^2) 降到 O(N log N)。核心是把x = mod(G'*u,2)换成逐层异或，每层处理 N/2 对，log2(N) 层搞定。网上有现成的快速编码函数，也可以自己写，十几行的事。

4.4 现象：LLR 计算出 NaN 或 Inf，译码结果乱跳

原因：f_llr里用了tanh(a/2)，当a很大时tanh饱和到 1，atanh(1)就是 Inf。信噪比高的时候 LLR 绝对值大，很容易触发。解决：用数值稳定的 f 函数实现，比如sign(a)*sign(b)*min(|a|,|b|)近似，或者加限幅。上面代码里加了isnan保护，但更稳妥的是直接用 min-sum 近似，精度损失很小，速度还快。

4.5 现象：换码长后性能反而变差

原因：A的构造没有随 N 重新计算，或者 GA 递推的层数写死了。N 变了，log2(N) 变了，递推层数必须跟着变。另外kron构造生成矩阵时，如果 N 不是 2 的幂，log2(N)不是整数，循环次数就错了。解决：在函数入口加断言mod(log2(N),1)==0，确保 N 是 2 的幂；A的构造每次根据当前 N 重新跑，不要复用上一个 N 的结果。

5. 从 SC 到 SCL：一个能明显改善 BER 的进阶改法和验证技巧

SC 译码的局限在于它每一步只保留一条路径，一旦早期判决错了，后面无法挽回。SCL（Successive Cancellation List）译码保留 L 条候选路径，在每一步扩展出 2L 条，然后按路径度量排序保留最好的 L 条，最后选度量最优的路径输出。这个改动在 Matlab 里不算大，但 BER 改善很明显：N=1024、K=512、L=8 时，Eb/N0=2.5 dB 的 BER 能从 SC 的 1e-2 降到 1e-3 量级。

实现 SCL 的关键是路径管理和度量更新。每条路径存一个u_hat向量和一个路径度量 PM，PM 初始为 0，每译一个比特，如果判决比特和 LLR 符号不一致，PM 增加|LLR|。扩展时对每条路径分别试 0 和 1，生成 2L 条，排序后保留 PM 最小的 L 条。Matlab 里用元胞数组存路径最方便，但速度慢；追求速度可以用矩阵批量操作，把 L 条路径的 LLR 拼成矩阵一起算。

function u_hat = polar_scl_decode(y, A, N, L) % SCL 译码, L 为列表大小 % 简化实现: 逐比特译码, 路径用元胞数组管理 n = log2(N); % 初始化 L 条路径 paths = cell(L, 1); for i = 1:L paths{i}.u = zeros(N, 1); paths{i}.pm = 0; paths{i}.llr = y; end for bit_idx = 1:N new_paths = {}; for p = 1:length(paths) % 对每条路径, 尝试判决 0 和 1 for decision = [0, 1] new_path = paths{p}; new_path.u(bit_idx) = decision; % 更新路径度量: 如果判决与 LLR 符号冲突, 累加 |LLR| llr_val = new_path.llr(bit_idx); if (decision == 0 && llr_val < 0) || (decision == 1 && llr_val >= 0) new_path.pm = new_path.pm + abs(llr_val); end new_paths{end+1} = new_path; end end % 按路径度量排序, 保留 L 条 [~, order] = sort(cellfun(@(x) x.pm, new_paths)); paths = new_paths(order(1:min(L, length(new_paths)))); end % 选 PM 最小的路径, 取信息位 [~, best] = min(cellfun(@(x) x.pm, paths)); u_hat = paths{best}.u(A); end

这段代码是简化版，目的是让你看清 SCL 的路径扩展和度量更新逻辑，实际工程里 LLR 的递归计算要复用 SC 的f/g函数，不能像上面这样直接用原始y。参数说明：L一般取 4、8、16，越大性能越好但复杂度线性增长；pm是路径度量，越小越好。验证 SCL 是否生效，最直接的方法是固定A和信噪比，对比 SC 和 SCL 的 BER，如果 SCL 没有改善，检查路径度量更新条件是不是写反了——判决 0 时 LLR 为负才加惩罚，判决 1 时 LLR 为正才加惩罚，这个符号很容易搞反。

我自己的习惯是：每次改完译码器，先跑一个固定噪声种子的单帧，把编码后的x、接收y、译码u_hat全打出来，人工对一遍。这比跑一万帧 BER 更快定位问题。极化码仿真这件事，坑大多不在算法本身，而在下标、符号、一致性这些细节上。希望帮到你。

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