% function [b, GU, WU, S0, MSWMFU] = mu_ic(H, A, Lxu, Lyu, cb)
%--------------------------------------------------------------------------
% Per-tonal (temporal dimension) multiuser interference channel
% receiver and per-user bits - Chris Adelico Ferrarin - 2023
%
% Inputs: H, A, Lxy, Lyu, cb
% Outputs: b, GU, WU, S0, MSWMFU
%
% Definitions:
%
% H: noise-whitened channel matrix [HUU ⋯ HU1] sum-Lyu x sum-Lxu
%                                  | ⋮   ⋱  ⋮ |
%                                  |H2U ⋯ H21|
%                                  [H1U ⋯ H11]
% A: Block Diag sq-root sum-Lxu x sum-Lxu discrete modulators, 
%    blkdiag([AU ... A1]);  The Au entries derive from each IC 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 input dimensions for each user U ... 1 in 1 x U row vector
% Lyu: # of output 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 x U with matrices indexed
%     from user U (e.g. GU(:,:,2) gives GU for user U-1).
% WU: unbiased feedforward linear equalizer sum-Lxu x sum-Lxu x U with
%     matrices indexed from user U (e.g. WU(:,:,2) gives WU for user U-1).
% S0: sub-channel channel gains sum-Lxu x sum-Lxu x U with matrices indexed
%     from user U (e.g. GU(:,:,2) gives S0 for user U-1).
% MSWMFU: unbiased mean-squared whitened matched filter, sum-Lxu x Ly x U
%         with matrices indexed from user U (e.g. MSWMFU(:,:,2) gives
%         MSWMFU for user U-1).
% 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_ic(H, A, Lxu, Lyu, cb)

% Initializing variables
[~ , U] = size(Lxu);
b_mat = zeros(U,U);
GU = [];
WU = [];
S0 = [];
MSWMFU = [];
endrow = sum(Lyu);

for u = 1:U
    % Extracting the channel matrix for this user's receiver
    lyu = Lyu(end - u + 1);
    startrow = endrow - lyu + 1;
    Hu = H(startrow:endrow, :);

    % Finding the MMSE solution for the MAC formed by this user's receiver
    [b_mac, GU_mac, WU_mac, S0_mac, MSWMFU_mac] = mu_mac(Hu, A, Lxu, cb);

    % Updating array of IC filter solutions
    b_mat(u,:) = b_mac;
    GU = cat(3, GU_mac, GU);
    WU = cat(3, WU_mac, WU);
    S0 = cat(3, S0_mac, S0);
    MSWMFU = cat(3, MSWMFU_mac, MSWMFU);

    endrow = endrow - lyu;
end

% Finding the minimum bit rate for each user across all MAC solutions
b = min(b_mat, [], 1);
