trials = 1

global B;
global var_noise;
global inv_Rh;
global inv_Rh_d;
global FIR_save;

% initialize rand() and randn() with time-of-day
%
  seed=floor(sum(100*clock))
  rand( 'state', seed);
  randn( 'state', seed);

% set FIR length
%
  FIR_length = 3;

% for Gaussian FIR
%
  avg_h = ones(FIR_length,1)
  var_h = 0.5
  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_Rh_d = inv_Rh * avg_h;

% set N and SNR
%
  N = 6;

for SNR = 6

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

  trial_count = 0;
  
  while( trial_count < trials)
  
    trial_count = trial_count + 1;
  
    % generate data
    %
    x4;
  
    % run EM for deterministic h
    %
    G=[];
  
    x4_em_init;
  
    wrong = (b ~= x');
    kwrong_where = kwrong_where + wrong;
  
    x4_em_iter;
    x4_em_iter;
  
    wrong = (b ~= x');
    dwrong_where = dwrong_where + wrong;
  
    % run EM for Gaussian h
    %
    G=1;
  
    x4_em_init;
  
    x4_em_iter;
    x4_em_iter;
  
    wrong = (b ~= x');
    gwrong_where = gwrong_where + wrong;
  end
  
  % display results
  %
  disp( '*wrong_where:');

  [dwrong_where gwrong_where kwrong_where]

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

%  [dwrong gwrong kwrong]

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

end % for SNR
