% function [snrGDFEu, GU, WU, S0, MSWMFU, b, bbar, snrGLEu] = computeGDFE(H, A, cb, Lx)
%--------------------------------------------------------------------------
% Inputs
% H: noise-whitened channel, Ly x Lx (one tone)
% A: any Lx x Lx square root of input autocorrelation matrix
%        This can be generalized non-square square-root
% cb: =1 if H is complex baseband; =2 if H is real baseband
% Lx: optional input of Lx when not equal to size of A
%     this is used to compute bits/dimension properly
%
% Outputs
% GU: unbiased feedback matrix
% WU: unbiased feedforward linear equalizer
% S0: sub-channel channel gains
% MSWMFU: unbiased mean-squared whitened matched filter
% b: bit distribution vector
% bbar: number of bits/dimension (real if cb=2, complex if cb=1)
% snrGDFEu - unbiased SNR in dB; assumes size of R_sqrt input
%    the user should recompute SNR if there is a cyclic prefix
% b = bit distribution over symbol dimensions (real cb=2; complex cb=1)
% bbar is sum(b)/Lx so total bits/(real cb=2; cplx cb=1) dimension
% snrGLEu is the linear GLE SNR
%
% Note: The R_sqrt need not be square non-singular, as long as
% R_sqrt*R_sqrt' = input autocorrelation matrix Rxx, but SNRLEu is off
%Thanks to Ethan Liang, corrected/updated J. Cioffi
%
%--------------------------------------------------------------------------
function [snrGDFEu, GU, WU, S0, MSWMFU, b , bbar, snrGLEu] = computeGDFE(H, A, cb, Lx)
%-------------------------------------------------------------------------
% Computing Ht: Ht = H*A
Ht = H*A;
%--------------------------------------------------------------------------

if (nargin<4)
  [~ , Lx , ~] = size(A);
end
[Ly,~,~]=size(H);
     
% Computing Rf, Rbinv, Gbar
Rf = Ht' * Ht;
Rbinv = Rf + eye(size(Rf));
Gbar = chol(Rbinv);
Ng=size(Gbar);
Ng=Ng(1);

% 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';
% Section to handle numerical is with large-dim determinants/geo-ave
if Ng < 10
   snrGDFEu = 10*log10(det(S0)^(1/(Lx))-1);
else
   snrGDFE=0;
   for n=1:Ng
    snrGDFE = snrGDFE + (10/(Lx))*log10(S0(n,n));
   end
   snrGDFEu=10*log10(10^(0.1*snrGDFE)-1);
end
b=(1/cb)*log2(diag(S0));
bbar=(1/Lx)*sum(b);
% This is linear equalizer computation
sGLE0=inv(diag(diag(inv(Rbinv))));
snrGLEu=10*log10(det(sGLE0)^(1/(Lx))-1);




