compute_c.cpp

本人学习“信息论”的时候的一个作业程序

#include "iostream.h"
#include "math.h"
#include "stdio.h"

void main()
{
	
	const int  R=4;//设置信源的符号个数
	const int  S=4;//设置信宿的符号个数
	const double  LIMIT=0.000001;//当C(n+1,n)迭代计算相邻两次的差达到此下限，即认为迭代结束
	
	double pba[R][S]={0.50,0.25,0.00,0.25,
		              0.00,1.00,0.00,0.00,
					  0.00,0.00,1.00,0.00,
					  0.25,0.00,0.25,0.50};

	double pab[S][R]={0},pb[S]={0},pa[R]={0},e[R]={0};
   	double sum1=0.0,sum2=0.0,sumei,c,tempc,difference;
	int i,j,k,refer;

   /******输入初始的信源概率分布*************************/
	do
	{
		sum1=0;
		cout<<"输入初始的信源的概率分布pa[i],有"<<R<<"个:"<<endl;
		for(i=0;i<R;i++)
		{
			cout<<"  pa["<<i<<"] = ";
			cin>>pa[i];
		}
		for(i=0;i<R;i++)
		{
			sum1+=pa[i];
		}
		if(fabs(sum1-1.0)>=0.0001)
		{
			cout<<"输入错误：你输入的pa[i]，全部垒加和不等于1，请重新输入。\n"<<endl;
			refer=1;
		}
		else
		{
			refer=0;
		}

	}while(refer);

	/*******打印信道矩阵***********************************/
	cout<<"这个固定信道的信道矩阵是："<<endl;
	for(i=0;i<R;i++)
	{
		for(j=0;j<S;j++)
		{
			printf(" %1.2f ",pba[i][j]);
			//cout<<pba[i][j]<<"  ";
		}
		printf("\n");
	}
	
	/************************************************************/
	/*****迭加算法开始*******************************************/
    
	do
	{
		sumei=0;
		c=0;
		tempc=0;

		/********************计算pab[i][j]****************/
		for(j=0;j<S;j++)
		{
			for(i=0;i<R;i++)
			{
				if(pba[i][j]!=0.00)   
				{
					for(k=0;k<R;k++)
					{
						pb[j]+=pa[k]*pba[k][j];
					}
					pab[j][i]=pa[i]*pba[i][j]/pb[j];
					pb[j]=0.0; 
				}
				
			}
		}
		/*****************计算信道容量C(n,n) ****************/
		for(j=0;j<S;j++)
		{
			for(i=0;i<R;i++)
			{
				if (pab[j][i]!=0.00)  
					tempc+=pa[i]*pba[i][j]*log(pab[j][i]/pa[i]);
			}
		}

		/*********************计算pa[i]********************/
		/*****先计算e[i]***********/
		for(i=0;i<R;i++)
		{
			for(j=0;j<S;j++)
			{
				if(pab[j][i]!=0.00)  
					sum2+=pba[i][j]*log(pab[j][i]);
			}
			e[i]=exp(sum2);
			sum2=0.0;
		}

		for(i=0;i<R;i++) {	sumei+=e[i];}       //再计算e[i]全部和
		for(i=0;i<R;i++) {	pa[i]=e[i]/sumei;}  //最后综合算出pa[i]

		/*****************计算信道容量C(n+1,n) ****************/
		for(j=0;j<S;j++)
		{
			for(i=0;i<R;i++)
			{
				if (pab[j][i]!=0.0)  
					c+=pa[i]*pba[i][j]*log(pab[j][i]/pa[i]);	
			}
		}
		difference=fabs(c-tempc);//计算相邻两次迭代的差值 

	}while(difference>=LIMIT);

    /******结果输出******************************************/
	cout<<endl;
	cout<<"这个固定信道的匹配信源是："<<endl;
	for(i=0;i<R;i++)
	{
		cout<<"  "<<pa[i]<<endl;
	}
           
	cout<<"匹配时的后验概率矩阵是："<<endl;
	for(j=0;j<S;j++)
	{
		for(i=0;i<R;i++)
			{
			   printf("  %6f  ",pab[j][i]);
			}
		printf("\n");   
	}
	c=c/log(2.0);//把容量转换为以比特为单位
	cout<<endl;
    cout<<"这个固定信道的信道容量C= "<<c<<"  比特"<<endl;	 //打印出信道容量
}