function [b, E, g, Marg] =dmt(Hds,Nds,E_total,Gap,PSDmask,const,usedtones);
%% ===========================================================================
%dmt - Computes bit distribution, and energy distribution for a dmt system.
% Parameter:    Hds   		The channel magnitude (one-sided, mW), where
%                               upstream and downstream bandwidth
%                               allocation must have been taken into  account. 
% Parameter:    Nds	     	The noise psd (mW)
% Parameter:    E_total	        The total energy (mW)
% Parameter:    Gap	        The effective Gap (mW)
% Parameter:    PSDmask	        The nominal transmit PSD (mW)
% Parameter:    const	        The constellation sizes used
% Parameter:    usedtones	The dmt tones used
%
% Returns:      b    	        Sub-channel bit rate's
% Returns:      E    	        Sub-channel Energy
% Returns:      g               Sub-channel normalized SNR's
% Returns:      Marg            Sub-channel margin's
% Example(s):
%
% Algorithmic details:
%
% Reference:
%    This code is based on code from Atul Salvekar @ Standford University
%% ===========================================================================

%% ===========================================================================
% Copyright (C):                                        
%       1999-2000 by Telia Research AB, Lulea, Sweden;                
%       2001-2002 by Forschungszentrum Telekommunikation Wien, Austria;
%                                                         All rights reserved.
% Project       : xDSL simulation tool
% Author(s)     : Daniel Bengtsson (Daniel.J.Bengtsson@Telia.se)
%               : Tomas Nordstrom (Tomas.Nordstrom@FTW.at)
%
% CVS:       $Id: dmt.m,v 1.4 2002/02/11 16:17:53 tono Exp $
%% ===========================================================================
% Change History
%      1999-xx-xx (DaB)  Import from a dmt bit-loading code from Stanford
%      2000-xx-xx (DaB)  Added margin calculations
%      2000-03-30 (ToNo) Moved around the file
%      2001-08-24 (ToNo) Added header and comments
%      2001-08-24 (ToNo) Fixed so that only used tones are part of
%                        the bit-loading
%      2002-01-29 (ToNo) Noted the incorrect margin calculation
%      2002-02-11 (ToNo) Fixed the case of no possible carriers at all
%% ===========================================================================

Ntotal=length(Hds);                  % Total number of sub-channels

if nargin<7,                         % If no "usedtones" given, use all
    usedtones=1:Ntotal;
end

b_min=min(const);                    % Limit constellation sizes
b_max=max(const);

% Sub-channel normalized SNR's // Clear out all non-used channels!
g=1e-14*ones(size(Hds));
g(usedtones)=(Hds(usedtones).^2 ./ Nds(usedtones)); 

[g_sorted,index] = sort(g);          % Sort SNR 
g_sorted=fliplr(g_sorted);           % end up with largest gain at 1
index=fliplr(index);

g=(Hds.^2 ./ Nds);                   % Sub-channel normalized SNR's,
                                     % that is, g= SNR./PSDmask

% Start from worst channel (with non-zero gain)
num_zero_g=length(find(g_sorted<10^-8));     
Nstar=Ntotal-num_zero_g;        

%% ===========================================================================
% Water-filling
if Nstar>0
while(1)
    K = 1/Nstar*(E_total+Gap*sum(1./g_sorted(1:Nstar)));
    E_min=K-Gap/g_sorted(Nstar);     % E_min occurs in the worst channel
    if (E_min<0)
        Nstar=Nstar-1;               % If negative reduce the no. of channels
    else
        break;                       % If all is positive, we are done
    end;
end;

Es=K-Gap./g_sorted(1:Nstar);         % Energies of subchannels 

count = 0;
E_extra = 0;
for k=1:Nstar                        % Sum up all excess energy
    if (Es(k) < PSDmask(index(k)))
      count = count + 1;
    else
      E_extra = E_extra + Es(k) - PSDmask(index(k));
    end;
end;
if (count == 0)		% count==0 means E_extra=0. avoid 0/0
    count = 1;
end;

E=zeros(1,Ntotal);
for k=1:Nstar		             % Redistribute excess energy
    if (Es(k)+E_extra/count < PSDmask(index(k)))
      E(index(k)) = Es(k) + E_extra/count;
    else
      E(index(k)) = PSDmask(index(k));
    end;
end;

%% ===========================================================================
% Bit rates of subchannels
b = log2(1+E.*g/Gap);             

% Make bit-rates integer and within limits
b=round(b);       
b(find(b>b_max))=b_max*ones(size(find(b>b_max)));  
b(find(b<b_min))=zeros(size(find(b<b_min)));

% Calculate the energy for each channel
E=(2.^b-1).*Gap./g;

% Calculate the signal margin for each sub-channel
tmp=find(b>0);
Marg=zeros(size(b));
Marg(tmp) = E(tmp).*g(tmp)./(2.^b(tmp)-1)./Gap;
% CHECK This seems not do the right thing
% It only calculates 1 for each channel where the b>0 ???????

else
    b = zeros(1,length(Hds));
    E = zeros(1,length(Hds));
    Marg = zeros(1,length(Hds));
end