% function [gn,En,bn,b_bar]=LC(p,SNRmfb,Ex_bar,Ntot,gap)
%
% Levin Campello's Method for a temporal channel
%
% Inputs
% p is the time-sampled pulse response
% SNRmfb is the SNRmfb, so EXbar*norm-p^2/sigma^2, in dB
% Ex_bar is the normalized energy
% Ntot is the total number of real subchannels, Ntot>2
%     The program user must adjust results for guard-extension losses
% gap is the gap in dB
%
% Outputs
% gn is channel gain
% En is the energy in the nth subchannel (PAM or QAM)
% bn is the bit in the nth subchannel (PAM or QAM)
% Nstar is the number of used subchannels
% b_bar is the bit rate per real dimension
%
% The first and last tones are PAM, the rest of them are QAM.
function [gn,En,bn,b_bar]=LC(p,SNRmfb,Ex_bar,Ntot,gap)
Noise_var=Ex_bar*(norm(p)^2)/(10^(SNRmfb/10));
gap=10^(gap/10);

% initialization
En=zeros(1,Ntot/2+1);
bn=zeros(1,Ntot/2+1);
gn=zeros(1,Ntot/2+1);
Hn = zeros(1,Ntot/2+1);
decision_table=zeros(1,Ntot/2+1);

% subchannel center frequencies
f=0:1/Ntot:1/2;

% find Hn vector
for i=1:length(p)
	Hn=Hn+p(i)*exp(j*2*pi*f*(i-1)); 
        % This value will be different depending if p represents 
        % p(1) + p(2)*D^-1 + ....  orp(1) + p(2)*D^+1....,
        % but we'll get same gn, thus same waterfilling result.
        % (Note that both have the same magnitude response!)
end

% find gn vector
gn=abs(Hn).^2/Noise_var;

%debugging purpose
%plot(gn)

%%%%%%%%%%%%%%%%%%%%%%%
% Now do LC %
%%%%%%%%%%%%%%%%%%%%%%%

%initialization

%used energy so far
E_so_far=0;
%decision table - QAM and PAM 
for n=2:Ntot/2
    if gn(n) ~= 0
        decision_table(n)=2*gap./gn(n);
    else
        decision_table(n)=inf;
    end
end
if gn(1) ~= 0
    decision_table(1)=3*gap/gn(1);
else
    decision_table(1)=inf;
end
if gn(Ntot/2+1) ~=0
    decision_table(Ntot/2+1)=3*gap/gn(Ntot/2+1);
else
    decision_table(Ntot/2+1)=inf;
end

%decision_table: debugging purpose

while(1)
    
    [y,index]=min(decision_table);
    E_so_far=E_so_far+y;
    
    if E_so_far > Ex_bar*Ntot
      
        break;
	
    else
      
        En(index)=En(index)+y;
        bn(index)=bn(index)+1;
	
	if (index ==1 | index == Ntot/2+1)
	  decision_table(index)=4*decision_table(index);
	else
	  decision_table(index)=2*decision_table(index);
	end

    end
    
end

% calculate b_bar
b_bar=1/Ntot*(sum(bn));

%>> [gn,En,bn,b_bar]=LC([1 0.9],10,1,8,0)
%
%gn =
%
%   19.9448   17.0320   10.0000    2.9680    0.0552
%
%
%En =
%
%    0.7521    1.7614    3.0000    2.0216         0
%
%
%bn =
%
%     2     4     4     2     0
%
%
%b_bar =
%
%    1.5000



    
        
        
