[信息论]信道容量迭代算法程序设计(基于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 1 p 2 ⋯ p n } \pmb{{p^{(0)}}}\left\{

\right\} p(0)p(0)⎩ ⎨ ⎧​p1​p2​⋯pn​​⎭ ⎬ ⎫​(初始化为均匀分布)，信道转移概率矩阵 P = { p 11 p 12 ⋯ p 1 m p 21 p 22 ⋯ p 2 m ⋮ ⋮ ⋱ ⋮ p n 1 p n 2 ⋯ p n m } \pmb{P}= \left\{

p11p12⋯p1mp21p22⋯p2m⋮⋮⋱⋮pn1pn2⋯pnm" role="presentation" style="position: relative;">p11p21⋮pn1p12p22⋮pn2⋯⋯⋱⋯p1mp2m⋮pnm$$\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\} PP=⎩ ⎨ ⎧​p11​p21​⋮pn1​​p12​p22​⋮pn2​​⋯⋯⋱⋯​p1m​p2m​⋮pnm​​⎭ ⎬ ⎫​,信道容量初值 $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
#include
#include
/* 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*******************************/
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162

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);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

测试数据与结果展示

测试数据

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

输入数据

P = { 1 0 1 0 0.5 0.5 0 1 0 1 } \pmb{P}= \left\{

\right\} PP=⎩ ⎨ ⎧​110.500​000.511​⎭ ⎬ ⎫​,$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)