% function [b, GU, WU, S0, MSWMFU] = mu_mac(H, A, Lxu , cb)
%--------------------------------------------------------------------------
% Per-tonal (temporal dimension) multiuser mac receiver and per-user bits
%
% Inputs: H, A , Uind , cb
% Outputs: b, GU, WU, S0, MSWMFU
%
% definitions
%
% H: noise-whitened channel matrix [HU ... H1] Ly x sum-Lxu
% A: Block Diag sq-root sum-Lxu x sum-Lxu discrete modulators, 
%    blkdiag([AU ... A1]);  The Au entries derive from each MAC user's 
%    Lxu x Lxu input autocorrelation matrix, where the trace of each such 
%    autocorrelation matrix is user u's energy/symbol. This is per-tone.
%    
% Lxu: # of dimensions for each user U ... 1 in 1 x U row vector
% cb: = 1 if complex baseband or 2 if real baseband channel
%
% GU: unbiased feedback matrix sum-Lxu x sum-Lxu
% WU: unbiased feedforward linear equalizer sum-Lxu x sum-Lxu
% S0: sub-channel channel gains sum-Lxu x sum-Lxu
% MSWMFU: unbiased mean-squared whitened matched filter, sum-Lxu x Ly
% b - user u's bits/symbol 1 x U
%    the user should recompute b if there is a cyclic prefix
% 
%
%--------------------------------------------------------------------------
function [b, GU, WU, S0, MSWMFU] = mu_mac(H, A, Lxu, cb)

%
[~ , U] = size(Lxu);
b=zeros(1,U);

% Computing Ht: Ht = H*A
Ht = H*A;
%--------------------------------------------------------------------------

     
% Computing Rf, Rbinv, Gbar
Rf = Ht' * Ht;
Rbinv = Rf + eye(size(Rf));
Gbar = chol(Rbinv);

% Computing the matrices of interest
G = inv(diag(diag(Gbar)))*Gbar;
S0 = diag(diag(Gbar))*diag(diag(Gbar));
W = inv(S0)*inv(G');

GU = eye(size(G)) + S0*pinv(S0-eye(size(G)))*(G-eye(size(G)));
WU = pinv(S0-eye(size(G)))*inv(G');
MSWMFU = WU*Ht';
index=0;
for u=1:U
  for l=1:Lxu(u)
  b(u) = b(u)+ (1/cb)*log2(S0(index+l,index+l));
  end
  index=index+Lxu(u);
end

