% for varying alpha_model

alpha_model = [0 0.25 0.5 0.75 0.9 0.95 0.96 0.97 0.975 0.98 0.99]

trials = 500

% disp('using N-FIR_length bits to compute BER');
disp('using N bits to compute BER');

disp('uncorrelated h''s');

disp('GM h');

global KNOWN;	KNOWN = 0;	% flags for known channel,
global DET;	DET = 0;	%  deterministic EM,
global GAUSS;	GAUSS = 0;	%  Gaussian EM,
global GM;	GM = 0;		%  Gauss-Markov EM,
global ML;	ML = 0;		%  direct ML,
global KALEH;	KALEH = 0;	%  Kalet-Valet

global GAUSS_H;	GAUSS_H = 0;	% flags for Gaussian h,
global GM_H;	GM_H = 0;	%  Gauss-Markov h

global pilot;		pilot=[];	% for pilot bits
global pilot_index;	pilot_index=[];

global B;
global var_noise;
global Rh;
global inv_Rh;
global inv_Rh_d;
global avg_h;
global FIR_save;
global FIR_g;
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.975
alpha2 = alpha*alpha;

alpha_save = alpha;
alpha2_save = alpha2;

% for Kaleh and Vallet simulation
%
%
global PS;      PS=[];  % state transition probability matrix
global Xt;              % X(t) vector values
global Nij;             % used to get i and j given n
global IJn;             % used to get n given i and j
global gamma;           % array for forward recursion
global beta;            % array for backward recursion
global LY;              % array for likelihood values
global gamma_scale;     % normalization for columns of gamma and beta

% initialize rand() and randn() with time-of-day
%
  seed=floor(sum(100*clock))
%
%  disp('Using seed from gm.04')
% 
%  seed = 210675

  rand( 'state', seed);
  randn( 'state', seed);

% set FIR length
%
  FIR_length = 3

% set N 
%
  N = 10

% for Gaussian FIR or the Gaussian part of Gauss-Markov
%
% Gauss-Markov
%
  GM_H = 1;

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

%  FIR_0 = ones(FIR_length,1)

%
%  avg_h = -ones(FIR_length,1)*(1-alpha)
%
%
% GAUSS_H = 1
%
% Gaussian with time-varying mean 1 -> 2
%
%  avg_h = ones(FIR_length,1) * (1 + [0:N-1]*((2-1)/(N-1)))
%
% Gaussian with time-varying mean -1 -> 1  % doesn't work for N=1
%
%  avg_h = ones(FIR_length,1) * (-1 + [0:N-1]*(2/(N-1)))
%
% Gaussian with mean 0 or 1
%
%  avg_h = ones(FIR_length,N)

  var_h = 0.001
  stdev_h = sqrt(var_h);

% uncorrelated h's
%
  T = eye(FIR_length,FIR_length)
%
% correlated h's
%
%  (h-g)' inv_Rh (h-g) / 2  =  (h-g)' T' T (h-g) / 2v   with inv_Rh = T' T / v
%
%   = z' z / 2v   with z = T (h-g)  where z is uncorrelated with mean 0
%
% So  h = inv_T z + g
%
% Choose T so that det(T) = 1, then det(Rh) == det(I/v)
%

% low correlation:
%
%  T = tril(ones(FIR_length,FIR_length))

  inv_T = inv(T)
  inv_Rh = (T'*T)/var_h

  Rh = inv(inv_Rh)

  if( GM_H)
    inv_Rh_d = inv_Rh*avg_h;	% inv(C)*d in the paper
  end

pilot = []; % zeros(N,1);
pilot_index = []; % [1:10:91, 100]

SNR = 10

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

  % run the Monte-Carlo simulations
  %
  trial_count = 0;
  
  x6_init;

  while( trial_count < trials)
  
    trial_count = trial_count + 1;

    % generate data
    %
    alpha = alpha_save;
    alpha2_save = alpha2_save;
    x6;

    % known channel case
    %
    G=[];
    KNOWN = 1;

    vit_C;

    b = get_ex;

    wrong = (b ~= x');
    %
    % kwrong = kwrong + sum(wrong(1:N-FIR_length));
    %
    kwrong = kwrong + sum(wrong);
  
    KNOWN = 0;


    %% direct ML

    ML = 1;

    for i=1:length(alpha_model)

      alpha = alpha_model(i);
      alpha2 = alpha*alpha;

      b = get_GM_ML;

      wrong = (b ~= x');
      %
      % gwrong(i) = gwrong(i) + sum(wrong(1:N-FIR_length));
      %
      gwrong(i) = gwrong(i) + sum(wrong);

    end % alpha_model

    ML = 0;

  end % trials
  
  % display results
  %
  % kwrong = kwrong / ((N-FIR_length)*trials);
  % gwrong = gwrong / ((N-FIR_length)*trials);
  %
  kwrong = kwrong / (N*trials);
  gwrong = gwrong / (N*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

