% using linear constraints on the FIR

trials = 20000

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

disp('using 2 EM iterations for all algorithms');

% 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 = 4

  PAM = 2

% set N
%
  N = 12;

% set FIR constraints
%
% double zero at DC:
%
  FIR_F=[1 1 1 1; 0 1 2 3]
  FIR_f=[0; 0]

% set orthogonal decomposition
%
  [FIR_U FIR_S FIR_V] = svd(FIR_F);

  FIR_rank = size(FIR_F,1); % assume full rank

  FIR_M = FIR_V(:,FIR_rank+1:FIR_length);

  FIR_M_ = FIR_V(:,1:FIR_rank);

  FIR_S1 = FIR_S(1:FIR_rank,1:FIR_rank);

  FIR_p_ = inv(FIR_S1)*(FIR_U'*FIR_f);

  FIR_M_p_ = FIR_M_ * FIR_p_;

% for Gaussian FIR
%
% references to _h here really mean _p for the constrained case
%
% FIR_length-2 is used below for the double-zero-at-DC constraint
%
  avg_h = ones(FIR_length-2,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-2,FIR_length-2))

% set inv_T and inv_Rh
%
  inv_T = inv(T);
  inv_Rh = (T'*T)/var_h

for SNR = 0:2:10

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

    % set Gaussian FIR
    %
    FIR_p = avg_h + inv_T * (stdev_h * randn(FIR_length-2,1));
    FIR = FIR_M * FIR_p + FIR_M_p_;

    % generate data
    %
    x4;

    % known channel case
    %
    x4_em_init_known;
  
    b_save = b;

    wrong = (b ~= x');
    kwrong_where = kwrong_where + wrong;
  
    % run EM for deterministic h
    %
    % use b from the known channel case

    for iter=1:2
      x4_em_iter;
    end

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

    % run EM for Gaussian h
    %
    % use b from the known channel case:
    %
      b = b_save;
      G = 1;

    for iter=1:2
      x4_em_c_iter;
    end

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

  end
  
  % display results
  %
  disp( 'det, Gauss, known');
  disp( 'wrong_where:');

  [dwrong_where gwrong_where kwrong_where]

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

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

end % for SNR
