[信息论]信道容量迭代算法程序设计(基于C++&Matlab实现)

算法分析

迭代法计算信道容量 $C$ 的步骤如下：

首先，记 $p(y_j|x_i)=p_{ij}$ , $p(x_i)=p_i$ , $p(x_i|y_j)={\varphi }_{ji}$ ,其中 $i=1,2,\cdots ,n;j=1,2,\cdots ,m$.算法中信息量的单位是奈特。

1.初始化参数

初始条件已知信源分布 $p^{(0)}{p}^{(0)}\left\{\begin{array}{c}{p}_1\\ p_2\\ \cdots \\ p_n\end{array}\right\}$(初始化为均匀分布)，信道转移概率矩阵 $PP=\left\{\begin{array}{cccc}{p}_{11}& p_{12}& \cdots & p_{1m}\\ p_{21}& p_{22}& \cdots & p_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ p_{n1}& p_{n2}& \cdots & p_{nm}\end{array}\right\}$,信道容量初值 $C^{(0)}{C}^{(0)}=-\infty$,信道容量相对误差门限 $\delta >0$, 状态参数 $k$ ,初值为 $0$.

2.迭代后验概率矩阵

$${\varphi }_{ji}^{(k)}=\frac{p_{ij}{p}_i^{(k)}}{{\sum }_i{p}_{ij}{p}_i^{(k)}}$$

3.迭代信源分布

由于算法中采用奈特作为信息量单位，同时已知下式
$$a^b=e^{blna}$$
因此迭代信源分布状态转移方程可以化简为：
$$p_i^{(k+1)}=\frac{e^{{\sum }_j{p}_{ij}ln{\varphi }_{ji}^{(k)}}}{{\sum }_i{e}^{{\sum }_j{p}_{ij}ln{\varphi }_{ji}^{(k)}}}=\frac{{\prod }_j{({\varphi }_{ji}^{(k)})}^{p_{ij}}}{{\sum }_i{\prod }_j{({\varphi }_{ji}^{(k)})}^{p_{ij}}},i=1,\cdots ,r$$
这样进行处理的好处是，后向概率中有些元素的值为 $0$ ,对于大部分编程语言来说，$ln0$ 的运算一般是输出 $nan$ 运算符，无法得到正确计算结果。

4. 迭代信道容量

由 $(3)$ 中的分析可知，信道容量迭代方程可以化简为
$$C^{(k+1)}=ln(\sum _i{e}^{{\sum }_j{p}_{ij}ln{\varphi }_{ji}^{(k)}})=ln(\sum _i\prod _j{({\varphi }_{ji}^{(k)})}^{p_{ij}})$$

5.判定信道容量是否达到 $\delta$ 精度范围要求

若 $\frac{|C^{(k+1)}-C^{(k)}|}{C^{(k+1)}}⩽\delta$ ,则停止迭代输出结果，否则继续迭代。

C++ Codes

/**
  ******************************************************************************
  * @file    Project/Template/main.cpp
  * @author  水声工程学院 LZH 2019051312
  * @version V1.0.0
  * @date    21-March-2022
  * @brief   Main program body
  ******************************************************************************
  */

/* Includes ------------------------------------------------------------------*/
#include<iostream>
#include<iomanip>
#include<cmath>
/* Private typedef -----------------------------------------------------------*/
/* Private define ------------------------------------------------------------*/
#define ll long long
#define ld long double
#define pb push_back
#define eps 1e-8  //Change this value to upgrade precision
using namespace std;
/* Private macro -------------------------------------------------------------*/
/* Private variables ---------------------------------------------------------*/
ld P_ij[105][105]={0};  //Transition probability matrix
ld Phi_ij[105][105]={0}; //Posterior probability matrix
ld p[105]; //Prior probability matrix
ld C0=-114514,C1=0; //C0:C(k) C1:C(k+1)
int n,m; //n rows and m columns
/* Private function prototypes -----------------------------------------------*/
ld C_update(int n,int m,ld (&phi_ij)[105][105],ld (&p_ij)[105][105]);
void p_update(int n,int m,ld (&phi_ij)[105][105],ld (&p_ij)[105][105]);
void Phi_update(int n,int m,ld (&phi_ij)[105][105],ld (&p_ij)[105][105]);
bool check(ld a,ld b);
/* Private functions -----------------------------------------------*/
/**
  * @brief  Get abs value of the input parameter.
  * @param  A long double format parameter.
  * @retval The abs value of the input parameter.
  */
ld abss(ld a){
   if (a<0)
	  return -a;
   return a;
}
/**
  * @brief  Main program.
  * @param  None
  * @retval None
  */
int main(){
	ios::sync_with_stdio(0); //accelerate cin and cout.
	int step=0; //Steps of current iteration.
	/*Get all parameters*/
	cin>>n>>m;
	for (int i=1; i<=n; i++){
		p[i]=1.0/(ld)n; //Initialized array p into evenly distributed state.
		for (int j=1; j<=m; j++)
			cin>>P_ij[i][j];
	}

	/*Iteration procedure*/
	while (!check(C1,C0) || !step){
		C0=C1;
		Phi_update(n,m,Phi_ij,P_ij); //Update phi_{k} first.
		p_update(n,m,Phi_ij,P_ij); //Update p_{k+1}
		C1=C_update(n,m,Phi_ij,P_ij);//Update C_{k+1}
		step++; //Transfer to next state
	}
	
	/*Print out all answers*/
	cout.setf(ios::fixed); 
	for (int i=1; i<=n; i++) //Output Prior probability matrix
		cout<<setprecision(15)<<p[i]<<' ';
	cout<<'\n';
	cout<<C1<<endl; //Output C(k+1)
	return 0;
}


/**
  * @brief  Iterate the value of C(k+1)
  * @param  n: total number of rows. Passed by value.
  * @param  m: total number of columns. Passed by value.
  * @param  phi_ij: The Posterior probability matrix. Passed by reference.
  * @param  p_ij: The Prior probability matrix. Passed by reference.
  * @retval C(k+1)
  */
ld C_update(int n,int m,ld (&phi_ij)[105][105],ld (&p_ij)[105][105]){ 
    ld sum=0,sum_tmp=0,ans;
	for (int i=1; i<=n; i++){
		sum_tmp=1;
		for (int j=1; j<=m; j++)
			sum_tmp*=pow(phi_ij[i][j],p_ij[i][j]);
		sum+=sum_tmp;
	}
	ans=log(sum);
	return ans;
}

/**
  * @brief  Iterate the sequence of p.
  * @param  n: total number of rows. Passed by value.
  * @param  m: total number of columns. Passed by value.
  * @param  phi_ij: The Posterior probability matrix. Passed by reference.
  * @param  p_ij: The Prior probability matrix. Passed by reference.
  * @retval None.
  */
void p_update(int n,int m,ld (&phi_ij)[105][105],ld (&p_ij)[105][105]){ 
	ld sum=0,sum_tmp=0;
	for (int i=1; i<=n; i++){
		sum_tmp=1;
		for (int j=1; j<=m; j++)
			sum_tmp*=pow(phi_ij[i][j],p_ij[i][j]);
		sum+=sum_tmp;
	}
	
	for (int i=1; i<=n; i++){
		sum_tmp=1;
		for (int j=1; j<=m; j++)
			sum_tmp*=pow(phi_ij[i][j],p_ij[i][j]);
		p[i]=sum_tmp/sum;
	}
	return;
}

/**
  * @brief  Iterate the Posterior probability matrix.
  * @param  n: total number of rows. Passed by value.
  * @param  m: total number of columns. Passed by value.
  * @param  phi_ij: The Posterior probability matrix. Passed by reference.
  * @param  p_ij: The Prior probability matrix. Passed by reference.
  * @retval None.
  */
void Phi_update(int n,int m,ld (&phi_ij)[105][105],ld (&p_ij)[105][105]){ 
	for (int i=1; i<=n; i++)
		for (int j=1; j<=m; j++){
			ld sum=0;
			for (int k=1; k<=n; k++)
				sum+=p_ij[k][j]*p[k];
			phi_ij[i][j]=p_ij[i][j]*p[i]/sum;
		}
	return;
}

/**
  * @brief  Check whether or not C(k+1) satisfies eps range.
  * @param  a: C(k+1). Passed by value.
  * @param  b: C(k). Passed by value.
  * @retval True or False.
  */
bool check(ld a,ld b){
	ld tot;
	tot=abss(a-b)/a;
	if (tot<=eps)
		return true;
	return false;
}
/**
  * @}
  */

/************************END OF FILE*******************************/

Matlab Codes

eps=1e-8;
owari=0;
C0=-114514;
C1=0;
n=input('Input row num=');
m=input('Input column num=');
Phi_ij=zeros(n,m);
for ii=1:n
    p(ii)=1/n;
    for jj=1:m
        P_ij(ii,jj)=input(''); %Input the Prior probability matrix
    end
end

while (~owari)
       C0=C1;
       %Iterate Matrix Phi
       for ii=1:n
           for jj=1:m
               sum=0;
               for kk=1:n
                   sum=sum+P_ij(kk,jj)*p(kk);
               end
               Phi_ij(ii,jj)=P_ij(ii,jj)*p(ii)/sum;
           end
       end
       %Iterate the sequence of p
       sum=0; sum_tmp=0;
       for ii=1:n
           sum_tmp=1;
           for jj=1:m
               sum_tmp=sum_tmp*(Phi_ij(ii,jj)^P_ij(ii,jj));
           end
           sum=sum+sum_tmp;
       end
       
       for ii=1:n
           sum_tmp=1;
           for jj=1:m
               sum_tmp=sum_tmp*(Phi_ij(ii,jj)^P_ij(ii,jj));
           end
           p(ii)=sum_tmp/sum;
       end
       %Iterate the value of C(k+1)
       sum=0; sum_tmp=0;
       for ii=1:n
           sum_tmp=1;
           for jj=1:m
               sum_tmp=sum_tmp*(Phi_ij(ii,jj)^P_ij(ii,jj));
           end
           sum=sum+sum_tmp;
       end
       C1=log(sum);
       %check whether or not C(k+1) satisfies eps range.
       tot=abs(C1-C0)/C1;
       if (tot<=eps)
           owari=1;
       end
end

%Output the results
fprintf('p(k+1)=[');
for ii=1:n
    fprintf('%20.15f ',p(ii));
end
str=[']'];
disp(str);
fprintf('C(k+1)=%20.15f\n',C1);

测试数据与结果展示

测试数据

测试数据使用教材 $P61例4.7$ 测试数据。

输入数据

$PP=\left\{\begin{array}{cc}1& 0\\ 1& 0\\ 0.5& 0.5\\ 0& 1\\ 0& 1\end{array}\right\}$,$eps=10^{-8}$, $p^{(0)}{p}^{(0)}=(0.2,0.2,0.2,0.2,0.2)$,$C^{(0)}{C}^{(0)}=-\infty$

输出数据

$p^{(k+1)}{p}^{(k+1)}=(0.25,0.25,0,0.25,0.25)$,$C^{(k+1)}{C}^{(k+1)}=1bit=0.6931471806nat$

测试结果(C++)

测试结果(Matlab)