www.pudn.com > MultiuserOFDM_0_1.zip > waterfillingPowerAllo.m

% Copyright (c) 2003-2004 The University of Texas 
% All Rights Reserved. 
%   
% This program is free software; you can redistribute it and/or modify 
% it under the terms of the GNU General Public License as published by 
% the Free Software Foundation; either version 2 of the License, or 
% (at your option) any later version. 
%   
% This program is distributed in the hope that it will be useful, 
% but WITHOUT ANY WARRANTY; without even the implied warranty of 
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
% GNU General Public License for more details. 
%   
% The GNU Public License is available in the file LICENSE, or you 
% can write to the Free Software Foundation, Inc., 59 Temple Place - 
% Suite 330, Boston, MA 02111-1307, USA, or you can find it on the 
% World Wide Web at http://www.fsf.org. 
%   
% Programmers:	Zukang Shen 
%  
% The authors are with the Department of Electrical and Computer 
% Engineering, The University of Texas at Austin, Austin, TX. 
% They can be reached at shen@ece.utexas.edu. 
 
function PowerAllo= waterfillingPowerAllo(Ptot,eigenVec) 
eigenVec=eigenVec+1e-10; 
eigenVec=eigenVec(:); 
[sortedEigVec,oriPosition]=sort(eigenVec); 
 
EigVecTemp=sortedEigVec; 
 
begin=1; 
while sum((EigVecTemp-EigVecTemp(1)*ones(length(EigVecTemp),1))./(EigVecTemp*EigVecTemp(1)))>Ptot 
    EigVecTemp=EigVecTemp(2:length(EigVecTemp)); 
    begin=begin+1; 
end 
 
P1=(Ptot-sum((EigVecTemp-EigVecTemp(1)*ones(length(EigVecTemp),1))./(EigVecTemp*EigVecTemp(1))))/length(EigVecTemp); 
for i=1:begin-1 
    PowerAlloOrdered(i)=0; 
end 
 
PowerAlloOrdered(begin)=P1; 
 
for i=2:length(EigVecTemp) 
    PowerAlloOrdered(begin+i-1)=P1+(EigVecTemp(i)-EigVecTemp(1))/(EigVecTemp(i)*EigVecTemp(1)); 
end 
 
for i=1:length(eigenVec) 
    PowerAllo(oriPosition(i))=PowerAlloOrdered(i); 
end 
 
sumCapa=0; 
for i=1:length(eigenVec) 
    sumCapa=sumCapa+log2(1+PowerAllo(i)*eigenVec(i)); 
end 
 
mu=P1+1/EigVecTemp(1); 
capacity=sum(log2(mu*EigVecTemp));