信道容量的迭代算法

实验目的

熟悉信道容量的迭代算法；
学习如何将复杂的公式转化为程序；

实验要求

已知：信源符号个数r，新宿符号份额数s、信道转移概率矩阵p；
输入：任意的一个信道转移概率矩阵。r、s、p在运行时从键盘输入。
输出：最佳信源分布P‘，信道容量C。

实验内容

信道容量的含义

一个信道，若给定输入分布Q(x) ，和转移概率分布P(y|x) ，则输出概率分布为
Θ(y)=∑xQ(x)P(y|x)<script type="math/tex" id="MathJax-Element-1">\Theta(y)=\sum_x{Q(x)P(y|x)}</script>
对DMC信道有Θ(y)=∑xQ(x)P(y|x)<script type="math/tex" id="MathJax-Element-2">\Theta(y)=\sum_x{Q(x)P(y|x)}</script>
当研究联合 空间XY时，则信道输出和输入之间的平均互信息量为

<script type="math/tex; mode=display" id="MathJax-Element-3">I(X;Y)=\sum_{k=0}^{K-1}\sum_{j=0}^{J-1}{Q_kP(j|k)log\frac{P(j|k)}{\sum_{k=0}^{K-1}Q_iP(j|i)}}</script>
<script type="math/tex; mode=display" id="MathJax-Element-4">=H(Y)-H(Y|X)</script>
<script type="math/tex; mode=display" id="MathJax-Element-5">=H(X)-H(X|Y)</script>
从上述公式中，当信道输入为X，接收端得到的平均信息量 I(X;Y)是信道输入概率和信道转移概率的函数。一般干扰越大，传输的信息量就越小。现在研究信道哦给定，即 {P(y|x)}给定的情况下，如何选择输入分布 Qk<script type="math/tex" id="MathJax-Element-6">{Q_k}</script>使传输的信息量 I(X;Y)达到最大。
信道容量的定义：
<script type="math/tex; mode=display" id="MathJax-Element-7">C=max{I(X;Y)}</script>
即C为改变输入分布时，使每个符号所能含有的平均互信息量的最大值。相应的输入分布称为最佳输入分布。

信道容量的迭代算法

算法流程：
1. 取出时分布 r(0)<script type="math/tex" id="MathJax-Element-8">r^{(0)}</script>，置k=0。
2. 由公式

<script type="math/tex; mode=display" id="MathJax-Element-9">Q^{(k)}(x|y)=\frac{P(y|x)r^{(k)}(x)} {\sum_xP(y|x^{'})r^{(k)}(x^{'})}</script>计算Q(k)(x|y)<script type="math/tex" id="MathJax-Element-10">Q^{(k)}(x|y)</script>
3. 由公式 <script type="math/tex; mode=display" id="MathJax-Element-11">r^{k+1}(x) = \frac{\prod_y{Q^{(k)}(x|y)^{P(y|x)}}}{\sum_x{'}\prod_y Q^{(k)(x^{'}|y)^{P(y|x)}}}</script>计算r(k+1)<script type="math/tex" id="MathJax-Element-12">r^{(k+1)}</script>
4. 由公式
<script type="math/tex; mode=display" id="MathJax-Element-13">C^{(k+1)} = I(r^{(k+1)},Q^{(k)})</script>
5. 若|C(k+1)−C(k)|≤δ<script type="math/tex" id="MathJax-Element-14">|C{(k+1)} - C^{(k)}| \leq \delta</script>，进入第7步。
6. 令k=k+1<script type="math/tex" id="MathJax-Element-15">k=k+1</script>, 转向第二步
7. 程序结束输出最佳输入分布r(k+1)<script type="math/tex" id="MathJax-Element-16">r^{(k+1)}</script>及C(k+1)<script type="math/tex" id="MathJax-Element-17">C^{(k+1)}</script>

算法实现

根据算法的基本原理和迭代的流程进行程序编写。
程序代码如下

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
void printMessage();
int main()
{
    int r, s, i, j, k = 0;
    double p[20]; //存放输入信源概率矩阵 
    double z[20];
    double q[20][20]; //存放信道转移概率矩阵 
    double F[20][20]; //存放(|)iipxy的概率分布矩阵的转置  
    double x, y, a;
    double epsilon = 1e-5;  //门限 
    double C = -1000.0;    //取初始迭代时的信道容量为一个较大的负数  
    printMessage();
    printf("!-------请输入信源符号个数:---------!\n");
    scanf("%d", &r);
    printf("!-------请输入信宿符号个数：---------!\n");
    scanf("%d", &s);
    printf("!-------请输入信道转移概率矩阵---------!\n\n");

    for (i = 0; i < r; i++)
    {
        for (j = 0; j < s; j++)
            scanf("%lf", &q[i][j]);
        printf("\n");
    }

    for (i = 0; i < r; i++)
        p[i] = (double)(1.0 / (double)r);
    //设初始信源分布为等概分布  
    do
    {
        k++;
        a = C;
        for (j = 0; j < s; j++)
        {
            x = 0.0;
            for (i = 0; i < r; i++)
                x = x + (p[i])*(q[i][j]);// x 为(|)iipxy的分母部分    
            if (x > 0)
                for (i = 0; i < r; i++)
                    F[i][j] = (p[i])*(q[i][j]) / x;// F为(|)iipxy的概率分布矩阵的转置    
            else
                for (i = 0; i < r; i++)
                    F[i][j] = 0.0;//(|)iipxy的分母部分为0时，令F=0 
        }
        y = 0.0;
        for (i = 0; i < r; i++)
        {
            z[i] = 0.0;
            for (j = 0; j < s; j++)
            {
                if (F[i][j] > 0)
                    z[i] = z[i] + (q[i][j] * (log(F[i][j]) / log(2.0)));//z[i]为ip的分子部分 
            }
            z[i] = (pow(2.0, z[i]));
            y = y + z[i];//z[i]为ip的分母部分 
        }
        for (i = 0; i < r; i++)
        {
            p[i] = z[i] / y;//更新输入信源概率矩阵 
        }
        C = (log(y) / log(2.0));//求信道容量单位为“bit” 
    } while (fabs((C - a) / C) > epsilon);
    printf("!----------迭代次数为:k=%d  --------!\n", k);//输出迭代次数  
    printf("!----------最佳信源分布为：--------!\n\n");
    for (i = 0; i < r; i++)
    {
        printf("%.3lf", p[i]); //输出信源的最佳分布，保留3位小数 
    }
    printf("\n");
    printf("!----------信道容量为:C=%.3lf bit --------!\n\n", C); //输出信道容量，保留3位小数 
    getchar();
    getchar();
}
void printMessage()
{
    printf("!---------                --------!\n\n");
    printf("!----------姓名：小猪嘎嘎     --------!\n\n");
    printf("!----------信道容量迭代算法 --------!\n\n");
    printf("!----------csdn   --------!\n\n");
}

运行过程需要输入信源符号个数，信宿的个数和对应的状态转移矩阵。最终通过迭代过程获得信道容量和对应的最佳输入分布。