一般信道容量的计算

#include<stdio.h>
#include<math.h>
#define MAX 9
int main()
{
/*定义数组及中间变量*/
double fun(int n,double g[MAX][MAX]);
double b[MAX]={0},d[MAX]={0},t[MAX][MAX]={{0}},g[MAX]={0},h[MAX]={{0}};
double a[MAX][MAX]={{0}},p[MAX][MAX]={{0}},f[MAX][MAX]={{0}};
double c=0,q=0,m=0,o=0;
int i=0,j=0,k=0,n=0;
/*输入计算所需数值*/
printf("the number of input&output channel symbol:");
scanf("%d",&n);
printf("\nthe channel transfer matrix is:\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
scanf("%lf",&p[i][j]);  //p[][]为信道传递矩阵的元素
if(p[i][j]>0)
b[i]+=p[i][j]*(log(p[i][j])/log(2));  //b[]为线性方程组的常数项
}
}
q=fun(n,p);  //利用克拉默法则求解方程组，q为系数行列式D的解
for(j=0;j<n;j++)  //此处循环得系数行列式D中第j列的元素用常数项代替后所得到的n阶行列式
{   for(k=0;k<n;k++)
{
for(i=0;i<n;i++)
{
if(k==j)
a[i][k]=b[i];
else
a[i][k]=p[i][k];
}
}

d[j]=fun(n,a)/q;  //求解βi
m+=pow(2,d[j]);
}
c=log(m)/log(2);  //求解信道容量
printf("the channel capacity is:%lf\n",c);
printf("Output Probability Distribution:\n");
for(i=0;i<n;i++)
{
g[i]=pow(2,d[i]-c);  //求解输出概率分布P(bj)
printf("%7.2lf\t",g[i]);
}
printf("\n");
for(i=0;i<n;i++)    //将信道传递矩阵转置，以求p(ai)
for(j=0;j<n;j++)
f[i][j]=p[j][i];
o=fun(n,f);
for(j=0;j<n;j++)    //此处for循环用于表示求解P[ai]、P[bj|ai]、P(bj)组成的线性方程组所需的各个矩阵
{   for(k=0;k<n;k++)
{
for(i=0;i<n;i++)
{
if(k==j)
t[i][k]=g[i];
else
t[i][k]=f[i][k];
}
}
h[j]=fun(n,t)/o;  //计算输入概率分布P(ai)
}
printf("Input probability distribution:\n");
for(i=0;i<n;i++)
printf("%7.2lf\t",h[i]);
printf("\n");
return 0;
}

double fun(int n,double g[MAX][MAX])       //求n阶矩阵的值
{
double e[MAX][MAX]={{0}};
double z=0,x=0;
int i=0,j=0,k=0;
if(n==1)
return g[0][0];
for(i=0;i<n;i++)                      //此循环用于表示余子式
{
for(k=0;k<n-1;k++)
{   for(j=0;j<n-1;j++)
{
if(k<i)
e[k][j]=g[k][j+1];
else
e[k][j]=g[k+1][j+1];
}
}
if(i%2==0)
x=1;
else
x=(-1);
z+=g[i][0]*fun(n-1,e)*x;   //z为n阶矩阵的值
}
return z;
}