% function [Bu, GU, S0, MSWMFunb , B] = mu_bc(H, AU, Lyu , cb)
%--------------------------------------------------------------------------
% Inputs: Hu, AU , Usize, cb
% Outputs: Bu, Gunb, Wunb, S0, MSWMFunb
%
% H: noise-whitened BC matrix [H1 ; ... ; HU] (with actual noise, not wcn)
%    sum-Ly x Lx x N
% AU: Block-row square-root discrete modulators, [A1 ... AU]
%     Lx  x (U * Lx) x N
% Lyu: # of (output, Lyu) 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 precoder matrices: (Lx U) x (Lx U) x N
%     For each of U users, this is Lx x Lx matrix on each tone
% S0: sub-channel dimensional channel SNRs: (Lx U) x (Lx U) x N
% MSWMFunb: users' unbiased diagonal mean-squared whitened matched matrices
%    For each of U cells and Ntones, this is an Lx x Lyu  matrix
% Bu - users bits/symbol 1 x U
%    the user should recompute SNR if there is a cyclic prefix
% B - the user bit distributions (U x N) in cell array
% 
%
%--------------------------------------------------------------------------
function [Bu, Gunb, S0, MSWMFunb, B] = mu_bc(H, AU, Lyu, cb)

%
[~,U] = size(Lyu);
Bu=zeros(1,U);
[Ly,Lx,N]=size(H);

% Computing Ht: Ht = H*A

%--------------------------------------------------------------------------
S0=cell(U,N);
% b=cell(U,N);
for u=1:U
    for n=1:N
     S0{u,n} = zeros(Lx,Lx);
%     b{u,n}=0;
    end
end
MSWMFunb=cell(U,N);
Gunb=cell(U,N);

for u=1:U
    for n=1:N
      MSWMFunb{u,n}=zeros(Lx,Lyu(u));
      Gunb{u,n}=zeros(Lx,Lx*sum(Lyu)); 
    end
end
     
% Use compute GDFE

for u=1:U
  for n=1:N
   [~,Gunbtemp(:,:,n),~,S0temp(:,:,n),MSWMFtemp(:,:,n),~,~,~]= computeGDFE( ...
       H(sum(Lyu(1:u))-Lyu(u)+1:sum(Lyu(1:u)),:,n), AU(:,:,n),cb);
   Gunb{u,n} = Gunbtemp(Lx*(sum(Lyu(1:u))-Lyu(u))+1:Lx*(sum(Lyu(1:u))),:,n);
   S0{u,n}=S0temp(Lx*(sum(Lyu(1:u))-Lyu(u))+1:Lx*(sum(Lyu(1:u))),...
       Lx*(sum(Lyu(1:u))-Lyu(u))+1:Lx*(sum(Lyu(1:u))),n);
   MSWMFunb{u,n}=MSWMFtemp(Lx*(sum(Lyu(1:u))-Lyu(u))+1:Lx*(sum(Lyu(1:u))),:,n);
   B{u,n}= (1/cb)*log2(det(S0{u,n}));
  end
   Bu(u)=sum([B{u,:}],2);
end







