% Gauss-Markov multi-user ISI problem
%
% call xgm1_init first
%
% R. Perry, Oct. 2001

trials = 300

for SNR_index = 1:size(SNR,2)

  disp( 'SNR');

  SNR(:,SNR_index)

  BER_known = zeros( Nusers, 1);
  BER_gm = zeros( Nusers, 1);

  for trial_count = 1:trials

    % generate each user's sequence
    %
    x = 2 * (rand( Nsyms, Nusers) > 0.5) - 1;

    % generate the FIR filters
    %
    FIR_0 = (stdev_v/sqrt(1-alpha2)) * randn( FIR_length, Nusers, Nchips);

    FIR_save = zeros( FIR_length, Nusers, Nchips, Nsyms);

    F = FIR_0;

    for t = 1:Nsyms
      F = alpha * F + stdev_v * randn( FIR_length, Nusers, Nchips);
      FIR_save(:,:,:,t) = F;
    end

    % normalize the FIR filters
    %
    F1 = FIR_save(:,1,:,:);
    F1_norm = sum(F1(:).*F1(:));
    for u = 2:Nusers
      F = FIR_save(:,u,:,:);
      F_norm = sum(F(:).*F(:));
      scale = sqrt( (F1_norm/F_norm) * 10^((SNR(u,SNR_index)-SNR(1,SNR_index))/10) );
      FIR_save(:,u,:,:) = F * scale;
      FIR_0(:,u,:) = FIR_0(:,u,:) * scale;
    end

    % apply the FIR filters
    %
    y = zeros( Nchips, Nsyms);

    for u = 1:Nusers
      xu = [flipud(states(1:Memory,1)); x(:,u)]'; % initial -1's in the channels
      for t = 1:Nsyms
        xt = xu(t+Memory:-1:t);
        for j = 1:Nchips
          y(j,t) = y(j,t) + xt * FIR_save(:,u,j,t);
        end
      end
    end

    % noise variance based on user #1
    %
    var_noise = 10^(-SNR(1,SNR_index)/10) * F1_norm/(Nsyms*Nchips);

    noise = sqrt(var_noise).*randn( Nchips, Nsyms);

    z = y + noise; % signal plus noise

    %%%
    %%
    % run Viterbi
    %
    if( one_trellis)

      Nv = 1;

      %%
      %
      TYPE = 30;       % multi-user, Gauss-Markov, known channel

      L = 1;

      vit_init( L, M, Nsyms, Nv);

      for t=1:Nsyms
        vit_update;
      end

      paths = vit_paths;

      ex_known = zeros( Nsyms, L, Nusers);

      for u=1:Nusers
        ex_known(:,:,u) = get_ex_multi( paths, u);
      end

      tmp = (x ~= reshape(ex_known(:,1,:),size(x)));
      BER_known = BER_known + sum(tmp)';

      %%
      %
      TYPE = 31;       % multi-user, Gauss-Markov cost

      L = 3; % get_Lmax;

      vit_init( L, M, Nsyms, Nv);

      for t=1:Nsyms
        vit_update;
      end

      paths = vit_paths;

      ex_gm = zeros( Nsyms, L, Nusers);

      for u=1:Nusers
        ex_gm(:,:,u) = get_ex_multi( paths, u);
      end

      tmp = (x ~= reshape(ex_gm(:,1,:),size(x)));
      BER_gm = BER_gm + sum(tmp)';

    %%
    %
    else % separate trellises

      Nv = Nusers;

      %%
      %
      TYPE = 40;       % multi-user, Gauss-Markov, known channel, separate trellises

      L = 3

      ex_known = zeros( Nsyms, L, Nusers);

      ex_p = zeros( L, Nusers);

      vit_init( L, M, Nsyms, Nv);

      Piter = 1

      for t=1:Nsyms
        Vw_save = Vw;
        for P = 1:Piter % iterate
          for u = 1:Nusers
            Vparam( 4 + u) = t - 1;     % reset Viterbi t
            Vw(:,:,u) = Vw_save(:,:,u); %  and cumulative costs
            vit_update( u);
            paths = vit_paths( u);
            ex_known(1:t,:,u) = get_ex( paths);
            p = exp( -(VW - min(VW)));   % change costs to probabilities
            ex_p(:,u) = p/sum(p);        %  and normalize so sum is 1.
          end
        end
      end

      tmp = (x ~= reshape(ex_known(:,1,:),size(x)));
      BER_known = BER_known + sum(tmp)';

    end

  end % trials

  BER_known = BER_known / (trials*Nsyms)

  BER_gm = BER_gm / (trials*Nsyms)

end % SNR
