function [gn,En,bn,b_bar,SNRdmt]=DMTLCra(P,SNRrcvr,Ex_bar,N,gap)

%
% EE379C 2008 Spring
%
% Levin Campello's Method with DMT
%
% Inputs
% P is the pulse response(an nu is set to length(p) - 1
% SNRrcvr is the channel output SNR for an i.i.d. input
% Ex_bar is the normalized energy
% N is the total number of real/complex subchannels, N>2
% gap is the gap in dB
%
% Outputs
% gn is the vector of channel gains
% En is the vector energy distribution (PAM or QAM)
% bn is the vector bit distribution (PAM or QAM)
% b_bar is the number of bits per dimension in the DMT symbol
%
% The first and last bins are PAM; the rest are QAM.
  
% dB into normal scale
Noise_var=Ex_bar*(norm(P)^2)/(10^(SNRrcvr/10));
gap=10^(gap/10);
nu=length(P)-1;

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

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

% find Hn vector
for i=1:length(P)
	Hn=Hn+P(i)*exp(j*2*pi*f*(i-1)); 
end

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

%debugging purpose
%plot(gn)

%%%%%%%%%%%%%%%%%%%%%%%
% Levin Campello Loading %
%%%%%%%%%%%%%%%%%%%%%%%

%initialization

%used energy so far
E_so_far=0;
%decision table - QAM and PAM 
decision_table(2:N/2)=2*gap./gn(2:N/2);
% Gap formula incremental energies.
if gn(1) ~= 0
    decision_table(1)=3*gap/gn(1);
else
    decision_table(1)=inf;
end
if gn(N/2+1) ~=0
    decision_table(N/2+1)=3*gap/gn(N/2+1);
else
    decision_table(N/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*N
      
        break;
	
    else
      
        En(index)=En(index)+y;
        bn(index)=bn(index)+1;
	
	if (index ==1 | index == N/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/(N+nu)*(sum(bn));
SNRdmt=10*log10(gap*(2^(2*b_bar)-1));