- #include <stdio.h>
- #include <math.h>
- #include <stdlib.h>
- void main()
- {
- register int i,j;
- register int k;
- int r,s; /*信道矩陣的行數和列數*/
- float *p=NULL;
- float *q=NULL;
- float *d=NULL;
- float *a=NULL;
- float **p_ji=NULL; /*二維數組p_yx存放信道矩陣*/
- float C,u,max,t,IU; /*C為信道容量*/
- printf("the number of input information source:\nr="); /*輸入行數和列數*/
- scanf("%d",&r);
- printf("the number of output information source:\ns=");
- scanf("%d",&s);
- p=(float *)calloc(r,sizeof(float)); /*初始化需要用到的數組*/
- q=(float *)calloc(r,sizeof(float));
- d=(float *)calloc(r,sizeof(float));
- a=(float *)calloc(r,sizeof(float));
- p_ji=(float **)calloc(r,sizeof(float));
- for(i=0;i<r;i++)
- p_ji[i]=(float *)calloc(r,sizeof(float));
-
- printf("Please matrix of channel:\n"); /*輸入信道矩陣*/
- for(i=0;i<r;i++)
- {
- for(j=0;j<s;j++)
- {
- printf("p_ji[%d][%d]=",i,j);
- scanf("%f",&p_ji[i][j]);
- }
- printf("\n");
- }
-
- for(i=0;i<r;i++) /*判斷輸入的信道矩陣是否符合概率分布*/
- {
- t=0.0;
- for(j=0;j<s;j++)
- t=t+p_ji[i][j];
- if(t<=0.999999||t>=1.000001)
- {
- printf("Error!");
- exit(-1);
- }
- }
-
- for(i=0;i<r;i++) /*取初始概率為均勻分布*/
- p[i]=1.0/r;
- for(j=0;j<s;j++) /*計算輸出信源概率分布q[j]*/
- {
- q[j]=0;
- for (i=0;i<r;i++)
- q[j]=q[j]+p[i]*p_ji[i][j];
- }
- for(i=0;i<r;i++) /*計算a[i]*/
- {
- d[i]=0;
- for(j=0;j<s;j++)
- {
- if(p_ji[i][j]==0)
- d[i]=d[i]+0;
- else
- d[i]=d[i]+p_ji[i][j]*log(p_ji[i][j]/q[j]);
- }
- a[i]=exp(d[i]);
- }
- u=0.0;
- for(i=0;i<r;i++) /*計算u*/
- u=u+p[i]*a[i];
- C=log(u)/log(2.0);
- max=a[0];
- for(i=0;i<r;i++)
- {
- if(a[i]>=max)
- max=a[i];
- }
- IU=log(max)/log(2.0); /*計算IU*/
-
- k=1;
- while((IU-C)>=10e-6) /*迭代計算*/
- { /*計算最佳輸入分布p[i]*/
- for(i=0;i<r;i++)
- p[i]=p[i]*a[i]/u;
- for(j=0;j<s;j++) /*計算輸出信源的概率分布q[j]*/
- {
- q[j]=0;
- for (i=0;i<r;i++)
- q[j]=q[j]+p[i]*p_ji[i][j];
- }
- for(i=0;i<r;i++) /*計算a[i]*/
- {
- d[i]=0;
- for(j=0;j<s;j++)
- {
- if(p_ji[i][j]==0)
- d[i]=d[i]+0;
- else
- d[i]=d[i]+p_ji[i][j]*log(p_ji[i][j]/q[j]);
- }
- a[i]=exp(d[i]);
- }
- u=0.0;
- for(i=0;i<r;i++) /*計算u*/
- u=u+p[i]*a[i];
- C=log(u)/log(2.0); /*計算IL*/
- max=a[0];
- for(i=0;i<r;i++)
- {
- if(a[i]>=max)
- max=a[i];
- }
- IU=log(max)/log(2.0); /*計算IU*/
- k=k+1;
- }
-
- printf("The matrix of channel is:\n"); /*輸出信道矩陣*/
- for(i=0;i<r;i++)
- {
- for(j=0;j<s;j++)
- printf("%.3f ",p_ji[i][j]);
- printf("\n");
- }
- printf("The iteration number is:\nk=%d\n",k); /*輸出迭代次數*/
- printf("The best input probability distribution is:\n"); /*計算輸出信源的概率分布q[j]*/
- for(i=0;i<s;i++)
- printf("%.3f ",q[i]);
- printf("\nThe capacity of the channel is:\nC=%fbit/symbol\n",C); /*輸出信道容量*/
- printf("The best input probability distribution is:\n"); /*輸出最佳輸入分布*/
- for(i=0;i<r;i++)
- printf("%.3f ",p[i]);
- }
|
QQ好友和群
QQ空間
騰訊微博
騰訊朋友
收藏
淘帖
頂
踩
