% with linear constraints

trials = 100

disp('using N bits to compute BER');
disp('using 2 EM iterations for all algorithms');

global B;
global var_noise;
global inv_Rh;
global inv_Rh_d;
global avg_h;
global FIR_save;
global FIR_length;
global FIR_length_p;	% underlying channel coefficients

global FIR_M;
global FIR_M_t;	% FIR_M'
global FIR_M_p_;

global DD;
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.8
alpha2 = alpha*alpha;


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

% set FIR length
%
  FIR_length = 4;

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

  FIR_length_p = FIR_length - size(FIR_F,1);

% 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_t = FIR_M';

  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 or the Gaussian part of Gauss-Markov
%
% references to _h here really mean _p for the constrained case
%
% Gauss-Markov
%
  avg_h = ones(FIR_length_p,1)*(1-alpha)	% so avg_h/(1-alpha) = all 1's
%
% Gaussian with time-varying mean 1 -> 2
%
%  avg_h = ones(FIR_length_p,1) * (1 + [0:N-1]*((2-1)/(N-1)))

  var_h = 0.1
  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)
%
  T = tril(ones(FIR_length_p,FIR_length_p))

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

  if( size(avg_h,2) == 1)
    inv_Rh_d = inv_Rh*avg_h;	% inv(C)*d in the paper
  end


% for constant noise variance and time-varying SNR
%
% for var_noise = 1
%
% for constant SNR and time-varying noise
%
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 FIR
    %
    FIR_save = zeros(FIR_length,N);

    if( size(avg_h,2) == 1) % Gauss-Markov
      %
      % initial FIR equal to the steady-state average
      %
        FIR_p = avg_h/(1-alpha);
        FIR_save(:,1) = FIR_M * FIR_p + FIR_M_p_;

      % succeeding FIR's have Markov statistics
      %
        for i=2:N
          FIR_p = alpha*FIR_p + avg_h + inv_T * (stdev_h * randn(FIR_length_p,1));
          FIR_save(:,i) = FIR_M * FIR_p + FIR_M_p_;
        end

    else    % Time-varying Gaussian

      for i=1:N
        FIR_p = avg_h(:,i) + inv_T * (stdev_h * randn(FIR_length_p,1));
        FIR_save(:,i) = FIR_M * FIR_p + FIR_M_p_;
      end

    end

    % generate data
    %
    x4;
  
    % known channel case
    %
    G=[];
  
    x4_em_init;

    b_save = b;

    wrong = (b ~= x');
    kwrong_where = kwrong_where + wrong;
  
    % run EM for unconstrained deterministic h
    %
    for iter=1:2
      x4_em_iter;
    end
  
    wrong = (b ~= x');
    dwrong_where = dwrong_where + wrong;
  
    % run EM for time-varying Gaussian or Gauss-Markov h
    %
    b = b_save;
  
    G=1;

    for iter=1:2
      x4_em_iter;
    end
  
    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);

  BER = [dwrong gwrong kwrong] / (N*trials);
  
% print var_noise or SNR depending on choices above
%
%  [N var_noise BER]
%
  [N SNR BER]

end % for SNR or var_noise
