C program implementation

/*****
Last Modification: 6 November 2006
****
****
N.B.: Use formula (51) on page 632 of Shannon's 1959 paper,
See also 
http://www.comelec.enst.fr/~boutros/coding/Allerton99_VialleBoutros_paper.pdf
*****/

#include "mathcom.h"

static int n;
static double RateFunction();


int main(argc, argv)
int argc;
char **argv;
{
 char chain[200], filename[200];
 double theta, rate;
 double a, b, c;
 int iter;
 double val, diff;
 double A, G, EL, factor, Pe, log_Pe, snrdb, snr;
 double snrdb1, snrdb2, snrstep;
 
 printf("######## Optimal Code (Spherical) Performance on AWGN Channel ######### \n");
 printf("Q(theta) formula (51) of Shannon 1959, page 632. \n");
 printf("Using the asymptotic rate formula (28), page 624 Shannon 1959\n");
 printf("Computing the cone half-angle for a given rate\n");
 printf("R=(1/n-1)*log2(sin(theta))+1/n*log2(cos(theta)*sqrt(2PIn)) \n\n");


 printf("Real Space dimension (from 20 up to 200000) [1000] : ");
 ReadString(chain);
 if(chain[0]) sscanf(chain,"%d", &n);
 else n=1000;

 if((n<20)||(n>200000))
   {
    fprintf(stderr,"%s: n out of range 100...200000\n", argv[0]);
    return(-1);
   }

 printf("Information Rate in bits/dimension (from 1/10 up to 8) [0.50] : ");
 ReadString(chain);
 if(chain[0]) sscanf(chain,"%lf", &rate);
 else rate=0.50;

 if((rate<0.1)||(rate>8.0))
   {
    fprintf(stderr,"%s: rate out of range 1/10...8\n", argv[0]);
    return(-1);
   }

 /* rate=1/10 <=> theta0=00.20 degrees */
 /* rate=8    <=> theta0=69.50 degrees */
 iter=0;
 a=0.20;  b=70.0; 
 do
   {
    c=(a+b)/2.0;
    val=RateFunction(c);
    if(val < rate) b=c; else a=c;
    c=(a+b)/2.0; ++iter;
    diff=b-a;
   }
 while((diff>1e-8)&&(iter<1000));

 printf("The cone half-angle for %1.4f bits/dimension is equal to %1.8f degrees \n", rate, c);

 theta=c*PI/180.0; /* we switch to radians */

 if(SQR(tan(theta)) >= (0.25*n) )
   {
    fprintf(stderr,"OmegaFunction(): Warning, theta too close to 90 or n is small !\n");
   }

 printf("The modulation alphabet is spherical, Es=energie per real symbol.\n");
 printf("n*P=squared radius of the sphere= n*2*Es=n*2*rate*Eb, P=2*rate*Eb \n");
 printf("A^2=P/N0=2*rate*Eb/N0, the SNR per bit is SNR=Eb/N0 \n");
 printf("Here, N0 is the noise variance per real component.\n\n");

 snrdb1=-1.50;
 printf("Start Eb/N0 in dB [%1.2f]: ", snrdb1);
 ReadString(chain);
 if(chain[0]) sscanf(chain,"%lf", &snrdb1);

 snrdb2=20.00;
 printf("End Eb/N0 in dB [%1.2f]: ", snrdb2);
 ReadString(chain);
 if(chain[0]) sscanf(chain,"%lf", &snrdb2);

 snrstep=0.05;
 printf("SNR step in dB [%1.2f]: ", snrstep);
 ReadString(chain);
 if(chain[0]) sscanf(chain,"%lf", &snrstep);


 sprintf(filename,"optimal_spherical_code_awgn_n=%d_rate=%1.2f_WER.dat", n, rate);
 printf("Output file name for WER versus SNR [%s]: ", filename);
 ReadString(chain);
 if(chain[0]) strcpy(filename, chain);

 sprintf(chain,"rm -f %s", filename);
 system(chain);

 printf("Eb/N0 (dB)  ||  Word Error Rate (WER)\n");
 Pe=1.0;
 for(snrdb=snrdb1; snrdb <= snrdb2+0.01; snrdb += snrstep)
   {
    if(Pe<1e-20) break;
    snr=exp10(0.1*snrdb);
    A=sqrt(2*rate*snr);
    if(atan(1.0/A) > theta)
      {
       /***
       fprintf(stderr,"Warning : theta=%1.4f but lower limit is cot^-1(1/A)=%1.4f ! \n", 
               theta*180.0/PI, atan(1.0/A)*180.0/PI);
       fprintf(stderr,"%s: Please increase the SNR for this rate=%1.4f.\n", argv[0], rate);
       ***/
       continue;
      }

    G=0.5*(A*cos(theta)+sqrt(A*A*cos(theta)*cos(theta)+4.0));
    EL=0.5*A*A-0.5*A*G*cos(theta)-log(G*sin(theta));
    factor=sqrt(n*PI)*sqrt(1.0+G*G)*sin(theta)
           *(A*G*sin(theta)*sin(theta)-cos(theta));
    Pe=exp(-n*EL)/factor;
    if(Pe<1.0) 
      {
       printf("%1.2f    %1.4e \n", snrdb, Pe);
       WritePlotFile(filename, snrdb, Pe, "%1.2f", "%1.4e");
      }
   }/* end of snrdb loop */

 return(0);
}/* end of main() */

/*---------------------------------------------------------------*/

static double RateFunction(theta)
double theta; /* in degrees, not radians */
{
  return( (1.0/n-1.0)*log2(sin(theta*PI/180.0))
          +1.0/n*log2(cos(theta*PI/180.0)*sqrt(2.0*PI*n)) );
}/* end of RateFunction() */

/*---------------------------------------------------------------*/