% check EM convergence vs. number of iterations

trials = 20

N_ITER = 5

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

% set FIR length and PAM
%
  FIR_length = 3

  PAM = 2

% 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

% set N and SNR
%
  N = 20

  SNR = 10

  % run the Monte-Carlo simulations
  %
  dnorm = zeros(N_ITER,1);
  gnorm = zeros(N_ITER,1);

  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));
  
    % generate data
    %
    x4;
  
    G=[];

    % run EM for deterministic h
    %
    x4_em_init_det;
  
    for iter = 1:N_ITER
      x4_em_iter;
      f = FIR - FIR_save;
      dnorm(iter) = dnorm(iter) + sum(f.*f);
    end

    % run EM for Gaussian h
    %
    x4_em_init_gauss;
  
    for iter = 1:N_ITER
      x4_em_iter;
      f = FIR - FIR_save;
      gnorm(iter) = gnorm(iter) + sum(f.*f);
    end

  end
  
  [ [1:N_ITER]' dnorm/trials gnorm/trials]
