编解码方案性能分析工具：外信息传递图（EXIT chart）及LDPC-EXIT代码参考

2004《Extrinsic Information Transfer Functions: Model and Erasure Channel Properties》学习笔记

  Extrinsic information transfer (EXIT) charts，在通信中用于预测迭代解码器的收敛行为，适用于诸如turbo码、LDPC（low-density parity-check ）码、RA（repeat–accumulate ）码的迭代译码分析。
  论文提出了一种解码模型并定义外信息转移函数，证明了其在BEC信道下的三条相关性质：

1. 条件熵和外信息转移函数曲线所围面积的关联→极限容量码的设计简化为曲线拟合问题
2. 外信息转移函数与码的Helleseth–Kløve–Levenshtein information function关联
3. Split information functions、split support weights和EXIT fuction的关联→线性码的EXIT函数与其对偶码的EXIT函数关联

辅以其他参考文献：
[1] $GF(q)$ 域上非规则 LDPC 码 EXIT 图分析方法研究
[2] EXIT 图分析多进制 LDPC 码

1. 密度演进和外信息转移图

均用于预测迭代译码的收敛性。

  DE跟踪概率密度函数，由仅适用于ECs（erasure channels）的外对数似然函数（L-value）的概率分布演化而来，从而适用于其他实际信道，用于计算码字的收敛门限。
  EXIT chart跟踪每结点平均互信息，依赖互信息这一统计信息的可靠性和在任何信道中的广泛适用性，对各种实际信道设计码字和迭代处理器意义重大。EXIT主要是因为DS跟踪概率密度函数计算量太大，所以跟踪互信息简化计算量。

收敛门限，某一码型迭代译码性能具有门限现象，即存在一个噪声方差门限（最大的信道参数），当码长区域无穷，信道噪声参数低于该门限时，译码比特错误概率能达到无限小。该参数值即为收敛门限，反应最低的信道质量要求。
${\sigma }^2=\frac{1}{2R\times E_b/N_0}$，其中R表示码率

EXIT图两个必要假设：

1. 从一个译码器到另一个译码器的外信息时高斯随机变量（统计特性）
2. 译码器的外信息相互独立

EXIT图即一个译码器的先验信息来自另外一个译码器的外信息。VND输出的外信息和码字之间的互信息为$I_{E1}$，其能用$I_{ch}$和$I_{a1}$表示。对于CND，其输出$I_{E2}$，其能用$I_{a2}$表示。EXIT图即在同一图中画出两个$I_E-I_A$曲线。

2. 解码模型（定义先验信息a、c和外信息e）

该部分定义解码模型并通过a posteriori probability (APP) decoding得到外信息的定义：思路是先模型，后定义向量a和c，再从已知向量中推导出d，并根据d和a的关系定义e。

大写公式字母-随机变量，小写公式字母-实际测量：

• $\underline{u}$为等概取值{0, 1}的信息bit流，长度为K。
• $\underline{x}$为长度为N的编码符号流。
• $\underline{y}$为$\underline{x}$经过communication channel后的含噪声信号。称为SC。
• $\underline{w}$为$\underline{v}$经过外部通道（先验通道）得到的先验信息信号。对于迭代解码器，$\underline{v}$并非来自Extrinsic channel而是迭代解码器中的一个子解码器，但用途还是先验信息。称为PC。
• $\underline{v}$，假设对$\underline{y}$执行迭代解码，解码器可分解为两个子解码器，内部解码器接收关于内部编码器输入位的外部信息，所以设置$\underline{v}=\underline{u}$，外部解码器接受外部编码器输出位的外部信息，所以设置$\underline{v}=\underline{x}$。

  外部通道无需任何现实设备辅助（Encoder 2 并非真实存在），其引入是为了更好的理解和进行性能分析。

  如果解码器的先验值来自BEC，那么EXIT函数下的面积是1减去条件熵，外部信道非BEC阻碍迭代解码的适用性。

  设$\underline{v},\underline{w},\underline{a}$, and $\underline{e}$皆为长度为m的向量，解码器使用 $\underline{y}$ 和 $\underline{w}$ 来计算$\underline{v}$的两个估计量：先验L-values（log-likelihood ratios）：$\underline{d}$和外部L-values：$\underline{e}$

对于无记忆信道（中文文献管这些叫比特似然值）：
  编码符号$w_i$附带关于随机变量$V_i$的先验信息$a_i=\log \frac{P\left(w_i\mid V_i=0\right)}{P\left(w_i\mid V_i=1\right)}$
  编码符号$y_i$携带关于随机变量$X_i$的先验信息$c_i=\log \frac{P\left(y_i\mid X_i=0\right)}{P\left(y_i\mid X_i=1\right)}$

posteriori probability decoder计算先验L-value：$$d_i=\log \frac{\mathrm{Pr}\left(V_i=0\mid \underline{y},\underline{w}\right)}{\mathrm{Pr}\left(V_i=1\mid \underline{y},\underline{w}\right)}$$

当采用MAP解码$\underline{u}$中bit流时，外部通道即完全噪声通道。
TODO：这里AAP decoingr和MAP decoding，区别在于是否最大化后验概率？

MAP解码下，内部解码器接收关于内部编码器输入位的外部信息$\underline{v}=\underline{u}$，所以设置。展开$d_i$的分子分母：$$\begin{array}{rl}\mathrm{Pr}(V_i& =0\mid \underline{y},\underline{w})\\ & =\sum _{\underline{u}:v_i(\underline{u})=0}P(\underline{u}\mid \underline{y},\underline{w})\\ =& \sum _{\underline{u}:v_i(\underline{u})=0}\frac{P(\underline{u})P(\underline{w}\mid \underline{u})P(\underline{y}\mid \underline{u},\underline{w})}{P(\underline{y},\underline{w})}\\ =& \sum _{\underline{u}:v_i(\underline{u})=0}\frac{P(\underline{u})P(\underline{w}\mid \underline{v}(\underline{u}))P(\underline{y}\mid \underline{x}(\underline{u}))}{P(\underline{y},\underline{w})}\\ =& \frac{P\left(w_i\mid V_i=0\right)}{P(\underline{y},\underline{w})}\\ & \cdot \sum _{\underline{u}:v_i(\underline{u})=0}P(\underline{u})P\left(\underline{w_{[i]}}\mid {\underline{v}}_{[i]}(\underline{u})\right)P(\underline{y}\mid \underline{x}(\underline{u}))\end{array}$$其中$\underline{v}(\underline{u})$ and $\underline{x}(\underline{u})$表示对应于$\underline{u}$的向量。$$\begin{array}{rl}{e}_i& =\log \frac{\mathrm{Pr}\left(V_i=0\mid \underline{y},{\underline{w}}_{[i]}\right)}{\mathrm{Pr}\left(V_i=1\mid \underline{y},{\underline{w}}_{[i]}\right)}\\ & =\log \frac{{\sum }_{\underline{u}:v_i(\underline{u})=0}P\left({\underline{w}}_{[i]}\mid {\underline{v}}_{[i]}(\underline{u})\right)P(\underline{y}\mid \underline{x}(\underline{u}))}{{\sum }_{\underline{u}:v_i(\underline{u})=1}P\left({\underline{w}}_{[i]}\mid {\underline{v}}_{[i]}(\underline{u})\right)P(\underline{y}\mid \underline{x}(\underline{u}))}.\end{array}$$可以推出：$d_i=a_i+e_i$。称$e_i$为$v_i$的外部$L-value$

3. 外信息转移函数

  迭代解码器有两个或多个子解码器，解码器之间交换外部对视似然函数值（或交换外部概率）。随后的分析使用互信息。

  来自一个解码器的$e_i$通过交织器作为先验值反馈给另外一个解码器，定义：

1. $I_A=\frac{1}{m}{\sum }_{i=1}^{m}I\left(V_i;A_i\right)$，进入解码器的平均先验信息
2. $I_E=\frac{1}{m}{\sum }_{i=1}^{m}I\left(V_i;E_i\right)$，译码器输出的平均外部信息

外信息转移函数图即$I_E$关于$I_A$的曲线，它们之间可以通过码字的参数和信道的参数联系起来。

从$e_i$表达式中可知其是关于$\underline{Y}$和${\underline{W}}_{[i]}$的函数，${\underline{W}}_{[i]}$ and ${\underline{A}}_{[i]}$可根据$a_i$定义式相互推导，所以$I\left(V_i;E_i\right)\le I\left(V_i;\underline{Y}{\underline{W}}_{[i]}\right)=I\left(V_i;\underline{Y}{\underline{A}}_{[i]}\right)=I\left(V_i;\underline{Y},{\underline{A}}_{[i]}\right)$，在具有外部消息传递的后验解码器中，该不等式取等$I\left(V_i;E_i\right)=I\left(V_i;\underline{Y}{\underline{A}}_{[i]}\right)$。

如果所有$V_i$分布相同，且外部信道无记忆且时不变，那么$I_A=I\left(V_1;A_1\right)\in [-1,1]$。

如果采用后验解码器，$I_E=\frac{1}{m}{\sum }_{i=1}^{m}I\left(V_i;\underline{Y}\underline{A_{[i]}}\right)$

7. LDPC码的EXIT分析

第一个下标表示先验信息a或者外信息e，第二个字母表示解码器v或者c

$I_{\mathbf{e}\mathbf{v}}=f^{\prime }\left(I_{av},I_{ch}\right)=f\left(I_{av},E_b/N_0\right)$表示VND输出的互信息时关于输入互信息和信道信息的函数f。
$I_{ec}=g\left(I_{ac}\right)$表示CND输出的互信息仅时输入互信息的函数。
  EXIT图即$f和g^{-1}$在同一图中，描述平均互信息的传递过程, 在不需要仿真误比特率曲线的前提下可直观地判断译码的收敛性, 得到特定码型的收敛门限值: 如果在 $E_b/N_0$ 一定的情况下, 上述两条曲线除了在（1，1）外均不相交, 表明该信噪比情况下译码器能够收敛; 减小 $E_b/N_0$, 两条曲线之间的距离会逐渐变小, 直到两曲线在没有其它交点的情况下相距最近, 此时得到的信噪比对应的噪声方差就是所求的收敛门限值。

主要难点在于公式的推导，即3中俩公式在特定码字中的化简过程。推导过程公式太多不好打，略去。
参考文献：LDPC码的EXIT图
规则码的EXIT图，参考自github[LDPC_tools-main]：

function out = J_sigma(sigma)
% this function is one of the steps in the EXIT calculation
% Ref. [1] Design of Low-Density Parity-Check Codes for Modulation and
% Detection -- Stephan ten Brink et al. Appendix


sigma_star = 1.6363;
aj1 = -0.0421061;
bj1 = 0.209252;
cj1 = -0.00640081;
aj2 = 0.00181491;
bj2 = -0.142675;
cj2 = -0.0822054;
dj2 = 0.0549608;

if sigma >= 0 && sigma <= sigma_star
    
    out = sum([aj1 bj1 cj1].* (sigma.^([3 2 1])));
    
elseif sigma > sigma_star
    
    out = 1 - exp(sum([aj2 bj2 cj2 dj2].* (sigma.^([3 2 1 0]))));
    
else
    
    out = 1;
    
end

function out = J_sigma_inv(ei)
% this function is one of the steps in the EXIT calculation
% Ref. [1] Design of Low-Density Parity-Check Codes for Modulation and
% Detection -- Stephan ten Brink et al. Appendix

ei_star = 0.3646;
as1 = 1.09542;
bs1 = 0.214217;
cs1 = 2.33727;
as2 = 0.706692;
bs2 = 0.386013;
cs2 = -1.75017;

if ei >= 0 && ei <= ei_star
   
    out = sum([as1 bs1 cs1].* (ei.^([2 1 0.5])));
    
elseif ei > ei_star && ei < 1
    
    out = - as2 * log(bs2*(1 - ei)) - cs2 * ei;
    
else
    
    out = 10;
    fprintf('Numerical error in the inverse J_sigma function\n');
    
end

function out = Iev_Iav(I_av, sigma_ch, dv)
% this function calculate the EXIT curve for variable nodes with degree dv
% Ref. [1] Channel Codes classical and modern -- William E. Ryan and Shu Lin 
% (9.42)

I_ev = zeros(1,length(I_av));
for ii = 1:length(I_av)
    
    tmp = J_sigma_inv(I_av(ii));
    J_arg = sqrt((dv - 1)*tmp^2+sigma_ch^2);
    I_ev(ii) = J_sigma(J_arg);
    
end

out = I_ev;

function out = Iec_Iac(I_ac,  dc)
% this function calculate the EXIT curve for check nodes with degree dc
% Ref. [1] Channel Codes classical and modern -- William E. Ryan and Shu Lin 
% (9.43)

I_ec = zeros(1,length(I_ac));

for ii = 1:length(I_ac)
    
    tmp = J_sigma_inv(1 - I_ac(ii));
    J_arg = sqrt((dc - 1)*tmp^2);
    I_ec(ii) = 1 - J_sigma(J_arg);
    
end

out = I_ec;

clear
close all
% this script calculate the EXIT chart for regular LDPC ensemble 
% Ref. [1] Channel Codes classical and modern -- William E. Ryan and Shu Lin 
% (9.42) (9.43)


dv = 3;
dc = 6;
code_rate = dv/dc;
ll = 0;
Ps = 1;
EbN0_dB = 1.1;
EbN0 = 10^(EbN0_dB/10);
sigma = sqrt(Ps/EbN0);
sigma_ch = sqrt(8 * code_rate * EbN0);


I_av = 0:0.05:1;
I_ev = Iev_Iav(I_av, sigma_ch, dv);


plot(I_av, I_ev)

I_ac = 0:0.05:1;
I_ec = Iec_Iac(I_ac, dc);

hold on
plot(I_ec, I_ac)