chcapacity2.cpp

信道容量迭代算法用C++ 通过给定的信道转移矩阵

#include<iostream.h>
#include<string>
#include<math.h>

double delta=1/10000;
void display(int ,int ,double**);
double **convert_to_fi(int ,int ,double **,double *);
double **log_fi(int,int,double **);
double *generate_newpr(int,int,double **,double **);
double channelcapacity(int,int,double **,double **);

void main()
{
	int r;
    int s;
	
	cout<<"请输入你的信源符号个数：\n";
	cin>>r;
	cout<<endl;
    cout<<"请输入你的信宿符号个数：\n";
	cin>>s;
	cout<<endl;
	int len=s*r;
	
	double **p;
    p=new double *[r];
	for(int i=0;i<r;i++)
	{
		p[i]=new double[s];
	}
	double **fi;
	fi=new double*[s];
	for(i=0;i<s;i++)
	{
		fi[i]=new double[r];
	}
	cout<<"下面输入信道转移概率矩阵：\n";
	for( i=0;i<r;i++)
	{
		cout<<"第"<<i+1<<"行：\n";
		for(int j=0;j<s;j++)
		{
			cout<<"第"<<j+1<<"个概率：";
			cin>>p[i][j];
			cout<<endl;        
		}
	}
	display(r,s,p);
	/*for(int ii=0;ii<r;ii++)
	{
		for(int jj=0;jj<s;jj++)
		{
			cout<<p[ii][jj]<<"   " ;
		}
		cout<<endl;
	}*/
	double *pr=new double[r];
    for(i=0;i<r;i++)
	{
		*(pr+i)=(double)1/r;
	}
	//double sum=0;
	//double *pxy=new double(r);
	fi=convert_to_fi(r,s,p,pr);
	display(s,r,fi);
	fi=log_fi(r,s,fi);
	display(s,r,fi);
	pr=generate_newpr(r,s,fi,p);
	//for( i=0;i<s;i++)
	//{
		for(int j=0;j<r;j++)
		{
		   cout<<pr[j]<<" ";
		}
	//	cout<<endl;
	//}
    double c=-1000000;
	double c_next=channelcapacity(r,s,fi,p);
	while((fabs(c_next-c)/c_next)>delta)
	{
        fi=convert_to_fi(r,s,p,pr);
        fi=log_fi(r,s,fi);
        pr=generate_newpr(r,s,fi,p);
		c=c_next;
        c_next=channelcapacity(r,s,fi,p);
		cout<<c_next<<endl;
	}
	cout<<c_next<<endl;
}
double **convert_to_fi(int r,int s,double **p,double *pr)
{
	double **fi;
	fi=new double*[s];
	for(int ii=0;ii<s;ii++)
	{
		fi[ii]=new double[r];
	}
	double sum=0;
	double *pxy=new double[r];
	for(int j=0;j<s;j++)
	{  
	   for(int i=0;i<r;i++)
	   {
		  sum=sum+p[i][j]*pr[i];
		  pxy[i]=p[i][j]*pr[i];
	   }
	   for( i=0;i<r;i++)
	   {
		   fi[j][i]=pxy[i]/sum;
	   }
       sum=0;
	}
    delete[]pxy;
	return fi;
	delete[]fi;
}
double **log_fi(int r,int s,double **fi)
{
    double **log_fi;
	log_fi=new double*[s];
	for(int i=0;i<s;i++)
	{
		log_fi[i]=new double[r];
	}
	for(i=0;i<s;i++)
	{
		for(int j=0;j<r;j++)
		{
			log_fi[i][j]=log(fi[i][j]);
		}
	}
	return log_fi;
	delete[]log_fi;
}
void display(int r,int s,double **metrix)
{
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<s;j++)
		{
		   cout<<metrix[i][j]<<" ";
		}
		cout<<endl;
	}
}
double *generate_newpr(int r,int s,double **loged_fi,double **p)
{
	double *new_pr;
	new_pr=new double[r];
    double sum1=0;
	double sum2=0;
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<s;j++)
		{
			sum1=sum1+p[i][j]*loged_fi[j][i];
		}
        sum1=exp(sum1);
		new_pr[i]=sum1;
		sum2=sum2+sum1;
		sum1=0;
	}
    for(i=0;i<r;i++)
	{
		new_pr[i]=new_pr[i]/sum2;
	}
	return new_pr;
	delete[]new_pr;
}

double channelcapacity(int r,int s,double **loged_fi,double **p)
{
	double *new_pr;
	new_pr=new double[r];
    double sum1=0;
	double sum2=0;
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<s;j++)
		{
			sum1=sum1+p[i][j]*loged_fi[j][i];
		}
        sum1=exp(sum1);
		new_pr[i]=sum1;
		sum2=sum2+sum1;
		sum1=0;
	}
	sum2=log(sum2);
	return sum2;
}