Water-filling example in Matlab/Octave

waterfilling.m

% Number of channels
N = 10;
N0 = 1; % Normalized noise level
SNR_dB = 10; % Signal-to-noise ratio in dB
P = 10**(SNR_dB / 10);


% The channel specific gains drawn from Gaussian distribution
g = abs(randn(N, 1));

% Bisection search for alpha
alpha_low = min(N0 / g);
alpha_high = (P + sum(N0 / g)) / N; % Initial high (upper bound)

stop_threshold = 1e-5; % Stop threshold

% Iterate while low/high bounds are futher than stop_threshold
while(abs(alpha_low - alpha_high) > stop_threshold)
    alpha = (alpha_low + alpha_high) / 2; % Test value in the middle of low/high

    p = 1 / alpha - N0 / g;

    p(p < 0) = 0; % Consider only positive powers

    if (sum(p) > P)
        alpha_low = alpha;
    else
        alpha_high = alpha;
    end
end

% Print the achievable rate in bits/s
rates = log2(1 + g*p/N0);
disp(['Achievable rate ', num2str(sum(rates(:))), ' bits/s'])