# Water-filling algorithm implemented using Matlab & Octave

In this post, I'm going to show you how to implement the water-filling algorithm using Matlab & Octave.

Matlab is a great numerical computations tool, which has been widely adopted by engineers and scientists for all sorts of numerical computations.

Water-filling is a generic method for allocating transmission power over (equalizing) communications channels in the absence of inter-channel interference (ICI). For an in-depth overview of the water-filling algorithm, check our post showing the algorithm behavior and derivation.

## Basic concept

We can write the throughput for $$N$$ orthogonal channels as the sum-rate over all the channels

$\sum_{n=1}^N \text{log}_2 (1 + \text{SNR}_n)$

We use the sum power budget ($$P$$), which limits the available power across the channels as

$\sum_{n=1}^N p_n = P$

Further, we must define the $$\text{SNR}_n$$ in terms of the $$p_n$$ as

$\text{SNR}(p_n) = {g_n p_n \over N_0}$

where $$g_n$$ denotes the channel gain for the $$n^\text{th}$$ channel, $$p_n$$ is the corresponding power allocation and $$N_0$$ is a fixed noise-level determined by the given SNR.

Now, we can write the capacity as

$C = \sum_{n=1}^N \text{log}_2 (1 + \text{SNR}(p^*_n))$

where $$p^*_n$$ is the sum-rate maximizing power allocation for channel $$n$$.  Now, we have the basic premise set up to come up with an algorithm that finds the optimal power allocation, namely, the water-filling algorithm.

## The algorithm

We'll skip the algorithm derivation and skip directly, our power allocation, for each channel $$n = 1,\ldots,N$$, is simply given by $p_n = \Big({1 \over {\alpha} } - {N_0 \over g_n}\Big)^{+}$

We can just search for such $$\alpha$$ that satisfies the power budget allocation $\sum_{n=1}^N p_n = P$

If you are interested in the complete derivation of the algorithm, check our previous post.

## Water-filling in Matlab / Octave

Now, we are ready to jump into the implementation of the water-filling algorithm using Matlab & Octave. We avoid using any Matlab-specific functionality everything here should apply to both Matlab & Octave.

We begin, by defining the system parameters


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


These are the basic parameters that define the SNR ranges and the system's complexity.

Next, we need to generate channel gains for all $$N$$ channels. We simply draw the from Gaussian distribution. There are more complex and realistic channel models that could be used, but those are not within the scope of this article.

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

We start implementing the actual algorithm, by initializing the bounds for our search space. Stop threshold is used to set the point when the algorithm has converged. Our test condition is $$| \alpha_\text{low} - \alpha_\text{high} | < \epsilon_\text{stop}$$

% 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

The complete search (bisection) loop is given by

% 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

We can see that adjusting $$\alpha$$ acts as a scaling coefficient for the power allocation. Thus, increasing $$\alpha$$ always lowers power allocation and vice versa. Keeping this in mind, we will implement a bisection search for the power allocation also known as the water-filling, where we search for the water level (power scaling), by bisecting the search range for $$\alpha$$.

We were able to write a custom algorithm for solving the channel capacity-achieving power allocation with a couple of lines of code, just by having a look at the KKT conditions and identifying the solution. This gives us an excellent place to start for more complex systems.

The complete source can be found on Github.