% vary var_h model

trials = 1000

disp('using N bits to compute BER');

var_h_model = [0.01 0.03 0.1 0.3 1.0]

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

global B;
global var_noise;
global Rh;
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;

% 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))
%
% seed from gauss.t.08
%
%  seed = 207608

  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
%
%  avg_h = ones(FIR_length,1)*(1-alpha)		% so avg_h/(1-alpha) = all 1's
%
% Gaussian with time-varying mean 1 -> 2
%
  avg_h = ones(FIR_length,1) * (1 + [0:N-1]*((2-1)/(N-1)))
%
% Gaussian with mean 0
%
%  avg_h = zeros(FIR_length,N)

  var_h = 0.1
  stdev_h = sqrt(var_h);

  var_h_save = var_h;
  stdev_h_save = stdev_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)
%
  T = tril(ones(FIR_length,FIR_length))

% set inv_T and inv_Rh

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

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

  inv_Rh_save = inv_Rh;
  Rh_save = Rh;

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


SNR = 10

% EM iteration counters
%
iter_avg_det = 0;
iter_max_det = 0;

iter_avg_g = zeros(length(var_h_model),1);
iter_max_g = zeros(length(var_h_model),1);

 
  dwrong = 0;
  gwrong = zeros(length(var_h_model),1);
  kwrong = 0;
  Nwrong = 0;

  trial_count = 0;
  
  while( trial_count < trials)
  
    trial_count = trial_count + 1;
  
    % generate data
    %
    var_h = var_h_save;
    stdev_h = stdev_h_save;
    inv_Rh = inv_Rh_save;
    Rh = Rh_save;

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

    vit_C;

    b = get_ex;

    kwrong = kwrong + sum(b ~= x');
  
    KNOWN = 0;

    %% run EM for deterministic h
    
    G=[];
    B=[];
    DET = 1;

    vit_C;
    b = get_ex;

    b_temp = b;
    iter = 0;
    while(1)
      iter = iter + 1;
        B = get_A( b, N, Memory);
        vit_C;
        b = get_ex;
      if( all(b == b_temp))
        break;
      end
      b_temp = b;
    end

    iter_avg_det = iter_avg_det + iter;

    if( iter > iter_max_det)
      iter_max_det = iter;
    end

    dwrong = dwrong + sum(b ~= x');
 
    DET = 0;

    %% direct ML
    
    G=[];
    ML = 1;

    vit_C;
    b = get_ex;

    Nwrong = Nwrong + sum(b ~= x');
 
    ML = 0;

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

    for i = 1:length(var_h_model)
  
      var_h = var_h_model(i);
  
      inv_Rh = (T'*T)/var_h;
  
      Rh = inv(inv_Rh);
  
      B=[];
  
      vit_C;
      b = get_ex;
  
      b_temp = b;
      iter = 0;
      while(1)
        iter = iter + 1;
          B = get_A( b, N, Memory);
          vit_C;
          b = get_ex;
        if( all(b == b_temp))
          break;
        end
        b_temp = b;
      end
  
      iter_avg_g(i) = iter_avg_g(i) + iter;
  
      if( iter > iter_max_g(i))
        iter_max_g(i) = iter;
      end
  
      gwrong(i) = gwrong(i) + sum(b ~= x');
  
    end

    GAUSS = 0;	% or GM

  end % trials
  
  % display results
  %
  dwrong = dwrong / (N*trials);
  kwrong = kwrong / (N*trials);
  gwrong = gwrong / (N*trials);
  Nwrong = Nwrong / (N*trials);

  disp('[var_h_model Det Gauss Known ML]');

  for i=1:length(var_h_model)

    fprintf( 1, '%g\t%g\t%g\t%g\t%g\n', [var_h_model(i) dwrong gwrong(i) kwrong Nwrong]);

  end

iter_avg_det = iter_avg_det / trials;
iter_avg_det
iter_max_det

iter_avg_g = iter_avg_g / trials;
iter_avg_g
iter_max_g
