番外 B： $\phi$ 反函数、数值近似与 GA 构造实现

番外 A 已经得到减号分支的均值递推：

\mu^- = \phi^{-1}\left(1-\left(1-\phi(\mu_0)\right)\left(1-\phi(\mu_1)\right)\right)

公式本身并不难，真正需要处理的是： $\phi(x)$ 由积分定义，通常没有方便的初等反函数；而一个长度为 $N$ 的 Polar 码需要反复计算它。下面从数值问题开始，最后给出一个可以独立编译的 C++ 构造示例。

先明确输入和输出

GA 构造的输入有三个：母码长 $N=2^n$ 、信息长度 $K$ ，以及设计信噪比 $\rho_{\mathrm{des}}$ 。这里的 $K$ 是进入 Polar 变换的比特数；如果消息向量先附加 CRC，CRC 比特也计入这个 $K$ 。码率为

R=\frac{K}{N}

输出不是一个无序集合，而是一条可靠性排序序列

\boldsymbol{\pi}=(\pi_0,\pi_1,\ldots,\pi_{N-1})

它满足

\mu_{\pi_0}\ge\mu_{\pi_1}\ge\cdots\ge\mu_{\pi_{N-1}}

其中 $\mu_i$ 是母码索引 $i$ 对应的 GA LLR 均值。信息集合最后定义为

\mathcal{A}=\{\pi_0,\ldots,\pi_{K-1}\},\qquad \mathcal{A}^c=[0,N-1]\setminus\mathcal{A}

因此，反函数计算只是中间步骤；最终目标是得到每个索引的均值，并据此排序。

为什么 $\phi^{-1}$ 要数值求解

$\phi(x)$ 的定义为

\phi(x)=1-\frac{1}{\sqrt{4\pi x}}\int_{-\infty}^{+\infty}\tanh\left(\frac{\ell}{2}\right)\exp\left(-\frac{(\ell-x)^2}{4x}\right)\,\mathrm{d}\ell,\qquad x>0

给定 $q$ ，反函数问题是寻找 $x\ge0$ ，使得

\phi(x)=q

积分、双曲正切和高斯密度组合在一起后，没有适合直接写成初等函数的反表达式。若每次递推都进行数值积分，再进行求根，计算量会很大。因此实际 GA 构造通常用一个分段拟合 $\widehat{\phi}(x)$ 代替积分函数，然后求这个近似函数的反函数。

$\phi(x)$ 随 $x$ 增大而下降，所以合法区间内的解至多一个。这一单调性使二分搜索成为可靠的选择：它不需要导数，也不会因为初始猜测不理想而跳到另一个解。

常用的分段近似

下面给出常用的工程拟合形式：

\widehat{\phi}(x)=\begin{cases}\exp\left(-0.4527x^{0.86}+0.0218\right),&0<x<10,\\[4pt]\sqrt{\frac{\pi}{x}}\left(1-\frac{10}{7x}\right)\exp\left(-\frac{x}{4}\right),&x\ge10.\end{cases}

并单独规定

\widehat{\phi}(0)=1

这是拟合式，不是积分定义本身。有限精度下，计算结果可能略微超出 $[0,1]$ ，所以进入递推和反函数前要做区间截断：

\operatorname{clip}(q)=\min(1,\max(0,q))

截断只负责保护浮点计算。它不能把一个已经不准确的近似变成精确值，因此仍要在合适的信噪比和码长下用仿真或密度进化检查构造结果。

第一段为什么可以直接写出近似反函数

当 $0<x<10$ 时，令 $q=\widehat{\phi}(x)$ ，有

q=\exp\left(-0.4527x^{0.86}+0.0218\right)

两边取自然对数：

\ln q=-0.4527x^{0.86}+0.0218

移项后得到

x^{0.86}=\frac{-\ln q+0.0218}{0.4527}

所以第一段的反函数可以近似写成

x\approx\left(\frac{-\ln q+0.0218}{0.4527}\right)^{1/0.86}

这个公式只能用于第一段对应的范围，不能把它外推到 $q$ 接近 0 的情形。实际实现可以先计算第二段起点 $\widehat{\phi}(10)\approx0.039436$ ：当 $q$ 大于这个值时使用显式近似；当 $q$ 更小时，在 $[10,+\infty)$ 上进行二分。

二分搜索的区间不变量

令

F(x)=\widehat{\phi}(x)-q

因为 $\widehat{\phi}(x)$ 单调下降，希望找到满足 $F(x)=0$ 的位置。对较小目标值，先取下界 $x_{\mathrm{lo}}=10$ ，再取一个上界 $x_{\mathrm{hi}}$ ，并不断扩大上界，直到满足

F(x_{\mathrm{lo}})\ge0,\qquad F(x_{\mathrm{hi}})\le0

这就是二分循环的不变量：真实解始终位于闭区间 $[x_{\mathrm{lo}},x_{\mathrm{hi}}]$ 内。

每轮取

x_{\mathrm{mid}}=\frac{x_{\mathrm{lo}}+x_{\mathrm{hi}}}{2}

如果 $\widehat{\phi}(x_{\mathrm{mid}})>q$ ，说明当前均值还太小，应令 $x_{\mathrm{lo}}=x_{\mathrm{mid}}$ ；否则令 $x_{\mathrm{hi}}=x_{\mathrm{mid}}$ 。区间长度每轮减半，经过 $B$ 轮后，误差上界不超过初始区间长度的 $2^{-B}$ 。

边界情况必须提前约定：

$q=1$ 对应 $x=0$ ；
$q=0$ 对应理想化的 $x=+\infty$ ，有限精度程序只能返回预设的最大均值；
$q$ 超出 $[0,1]$ 时先截断，再进行反函数计算。

最大均值是数值饱和值，不是数学上的无穷大。它需要和后续均值相加、LLR 计算以及浮点类型的表示范围一起选择。

一个二分搜索的具体例子

取目标值 $q=0.01$ 。因为

\widehat{\phi}(10)\approx0.039436>0.01,\qquad \widehat{\phi}(20)\approx0.002480<0.01

所以初始区间 $[10,20]$ 已经满足不变量。前几轮如下：

轮次	$x_{\mathrm{lo}}$	$x_{\mathrm{mid}}$	$x_{\mathrm{hi}}$	$\widehat{\phi}(x_{\mathrm{mid}})$	更新
1	10.0000	15.0000	20.0000	0.009738	收紧上界
2	10.0000	12.5000	15.0000	0.019509	收紧下界
3	12.5000	13.7500	15.0000	0.013769	收紧下界
4	13.7500	14.3750	15.0000	0.011576	收紧下界
5	14.3750	14.6875	15.0000	0.010617	收紧下界

继续迭代约 80 轮，得到

\widehat{\phi}^{-1}(0.01)\approx14.903850

这个例子体现了二分搜索的关键：不需要猜测精确答案，只要维持“下界函数值不小于目标、上界函数值不大于目标”的关系，区间就会稳定收缩。

从一个均值递推到全部母码位置

现在说明代码中的分层计算。初始时只有一个均值 $\mu_{\mathrm{ch}}$ 。每进行一层，就把当前序列中的每个均值 $\mu_j$ 变成两个结果：

\mu_{\mathrm{next},2j}=\phi^{-1}\left(1-\left(1-\phi(\mu_j)\right)^2\right),\qquad \mu_{\mathrm{next},2j+1}=2\mu_j

这里 $j$ 是当前层结果序列中的 0-based 位置， $2j$ 和 $2j+1$ 是下一层序列中的位置。它们只是数组位置，不代表额外的通信对象。

以番外 A 的 $N=4$ 、初始均值 2 为例：

\begin{aligned} \text{第 0 层}:&\quad (2),\\ \text{第 1 层}:&\quad (0.823364,4),\\ \text{第 2 层}:&\quad (0.209864,1.646729,2.282073,8). \end{aligned}

第 2 层的顺序是递推顺序。若采用 $\mathbf{G}_N=\mathbf{B}_N\mathbf{F}_2^{\otimes n}$ ，需要用比特反转把递推序列放到母码索引上。对 $N=4$ ，映射后为

(\mu_0,\mu_1,\mu_2,\mu_3)=(0.209864,2.282073,1.646729,8)

再按均值降序排序，就得到

\boldsymbol{\pi}=(3,1,2,0)

这一步要特别区分两个概念：均值向量的下标是母码位置，排序序列的元素是这些位置的编号。不能把“第几个排序结果”误当成“母码位置”。

一段完整的 C++ 示例

下面的示例只负责 GA 构造：输入 $N$ 、 $K$ 和 $\rho_{\mathrm{des}}$ ，输出从可靠到不可靠的索引序列。索引全部从 0 开始； $\rho_{\mathrm{des}}$ 以 dB 输入，函数内部先转换为线性值。这里的 LLR 约定为 $L=\ln(W(Y\mid0)/W(Y\mid1))$ ，发送 0 时正 LLR 表示支持正确方向，均值越大表示子信道越可靠。代码没有包含编码、译码、rate matching 或 AWGN 仿真，这些模块应使用各自的 $\rho_{\mathrm{sim}}$ 和参数约定。

代码中的 layer[j] 表示某一层递推结果序列的第 $j$ 个均值；next[2 j] 是减号结果，next[2 j + 1] 是加号结果。每轮处理后，layer 的长度翻倍。最后的 bitReverse 只负责把递推顺序映射到采用比特反转约定的母码索引。

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>

namespace polar_ga_example {

constexpr double kPi = 3.14159265358979323846;
constexpr double kMaxMean = 1.0e4;

double clip01(double value) {
    return std::min(1.0, std::max(0.0, value));
}

double phiApprox(double x) {
    if (x <= 0.0) {
        return 1.0; // 零均值没有方向信息，phi(0)按定义取1。
    }

    double value = 0.0;
    if (x < 10.0) {
        value = std::exp(-0.4527 * std::pow(x, 0.86) + 0.0218);
    } else {
        value = std::sqrt(kPi / x)
              * (1.0 - 10.0 / (7.0 * x))
              * std::exp(-x / 4.0);
    }
    return clip01(value); // 拟合和舍入误差不能把结果带出[0,1]。
}

double phiInverse(double target) {
    target = clip01(target);
    if (target >= 1.0) {
        return 0.0; // phi(0)=1。
    }
    if (target <= 0.0) {
        return kMaxMean; // phi^{-1}(0)为无穷大，程序用饱和值表示。
    }

    const double phiAtTen = phiApprox(10.0);
    if (target >= phiAtTen) {
        // 目标落在第一段，直接对第一段拟合式取反函数。
        const double numerator = -std::log(target) + 0.0218;
        const double estimate =
            std::pow(numerator / 0.4527, 1.0 / 0.86);
        return std::min(10.0, std::max(0.0, estimate));
    }

    // 二分不变量：phi(lo)>=target，phi(hi)<=target，解始终在[lo,hi]内。
    double lo = 10.0;
    double hi = 20.0;
    while (phiApprox(hi) > target && hi < kMaxMean) {
        hi = std::min(kMaxMean, 2.0 * hi);
    }

    // 80轮足以把双精度区间压得很窄。
    for (int iteration = 0; iteration < 80; ++iteration) {
        const double mid = 0.5 * (lo + hi);
        if (phiApprox(mid) > target) {
            lo = mid; // phi仍偏大，均值需要向右移动。
        } else {
            hi = mid; // phi已不大于目标，均值需要向左收紧。
        }
    }
    return 0.5 * (lo + hi);
}

std::size_t bitReverse(std::size_t index, std::size_t width) {
    std::size_t reversed = 0;
    for (std::size_t bit = 0; bit < width; ++bit) {
        reversed = (reversed << 1) | (index & 1U);
        index >>= 1;
    }
    return reversed;
}

std::vector<double> buildMeans(std::size_t N, double rootMean) {
    std::size_t depth = 0;
    for (std::size_t size = N; size > 1; size >>= 1) {
        ++depth;
    }

    std::vector<double> layer(1, rootMean);
    for (std::size_t level = 0; level < depth; ++level) {
        std::vector<double> next(2 * layer.size(), 0.0);
        for (std::size_t j = 0; j < layer.size(); ++j) {
            const double p = phiApprox(layer[j]);
            const double minusTarget = 1.0 - (1.0 - p) * (1.0 - p);

            // 偶数位置保存减号结果，奇数位置保存加号结果。
            next[2 * j] = phiInverse(minusTarget);
            next[2 * j + 1] =
                std::min(kMaxMean, 2.0 * layer[j]);
        }
        layer.swap(next);
    }

    // layer是递推顺序；这里映射为母码索引i的均值mu[i]。
    std::vector<double> means(N, 0.0);
    for (std::size_t i = 0; i < N; ++i) {
        means[i] = layer[bitReverse(i, depth)];
    }
    return means;
}

std::vector<std::size_t> gaReliabilityOrder(
    std::size_t N, std::size_t K, double rhoDesDb) {
    if (N == 0 || (N & (N - 1)) != 0 || K == 0 || K > N) {
        throw std::invalid_argument(
            "N must be a power of two and 0 < K <= N");
    }

    const double rate = static_cast<double>(K) / N;
    const double rhoDesLin = std::pow(10.0, rhoDesDb / 10.0);
    // Es=1且rhoDesLin=Eb/N0时，BI-AWGN初始均值为4R(Eb/N0)。
    const double rootMean = 4.0 * rate * rhoDesLin;
    const std::vector<double> means = buildMeans(N, rootMean);

    std::vector<std::size_t> order(N);
    for (std::size_t i = 0; i < N; ++i) {
        order[i] = i;
    }

    // 先按均值降序；均值相同时按母码索引升序，保证结果确定。
    std::sort(order.begin(), order.end(),
        [&means](std::size_t a, std::size_t b) {
            if (means[a] != means[b]) {
                return means[a] > means[b];
            }
            return a < b;
        });

    return order; // order[0]最可靠，order[N-1]最不可靠。
}

} // namespace polar_ga_example

代码中的几个边界值得单独记住。phiApprox(0) 对应完全没有方向信息；phiInverse(1) 必须返回 0；目标值接近 0 时，均值会很大，程序只能使用 kMaxMean 表示有限精度下的饱和状态。二分搜索每轮都维持同一个区间不变量，排序则明确采用“均值降序、索引升序”的规则。

如果取 $N=4$ 、 $K=2$ 、 $\rho_{\mathrm{des}}=0$ dB，则 $R=1/2$ ，初始均值为

\mu_{\mathrm{ch}}=4\times\frac12\times10^0=2

代码得到的母码均值约为

(\mu_0,\mu_1,\mu_2,\mu_3)=(0.209864,2.282073,1.646729,8)

因此排序结果为 $\boldsymbol{\pi}=(3,1,2,0)$ ，与番外 A 的手算完全一致。若取前 $K=2$ 个位置，信息集合就是 $\mathcal{A}=\{3,1\}$ 。

如何验收构造结果

一个可用的构造实现至少应检查以下关系：

$N$ 是 2 的整数次幂，且 $0<K\le N$ 。
输出序列长度为 $N$ ，每个索引都属于 $[0,N-1]$ 。
输出序列没有重复索引。
对所有相邻位置，排序规则满足 $\mu_{\pi_j}\ge\mu_{\pi_{j+1}}$ 。
改变 $\rho_{\mathrm{des}}$ 后重新构造，不能把结果继续误认为与信道无关的固定顺序。
用 $\rho_{\mathrm{sim}}$ 做 AWGN 仿真时，记录实际采用的码率、噪声方差和统计对象；不要把构造信噪比直接当成仿真信噪比。

如果只想检查反函数，可以选取若干 $x$ ，先计算 $q=\widehat{\phi}(x)$ ，再检查 $\widehat{\phi}(\widehat{\phi}^{-1}(q))$ 是否接近 $q$ 。当 $q$ 很小时，优先检查上界扩大和均值饱和逻辑；当排序出现重复或越界，优先检查比特反转映射和排序容器的初始化。

小结

$\phi^{-1}$ 的实现依赖三个事实： $\phi(x)$ 单调下降，第一段拟合式可以显式求反，以及更小目标值可以用带区间不变量的二分搜索处理。完成每个位置的均值计算后，还必须经过比特反转映射，再按均值降序、索引升序生成可靠性排序 $\boldsymbol{\pi}$ 。

GA 的优点是把分布进化压缩成均值递推，计算简单、适合构造；它的局限是减号分支只保留了高斯近似下的一个统计描述。需要更高精度时，应回到密度进化或通过固定参数的仿真验证构造结果。

参考文献

E. Arıkan, “Channel polarization: A method for constructing capacity-achieving codes for symmetric binary-input memoryless channels,” IEEE Transactions on Information Theory, vol. 55, no. 7, pp. 3051–3073, Jul. 2009, doi: 10.1109/TIT.2009.2021379.
D. Trifonov, “Efficient design and decoding of polar codes,” IEEE Transactions on Communications, vol. 60, no. 11, pp. 3221–3227, Nov. 2012, doi: 10.1109/TCOMM.2012.090512.110070.

番外 B：\phi 反函数、数值近似与 GA 构造实现