/*
 * Tabulate bit-error-rate for example from section 1 of
 * "Time Recursive EM Examples", R. Perry, June 1998
 */

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

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

SS = [3 6 9 12];	/* SNR values */

BER = zeros(5,4);

/*
 * BER rows will be:
 *  known
 *  block
 *  K=100
 *  K=50
 *  K=1
 *
 * The case of known channel parameter c is simulated by
 * initializing b to the original generated values so that
 * the M-step will get the best possible estimate of c.
 */

/*
 * This is coded inline rather than using loops to get around
 * some memory management quirks (not quite bugs!) in M.
 */

J=1
print(J)
get_y( SS(J), pp, cc, NN)
b = b_save
iter(1)
BER(1,J) = sum(b!=b_save) / NN
b = ones( 1, NN)
iter(3)
BER(2,J) = sum(b!=b_save) / NN
K=100
iter_r
BER(3,J) = sum(b!=b_save) / NN
K=50
iter_r
BER(4,J) = sum(b!=b_save) / NN
K=1
iter_r
BER(5,J) = sum(b!=b_save) / NN

J=2
print(J)
get_y( SS(J), pp, cc, NN)
b = b_save
iter(1)
BER(1,J) = sum(b!=b_save) / NN
b = ones( 1, NN)
iter(3)
BER(2,J) = sum(b!=b_save) / NN
K=100
iter_r
BER(3,J) = sum(b!=b_save) / NN
K=50
iter_r
BER(4,J) = sum(b!=b_save) / NN
K=1
iter_r
BER(5,J) = sum(b!=b_save) / NN

J=3
print(J)
get_y( SS(J), pp, cc, NN)
b = b_save
iter(1)
BER(1,J) = sum(b!=b_save) / NN
b = ones( 1, NN)
iter(3)
BER(2,J) = sum(b!=b_save) / NN
K=100
iter_r
BER(3,J) = sum(b!=b_save) / NN
K=50
iter_r
BER(4,J) = sum(b!=b_save) / NN
K=1
iter_r
BER(5,J) = sum(b!=b_save) / NN

J=4
print(J)
get_y( SS(J), pp, cc, NN)
b = b_save
iter(1)
BER(1,J) = sum(b!=b_save) / NN
b = ones( 1, NN)
iter(3)
BER(2,J) = sum(b!=b_save) / NN
K=100
iter_r
BER(3,J) = sum(b!=b_save) / NN
K=50
iter_r
BER(4,J) = sum(b!=b_save) / NN
K=1
iter_r
BER(5,J) = sum(b!=b_save) / NN

