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

%disp('initial FIR = [1 0 -1]''');
%disp('DET and GM init [0 1 0]')
%disp('DET and GM init using avg_h')
%disp('DET and GM init using ML b')
% disp('DET and GM init using FIR_save(:,1)')
% disp('Det init using avg_h')
% disp('Gauss init all 0''s');

disp('EM init using pilot, alpha and h_0')

% disp('EM init using ML results');

% disp('high correlation');
disp('uncorrelated h''s');

disp('GM h');

%disp('Gauss EM 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 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)

%
%  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.01
  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


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

% for constant SNR and time-varying noise
%
SNR_array = [0:3:15]

trials_array = 1 % [10 15 20 40 80 200] * 10

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

EM_local_min = zeros(length(SNR_array),1);
ML_local_min = zeros(length(SNR_array),1);

EM_equal_ML = zeros(length(SNR_array),1);

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

% for test with Gauss-Markov h and Gaussian model
%
% disp('switching to Gaussian model');
% 
% avg_h_save = avg_h;
% avg_h_g = avg_h * (ones(1,N)/(1-alpha))
% Rh_save = Rh;
% Rh_g = Rh/(1-alpha2)
% % inv_Rh_save = inv_Rh;
% % inv_Rh_g = inv(Rh_g)


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);
  % N2wrong_where = zeros(N,1);

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

  trial_count = 0;
  
  trials = trials_array(SNR_index);

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

    % if( mod(trial_count,100) == 0)
    %  trial_count
    % end

    % generate data
    %
    % for test with Gauss-Markov h and Gaussian model
    % 
	% GM_H = 1;
	% GAUSS_H = 0;
	% avg_h = avg_h_save;

    x4;

	% GM_H = 0;
	% GAUSS_H = 1;
	% avg_h = avg_h_g;

    % set pilot bits
    %
    % pilot(pilot_index) = x(pilot_index);
 
    % 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;

 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;

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

    if( GM_H)
      r = check_GM_ML_min( b);
    else
      r = check_Gauss_ML_min( b);
    end
   
    if( r)
      ML_local_min(SNR_index) = ML_local_min(SNR_index) + 1;
    end

%%% testing Gauss ML with GM h
%%
%
%%    avg_h = avg_h_g;
%%    Rh = Rh_g;
%%
%%    G=[];
%%    vit_C;
%%    b = get_ex;
%%
%%    avg_h = avg_h_save;
%%    Rh = Rh_save;
%%
%%    wrong = (b ~= x');
%%    N2wrong_where = N2wrong_where + wrong;
%
%%
%%%
    ML = 0;
 end

    %% run EM for deterministic h
    
 if(0)
    DET = 1;
    G=[];

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

    % init using channel initial conditions and pilot bits
    %
    % b = pilot;

    % init using ML results
    %
    % b = b_ML;

    %%% 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_C( b);
        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;
 end

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

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

    % init using channel initial conditions and pilot bits
    %
    % b = pilot;

    % init using ML results
    %
    % b = b_ML;

    %%% 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 = zeros(N,1);
    %
    % b = -x';

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

    b_temp = b;
    iter = 0;
    while(1)
      %%% get_Gauss_ML_cost(b)
      iter = iter + 1;
        B = get_A_C( b);
        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

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

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

    if( r)
      EM_local_min(SNR_index) = EM_local_min(SNR_index) + 1;
    % else
    %  disp('EM != local min');
    %  return;
    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;

    GM = 0;	% GM or GAUSS

 if(0)
    % run Kaleh algorithm
    %
    % initialize with average channel:
    %
    if( GM_H)
      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 ] % N2wrong_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));
  Nwrong = sum(Nwrong_where(1:N-FIR_length));
  % N2wrong = sum(N2wrong_where(1:N-2*FIR_length));

  % Kwrong = sum(Kwrong_where);   % for Kaleh

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

%  [dwrong gwrong kwrong]

  if( isempty(pilot))
    p = 0;
  else
    p = sum(pilot(1:N-FIR_length) ~= 0);
  end

  % BER = [dwrong gwrong kwrong Nwrong N2wrong] / ((N-FIR_length-p)*trials);

  % BER = [dwrong gwrong kwrong Nwrong] / ((N-2*FIR_length-p)*trials);

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

  % BER = [dwrong gwrong kwrong] / ((N-2*FIR_length)*trials);

  % BER = [dwrong gwrong kwrong Kwrong] / ((N-2*FIR_length-p)*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_array';
iter_avg_det
iter_max_det

iter_avg_g = iter_avg_g ./ trials_array';
iter_avg_g
iter_max_g

% fprintf( 1, 'EM_local_min %% = %g\n', 100*EM_local_min/(trials*length(SNR_array)) );

disp('EM local min %');

100*EM_local_min./trials_array'

disp('ML local min %');

100*ML_local_min./trials_array'

disp('EM_equal_ML %');

100*EM_equal_ML./trials_array'
