wongpowerallo.m

ofdmaz自适应资源分配

% Author: Ian C. Wong 
% Copyright (C) 2004   Ian C. Wong
% 
% 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.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
% 
% You may reach the author at [e-mail].
% Or visit his website at www.ece.utexas.edu/~iwong
function [p]= wongpowerallo(channel,suballo,N,K,Ptot,noise,gamma)
global oops
Nk = sum(suballo');
H = channel.*suballo/noise;
[Hsorted,I] = sort(H');
Hsorted = Hsorted';
n = [2:N];
k = [2:K];
V = zeros(1,K);
W = ones(1,K);
for k = 1 : K,
    id = find(Hsorted(k,:));
    Hk1 = Hsorted(k,id(1));
    if k == 1,
        H11 = Hk1;
    end;
    if length(id) > 1,   
        for n = id(2) : N,
            V(k) = V(k) + (Hsorted(k,n) - Hk1)/(Hk1*Hsorted(k,n));
            W(k) = W(k)*(Hsorted(k,n)/Hk1)^(1/(Nk(k)));
            a(k) = -Nk(1)/Nk(k)*Hk1*W(k)/(H11*W(1));
            b(k) = Nk(1)/(H11*W(1))*(W(k) - W(1) + H11*V(1)*W(1)/Nk(1) - Hk1*V(k)*W(k)/Nk(k));
        end;
    else,
        V(k) = 0;
        W(k) = 1;
        a(k) = -Nk(1)/Nk(k)*Hk1*W(k)/(H11*W(1));
        b(k) = Nk(1)/(H11*W(1))*(W(k) - W(1) + H11*V(1)*W(1)/Nk(1) - Hk1*V(k)*W(k)/Nk(k));
    end;
end;
a(1) = 0;
b(1) = Ptot;
A = diag(a) + not([0; ones(K-1,1)]*[0 ones(1,K-1)]);
% p = (A\b')'
p = zeros(length(a),1);
bcut = b(2:length(b));
acut = a(2:length(a));
p(1) = (Ptot-sum(bcut./acut))/(1-sum(1./acut));
p(2:length(a)) = ((bcut)-p(1))./acut;

if any(p < zeros(size(p)) | ~isreal(p)),
    oops = oops + 1;
    p = 1/K*ones(1,K);
end;