function?[gnen_barbn_barNstarb_bar]=waterfill（PSNRmfbEx_barNtotgap）

%?EE379C?2001-2002?Spring
%?EE379C?1999-2000?Spring?
%?EE379C?1997-1998?Spring?
%?Written?by?Atul?Salvekar?Edited?by?Wonjong?Rhee
%
%?function?[gnen_barbn_barNstarb_bar]=waterfill（PSNRmfbEx_barNgap）
%
%?P?is?the?pulse?response
%?SNRmfb?is?the?SNRmfb?in?dB
%?Ex_bar?is?the?normalized?energy
%?Ntot?is?the?total?number?of?real/complex?subchannels?Ntot>2
%?gap?is?the?gap?in?dB
%
%?gn?is?channel?gain
%?en_bar?is?the?energy/dim?in?the?nth?subchannel
%?bn_bar?is?the?bit/dim?in?the?nth?subchannel
%?Nstar?is?the?number?of?subchannel?used
%?b_bar?is?the?bit?rate

%?dB?into?normal?scale
Noise_var=Ex_bar*（norm（P）^2）/（10^（SNRmfb/10））;
gap=10^（gap/10）;

%?initialization
en=zeros（1Ntot）;
bn=zeros（1Ntot）;
gn=zeros（1Ntot）;
Hn?=?zeros（1Ntot）;

%?subchannel?center?frequencies
f=-1/2+1/Ntot: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?+?....??or?P（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;
%plot（gn）

%%%%%%%%%%%%%%%%%%%%%%%
%?Now?do?waterfilling?%
%%%%%%%%%%%%%%%%%%%%%%%

%sort
[gn_sorted?Index]=sort（gn）;??%?sort?gain?and?get?Index

gn_sorted?=?fliplr（gn_sorted）;%?flip