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

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

% set FIR
%
  FIR_p = [ 1; 1]

  FIR_save = FIR_M * FIR_p + FIR_M_p_
  
for SNR = 0:2:10

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

  trial_count = 0;
  
  while( trial_count < trials)
  
    trial_count = trial_count + 1;
  
    % generate data
    %
    FIR = FIR_save;
    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 constrained deterministic h
    %
    % use b from the known channel case

    b = b_save;

    for iter=1:2
      x4_em_c_iter;
    end

    wrong = (b ~= x');
    dcwrong_where = dcwrong_where + wrong;
  
  end
  
  % display results
  %
  disp( 'det, det_c, known');
  disp( 'wrong_where:');

  [dwrong_where dcwrong_where kwrong_where]

%  dwrong = sum(dwrong_where(1:N-FIR_length));
%  dcwrong = sum(dcwrong_where(1:N-FIR_length));
%  kwrong = sum(kwrong_where(1:N-FIR_length));
%
%  BER = [dwrong dcwrong kwrong] / ((N-FIR_length)*trials);

  dwrong = sum(dwrong_where);
  dcwrong = sum(dcwrong_where);
  kwrong = sum(kwrong_where);

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

end % for SNR
