function [gn,en_bar,bn_bar,Nstar,b_bar_check,margin]=DMTma(P,SNRrcvr,Ex_bar,b_bar,N,gap)

% EE379C 2008 Spring,
% J. Cioffi - thanks to Seong Taek Chung
%
% P is the pulse response where v=length(P) -1
% SNRrcvr is the SNR (in dB) at te receiver input if the input is i.i.d
% b_bar is the target normalized bit rate
% Ex_bar is the normalized energy
% N is the total number of real/complex subchannels, N>2
% gap is the gap in dB
%
% gn is the vector of channel gains
% en_bar is the vector of energy/dim for all subchannels
% bn_bar is the vector of bit/dim for all subchannels
% Nstar is the number of used subchannel

% 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);
bn=zeros(1,N);
gn=zeros(1,N);
Hn = zeros(1,N);

% subchannel center frequencies
f=-1/2+1/N: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;

%%%%%%%%%%%%%%%%%%%%%%%
% MA waterfilling %
%%%%%%%%%%%%%%%%%%%%%%%


%sort
[gn_sorted, Index]=sort(gn);  % sort gain, and get Index

gn_sorted = fliplr(gn_sorted);% flip left/right to get the largest 
                              % gain in leftside
Index = fliplr(Index);        % also flip index  

num_zero_gn = length(find(gn_sorted == 0)); %number of zero gain subchannels
Nstar=N - num_zero_gn;    
 	% Number of used channels, 
 	% start from N - (number of zero gain subchannels)
b=(N+nu)*b_bar;
Klog2tilde=2*b-(1/log(2))*sum(log(gn_sorted));

while(1) 
 	Klog2=(1/Nstar)*Klog2tilde + log(gap)/log(2);
 	En_min=2^(Klog2)-gap/gn_sorted(Nstar);	% En_min occurs in the worst channel
 	if (En_min<0)
            Klog2tilde=Klog2tilde+(1/log(2))*log(gn_sorted(Nstar));
    		Nstar=Nstar-1;  % If negative En, continue with less channels
 	else 
    		break;       % If all En positive, done.
 	end
end
K=2^(Klog2);
En=K-gap./gn_sorted(1:Nstar); 		% Calculate En
Bn=.5*log2(K*gn_sorted(1:Nstar)/gap); 	% Calculate bn

bn(Index(1:Nstar))=Bn;		% return values in original index
en(Index(1:Nstar))=En;		% return values in original index

en_bar=en;
bn_bar=bn;

% calculate b_bar
b_bar_check=1/(N+nu)*(sum(bn));

%calculate margin
margin=10*log10(N*Ex_bar/sum(en));