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

disp('uncorrelated h''s');

disp('Gauss EM and ML')

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.995
alpha2 = alpha*alpha;

% 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)
%
%
 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.1
  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

%%% high correlation:
%%
%%  inv_T = tril(ones(FIR_length,FIR_length))/20;
%%  inv_T(:,1) = 1;
%%
%%  T = inv(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


SNR = 10

trials = 10000

% EM iteration counters
%
iter_avg_g = zeros(N+1,1);
iter_max_g = zeros(N+1,1);

EM_equal_ML = zeros(N+1,1);

EM_local_min = zeros(N+1,1);

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

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

    % generate data
    %
    x4;

    % set pilot bits
    %
    % pilot(pilot_index) = x(pilot_index);
 
 if(1)
    %% direct ML

    ML = 1;

    % for Gauss-Markov h
    %
    if( GM_H)

      b = get_GM_ML;

    else
    %
    %% for Gaussian h
    
      G=[];
      vit_C;
      b = get_ex;

      %%% for p=1:N, if( pilot(p) & (b(p) ~= x(p))) disp('pilot mismatch!'); end; end;
    end

    b_ML = b;

    ML = 0;
 end

    %% run EM for time-varying Gaussian or Gauss-Markov h
    
    GAUSS = 1;	% GM or GAUSS

    % init using ML results
    %
     b_ML_delta = b_ML;

   for b_ML_index = 1:N+1

    if( b_ML_index > 1)
      % b_ML_delta(b_ML_index-1) = - b_ML_delta(b_ML_index-1);
      b_ML_delta(N-b_ML_index+2) = - b_ML_delta(N-b_ML_index+2);
    end

    b = b_ML_delta;

    b_temp = b;
    iter = 0;
    while(1)
        iter = iter + 1;
        B = get_A_C( b);
        vit_C;
        b = get_ex;
      if( all(b == b_temp))
        break;
      end
      b_temp = b;
    end

    if( all( b == b_ML))
      EM_equal_ML(b_ML_index) = EM_equal_ML(b_ML_index) + 1;
    end

    if( GM)
      r = check_GM_ML_min( b);
    else
      r = check_Gauss_ML_min( b);
    end

    if( r)
      EM_local_min(b_ML_index) = EM_local_min(b_ML_index) + 1;
    % else
    %  disp('EM != local min');
    %  return;
    end

    iter_avg_g(b_ML_index) = iter_avg_g(b_ML_index) + iter;

    if( iter > iter_max_g(b_ML_index))
      iter_max_g(b_ML_index) = iter;
    end

   end % b_ML_index

    GAUSS = 0;	% GM or GAUSS

  end % trials
  
iter_avg_g = iter_avg_g ./ trials
iter_max_g

disp('EM_local_min %');

EM_local_min = 100*EM_local_min./trials

disp('EM_equal_ML %');

EM_equal_ML = 100*EM_equal_ML./trials
