简介这份资源是面向通信工程、信息论学习者与科研人员的极化码基础仿真代码包围绕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、2048K 是信息比特数决定码率 R K/NA 是信息位集合也就是「好信道」的下标集合大小等于 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 调制, 码率 RK/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 dBA返回的是排序后的下标编码时把信息比特放在这些位置上。注意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)在 N1024 时矩阵乘是 1024×1024跑一次几毫秒仿真几万帧就慢了后面可以换成蝶形快速编码。2.3 SC 译码递归函数怎么写、LLR 怎么传SCSuccessive Cancellation译码是极化码的经典译码算法核心是递归的f和g函数。f函数处理两个 LLR 的软信息合并g函数利用已译出的比特做干扰消除。Matlab 里用递归实现最直观但要注意递归深度是 log2(N)N1024 时深度 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(half1: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(half1: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负值倾向 1A在递归里其实没直接用因为冻结比特在叶子节点会被强制置 0——上面代码为了简洁没写这个判断实际用时要在叶子节点加一句「如果当前下标不在 A 中强制判 0」。这个坑很隐蔽不加的话冻结比特会随机判BER 直接崩。3. 把 BER 曲线跑出来仿真框架、信噪比扫描和结果验证3.1 仿真主循环帧数、Eb/N0 范围和停止条件跑通编码和译码只是第一步真正要验证性能得写一个仿真主循环扫 Eb/N0 从 0 到 4 dB每个点跑足够多的帧统计误比特率。帧数怎么定经验值是每个信噪比点至少 1000 帧或者错误比特数达到 100 就停。N1024、K512 时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~5000EbN0_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 译码下N1024、K512 时Eb/N02 dB 附近 BER 大概在 1e-2 量级3 dB 到 1e-3 左右如果你跑出来 2 dB 就 1e-5大概率是信息位集合构造错了或者噪声加错了。第三换码长验证N 翻倍时相同码率下 BER 应该更好如果 N2048 比 N1024 还差检查A的构造是否随 N 重新计算。常见异常有三种。一是 BER 恒为 0.5通常是A在编码和译码端不一致或者冻结比特没强制置 0。二是 BER 曲线在高信噪比出现「地板」降不下去可能是 SC 译码本身的局限——SC 译码不是最优的码长较短时性能受限换 SCL 译码会改善。三是仿真跑得特别慢N1024 时矩阵乘编码是瓶颈换成蝶形快速编码能快一个数量级。这些现象背后都有明确原因不要靠调参数碰运气。4. 避坑与排查极化码 Matlab 仿真里最容易翻车的五件事4.1 现象BER 一直是 0.5译码完全无效原因最常见的是信息位集合A在编码和译码两端不一致。比如编码用 GA 构造的A译码却用了默认的全 1 或者别的排序。另一个原因是冻结比特没有在译码叶子节点强制置 0SC 译码把冻结位当信息位判了。解决把A存下来编码和译码传同一个变量在sc_recursive的n0分支加判断如果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 矩阵乘N1024 时每次编码 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饱和到 1atanh(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 译码的局限在于它每一步只保留一条路径一旦早期判决错了后面无法挽回。SCLSuccessive Cancellation List译码保留 L 条候选路径在每一步扩展出 2L 条然后按路径度量排序保留最好的 L 条最后选度量最优的路径输出。这个改动在 Matlab 里不算大但 BER 改善很明显N1024、K512、L8 时Eb/N02.5 dB 的 BER 能从 SC 的 1e-2 降到 1e-3 量级。实现 SCL 的关键是路径管理和度量更新。每条路径存一个u_hat向量和一个路径度量 PMPM 初始为 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{end1} 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 更快定位问题。极化码仿真这件事坑大多不在算法本身而在下标、符号、一致性这些细节上。希望帮到你。本文还有配套的精品资源点击获取
