/*
 * Tabulate bit-error-rate for example from section 2 of
 * "Sequence Estimation EM Examples", R. Perry, June 1998
 */

SNR_min = -6;	/* min, max, and number of SNR points to use */
SNR_max = 10;
Npts = 17;

Ntrials = 10;	/* number of trials for each SNR value */

NN = 100;	/* number of samples in each trial */
pp = 0.8;	/* P(b==1) */
cc = 2;		/* scale factor */

ds = (SNR_max - SNR_min) / (Npts - 1);	/* delta SNR */

known = BER = SS = zeros(1,Npts);	/* arrays for SNR values and BER results */

verbose = 0;	/* to avoid messages from iter() */

(void) cleanup(1000);	/* run faster, but use more memory */

for( i = 1; i <= Npts; ++i) {	/* for each SNR */

  SS(i) = SNR_min + (i-1)*ds;

  printf( "i = %g, SNR = %g\n", i, SS(i));

  for( j = 0; j < Ntrials; ++j) {	/* for each trial */

    get_y( SS(i), pp, cc, NN);

   /*
    * Simulate the case of known channel parameter c
    * by initializing b to the original generated values
    * so that the M-step will get the best possible
    * estimate of c.
    */
    b = b_save;
    iter(1);
    count_known = sum(b!=b_save);	/* number of bit errors */
    known(i) = known(i) + count_known;

   /*
    * Now do the regular case with b initialized to all 1's
    */
    b = ones( 1, NN);
    iter(3);
    count = sum(b!=b_save);		/* number of bit errors */
    BER(i) = BER(i) + count;

    printf( "  j = %g, count_known = %g, count = %g\n", j, count_known, count);
  }

  known(i) = known(i) / (NN * Ntrials);	/* normalize by total number of bits */
  BER(i) = BER(i) / (NN * Ntrials);
}
