trials = 50000

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

disp('using uncorrelated h''s');

% disp('skipping small norm FIR''s');

% disp('using Ghosh and Webers fixed FIR from Fig. 2')

% disp('using poor initialization');

% initialize rand() and randn() with time-of-day
%
  seed=floor(sum(100*clock))
%
%  disp('seed from ml.01')
%  seed = 210886

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

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 Rh;
global inv_Rh;
global inv_Rh_d;
global avg_h;
global FIR_save;
global FIR_length;

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

% 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 var_noise;	% noise variance
global gamma_scale;     % normalization for columns of gamma and beta

% set FIR length and PAM
%
  FIR_length = 3

  PAM = 2

% for Gaussian FIR
%
  avg_h = ones(FIR_length,1)
  var_h = 0.05
  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
%
%  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))/10;
%%  inv_T(:,1) = 0.9;
%%  inv_T(1,1) = 1
%%
%%  T = inv(inv_T)
%%
%%  inv_Rh = (T'*T)/var_h

Rh = inv(inv_Rh)

% set N and SNR
%
  N = 10 + FIR_length

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

% FIR_skip = 0;

for SNR_index = 1:length(SNR_array)

  % run the Monte-Carlo simulations
  %
  SNR = SNR_array(SNR_index);

  % run the Monte-Carlo simulations
  %
  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;
  
    % set FIR
    %
    % uniform FIR over [0,2]
    %
    % FIR = 2*rand(FIR_length,1);
    %
    % Gaussian FIR with mean 1
    %
     FIR = avg_h + inv_T * (stdev_h * randn(FIR_length,1));
     %
     % skip FIR's with small norm, 31 May 2000
     %
     % s = sum(FIR.*FIR);
     % while( s < 1)
     %   FIR_skip = FIR_skip + 1;
     %   FIR = avg_h + inv_T * (stdev_h * randn(FIR_length,1));
     %   s = sum(FIR.*FIR);
     % end
 
    % force to have zero at DC
    %
    % FIR(FIR_length) = -sum(FIR(1:FIR_length-1));

    % try Ghosh and Weber's fixed FIR
    %
    % FIR = [0.407 0.815 0.407]';

    % generate data
    %
    x4;

    % known channel case
    %
    KNOWN = 1;
    G = [];
    vit_C;
    b = get_ex;
  
    b_save = b;

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

  if(0)
    % get ML solution
    %
    ML = 1;

    b = get_ML;

    b_ML = b;

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

    ML = 0;
  end

if(0)
    % run EM for deterministic h
    %
    % use true symbols
    %
    %  b=x';
    %
    % or just use b from the known channel case
    % 
    % or reinitialize:
    %
    %  x4_em_init_det;
    %
    DET = 1;
    G = [];

    % init using b_ML
    %
    % b = b_ML;

    % reasonable init:
    %
    FIR = avg_h;
    trans = FSM_trans( states, FIR);
    vit_C;
    b = get_ex;

    % Ghosh and Weber init:
    %
    % [0 ... 0 1 0 ... 0]
    %
    % FIR = zeros(FIR_length,1);
    % FIR( (size(FIR,1)+1)/2) = 1;
    % trans = FSM_trans( states, FIR);
    % vit_C;
    % b = get_ex;

    % 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);
	FIR = pinv(B)*(z');
	trans = FSM_trans( states, FIR);
	vit_C;
	b = get_ex;
      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 Gaussian h
    %
    % use true symbols
    %
    %  b=x';
    %  G = 1;
    %
    % or just use b from the known channel case:
    %
    % or reinitialize:
    %
    %  x4_em_init_gauss;
    %
    GAUSS = 1;

    % init using b_ML
    %
    % b = b_ML;

    % reasonable init:
    %
    FIR = avg_h;
    trans = FSM_trans( states, FIR);
    G = Rh;
    vit_C;
    b = get_ex;

    % bad init:
    %
    % b = b_bad;

    b_temp = b;
    iter = 0;
    while(1)
      iter = iter + 1;
        B = get_A( b, N, Memory);
        G = inv( (B'*B)/var_noise + inv_Rh);
        FIR = G * (B'*(z')/var_noise + inv_Rh*avg_h);
        trans = FSM_trans( states, FIR);
        vit_C;
        b = get_ex;
      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;
end

  if(0) % skip Kaleh
    % run Kaleh algorithm
    %
    % with bad init:
    %
    %  FIR = avg_h;
    %
    % initialize with true channel:
    %
     FIR = FIR_save;

    G = [];

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

    for iter=1:4
      kaleh_iter;
    end

    % run Viterbi algorithm
    %
    trans = FSM_trans( states, FIR);
    vit_C;
    b = get_ex;

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

  end % skip Kaleh

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

  [dwrong_where gwrong_where kwrong_where Nwrong_where]

%  [dwrong_where gwrong_where kwrong_where Kwrong_where]	% for Kaleh

%  dwrong = sum(dwrong_where);
%  gwrong = sum(gwrong_where);
%  kwrong = sum(kwrong_where);
%  Nwrong = sum(Nwrong_where);
%
%  Kwrong = sum(Kwrong_where);	% for Kaleh

  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));
%
%  Kwrong = sum(Kwrong_where(1:N-FIR_length));	% for Kaleh
%
%  [dwrong gwrong kwrong]

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

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

%  BER = [dwrong gwrong kwrong Kwrong] / (N*trials);	% for Kaleh
%
%  BER = [dwrong gwrong kwrong Kwrong] / ((N-FIR_length)*trials);	% for Kaleh
  
  [N SNR BER]

end % for SNR

disp( 'N, SNR, det, gauss, 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

% FIR_skip
