% for varying alpha_model

trials = 10000

alpha_model = [0 0.25 0.5 0.7 0.75 0.8 0.85 0.9 0.95]

disp('using N-FIR_length bits to compute BER');
disp('using 4 EM iterations for all algorithms');

global B;
global var_noise;
global inv_Rh;
global inv_Rh_d;
global avg_h;
global FIR_save;
global FIR_length;

global GG;
global qq;

% alpha is the 1st-order Markov factor such that alpha = exp(-wT)
% where B = w/pi is the Doppler spread
%
global alpha;
global alpha2;

% 1/(1-0.95) = 20
% 1/(1-0.95*0.95) = 10.256
%
alpha = 0.8
alpha2 = alpha*alpha;

alpha_save = alpha;
alpha2_save = alpha2;

global KALEH;		% flag used in cost_init
KALEH = 0;

% initialize rand() and randn() with time-of-day
%
  seed=floor(sum(100*clock))
  rand( 'state', seed);
  randn( 'state', seed);

% set FIR length
%
  FIR_length = 3;

% set N 
%
  N = 10;

% Gauss-Markov
%
  avg_h = ones(FIR_length,1)*(1-alpha)		% so avg_h/(1-alpha) = all 1's

  var_h = 0.1
  stdev_h = sqrt(var_h);

  T = tril(ones(FIR_length,FIR_length))

% set inv_T, inv_Rh, and inv_Rh_d
%
  inv_T = inv(T);
  inv_Rh = (T'*T)/var_h

  if( size(avg_h,2) == 1)
    inv_Rh_d = inv_Rh*avg_h;	% inv(C)*d in the paper
  end

SNR = 6

kwrong = 0;
gwrong = zeros(size(alpha_model));

  % run the Monte-Carlo simulations
  %
  trial_count = 0;
  
  while( trial_count < trials)
  
    trial_count = trial_count + 1;
  
    % generate data
    %
    alpha = alpha_save;
    alpha2_save = alpha2_save;
    x4;
  
    % known channel case
    %
    G=[];
  
    x4_em_init;

    b_save = b;

    wrong = (b ~= x');
    kwrong = kwrong + sum(wrong(1:N-FIR_length));
  
    % run EM for Gauss-Markov h
    %
    for i=1:length(alpha_model)

      alpha = alpha_model(i);
      alpha2 = alpha*alpha;
      b = b_save;
      G=1;

      for iter=1:4
        x4_em_iter;
      end
  
      wrong = (b ~= x');
      gwrong(i) = gwrong(i) + sum(wrong(1:N-FIR_length));

    end % for alpha_model(i)

  end % trials
  
  % display results
  %
  kwrong = kwrong / ((N-FIR_length)*trials);
  gwrong = gwrong / ((N-FIR_length)*trials);

  disp('[alpha_model Gauss GM Known]');

  for i=1:length(alpha_model)

    fprintf( 1, '%g %g %g %g\n', [alpha_model(i) gwrong(1) gwrong(i) kwrong]);

  end

