trials = 1

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

disp('uncorrelated h''s')

% disp('initial FIR using steady-state statistics')

%disp('DET and GM init [0 1 0]')

disp('DET and Gauss init using avg_h')

% disp('DET and GM init using FIR_save(:,1)')
% disp('GM update using alpha and avg_h')

% disp('using bad EM init');
% disp('high correlation');

% disp('GM h and ML')

% disp('DET FIR init using [0 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 pilot;	pilot=[];	% for pilot bits

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.95
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.25')
% 
%  seed = 214135

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

  avg_h = [1 0 -1]' * (-1 + [0:N-1]*(2/(N-1)));
  avg_h((FIR_length+1)/2,:) = 1

% Gaussian with mean 0
%
%  avg_h = zeros(FIR_length,N)

  var_h = 0.05
  stdev_h = sqrt(var_h);

% uncorrelated h's
%
  T = eye(FIR_length,FIR_length)
  inv_T = T
  inv_Rh = T/var_h

%
% 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( size(avg_h,2) == 1)
    inv_Rh_d = inv_Rh*avg_h;	% inv(C)*d in the paper
  end


% just testing x4
%
  SNR = 10;
  x4;
  return;

% for constant SNR and time-varying noise
%
SNR_array = [0:5:20]

% EM iteration counters
%
iter_avg_det = zeros(length(SNR_array),1);
iter_max_det = zeros(length(SNR_array),1);

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


% pilot = zeros(1,N);

for SNR_index = 1:length(SNR_array)

  % run the Monte-Carlo simulations
  %
  SNR = SNR_array(SNR_index);
 
  dwrong_where = zeros(N,1);
  gwrong_where = zeros(N,1);
  kwrong_where = zeros(N,1);
  Nwrong_where = zeros(N,1);

  % Kwrong_where = zeros(N,1);	% for kaleh

  trial_count = 0;
  
  while( trial_count < trials)
  
    trial_count = trial_count + 1;
  
    % generate data
    %
    x4;

    % set pilot bits
    %
    % pilot(N) = x(N);
    % pilot(N/2) = x(N/2);
    % pilot(1) = x(1);
 
    % known channel case
    %
    G=[];
    KNOWN = 1;

    vit_C;

    b = get_ex;

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

    b_save = b;

    wrong = (b ~= x');
    kwrong_where = kwrong_where + wrong;
  
    KNOWN = 0;

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

    % reasonable init:
    %
    B=[];
    vit_C;
    b = get_ex;

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

    % bad init:
    %
    % b_bad = 2*(rand(N,1) > 0.5) - 1; % random [-1,+1]
    % b = b_bad;

    b_temp = b;
    iter = 0;
    while(1)
      iter = iter + 1;
        B = get_A( b, N, Memory);
        vit_C;
        b = get_ex;
        %%% for p=1:N, if( pilot(p) & (b(p) ~= x(p))) disp('pilot mismatch!'); end; end;
      if( all(b == b_temp))
        break;
      end
      b_temp = b;
    end

    iter_avg_det(SNR_index) = iter_avg_det(SNR_index) + iter;

    if( iter > iter_max_det(SNR_index))
      iter_max_det(SNR_index) = iter;
    end

    wrong = (b ~= x');
    dwrong_where = dwrong_where + wrong;
 
    DET = 0;

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

    % reasonable init:
    %
    B=[];
    vit_C;
    b = get_ex;

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

    % bad init:
    %
    % b = b_bad;

    % init with known channel results
    %
    % b = b_save;

    b_temp = b;
    iter = 0;
    while(1)
      iter = iter + 1;
        B = get_A( b, N, Memory);
        vit_C;
        b = get_ex;
        %%% for p=1:N, if( pilot(p) & (b(p) ~= x(p))) disp('pilot mismatch!'); end; end;
      if( all(b == b_temp))
        break;
      end
      b_temp = b;
    end

    iter_avg_g(SNR_index) = iter_avg_g(SNR_index) + iter;

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

    wrong = (b ~= x');
    gwrong_where = gwrong_where + wrong;

    GAUSS = 0;	% GM or GAUSS

% if(0)
    %% direct ML

    ML = 1;

    % using Gauss-Markov model
    %
    if( size(avg_h,2) == 1)

      b = get_ML;

    else
    %
    %% using Gaussian model
    
      G=[];

      vit_C;
      b = get_ex;

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

    wrong = (b ~= x');
    Nwrong_where = Nwrong_where + wrong;

    ML = 0;
% end

 if(0)
    % run Kaleh algorithm
    %
    % initialize with average channel:
    %
    if( size(avg_h,2) == 1)
      FIR = avg_h/(1-alpha);	% Gauss-Markov
    else
      FIR = (sum(avg_h')/N)';	% Gaussian
    end

    KALEH = 1;
    G = [];

    if( isempty(PS))
      kaleh_init;       % one-time Kaleh initialization
    end

    for iter=1:4
      kaleh_iter;
    end

    % run Viterbi algorithm
    %
    vit_C;
    b = get_ex;

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

    wrong = (b ~= x');
    Kwrong_where = Kwrong_where + wrong;

    KALEH = 0;
  end

  end % trials
  
  % display results
  %
  disp( '*wrong_where:');

  [dwrong_where gwrong_where kwrong_where Nwrong_where] % Kwrong_where]

  dwrong = sum(dwrong_where);
  gwrong = sum(gwrong_where);
  kwrong = sum(kwrong_where);
  Nwrong = sum(Nwrong_where);

  % dwrong = sum(dwrong_where(1:N-FIR_length));
  % gwrong = sum(gwrong_where(1:N-FIR_length));
  % kwrong = sum(kwrong_where(1:N-FIR_length));

  % Kwrong = sum(Kwrong_where);   % for Kaleh

  % Kwrong = sum(Kwrong_where(1:N-FIR_length));   % for Kaleh

%  [dwrong gwrong kwrong]

  BER = [dwrong gwrong kwrong Nwrong] / (N*trials);

  % BER = [dwrong gwrong kwrong Kwrong] / ((N-FIR_length)*trials);        % for Kaleh
  
% print var_noise or SNR depending on choices above
%
%  [N var_noise BER]
%
  [N SNR BER]

end % for SNR or var_noise

disp( 'N, SNR, det, g, known, ML');

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
