RSS Feed

Blahut-Arimoto algorithm in Matlab

0 Comments
Posted by wkshum on 13/08/2012 at 11:53 am

For a discrete memoryless channel p(y|x), the capacity is defined as

    \[ C = \max_{r(x)} I(X;Y) = \max_{r(x)} \sum_{x} \sum_{y} r(x)p(y|x) \log \frac{r(x)p(y|x)}{r(x) \sum_{x'} r(x') p(y|x')} \]

where X and Y denote the input and output variables of the channel respectively, and the maximization is taken over all input distributions r(x).

Given a channel transition matrix whose (i,j)-entry is the conditional probability p(j|i), the Blahut-Arimoto algorithm computes the capacity of the discrete memoryless channel, and the input distribution r(x) that attains the maximum.

 

Reference: Chapter 9 in the book Information Theory and Network Coding by Raymond Yeung.

function [C r] = BlahutArimoto(p)
% Capacity of discrete memoryless channel
% by Blahut-Arimoto algorithm
% Input
% p: m x n matrix
% p is the transition matrix for a channel with m inputs and n outputs
% 
% The input matrix p should contain no zero row and no zero column.
%
% p(i,j) is the condition probability that the channel output
% is j given that the input is i
% (i=1,2,...,m and j = 1,2,...,n)
%
%
% Output 
% capacity : capacity in bits
% r: channel input distribution which achieves capacity
%
% For example, the transition matrix for the erasure channel is
% can be calculated as
% e = 0.5;
% p = [1-e e 0; 0 e 1-e]; % conditional prob. for erasure channel
% The capacity can be calculated by BlahutArimoto(p), and is equal to 1-e
%
% Check that the entries of input matrix p are non-negative
if ~isempty(find(p < 0))
 disp('Error: some entry in the input matrix is negative')
 C = 0; return;
end
% Check that the input matrix p does not have zero column
column_sum = sum(p);
if ~isempty(find(column_sum == 0))
 disp('Error: there is a zero column in the input matrix');
 C= 0; return;
end
% Check that the input matrix p does not have zero row
row_sum = sum(p,2);
if ~isempty(find(row_sum == 0))
 disp('Error: there is a zero row in the input matrix');
 C= 0; return;
else
 p = diag(sum(p,2))^(-1) * p; % Make sure that the row sums are 1
end
[m n] = size(p);
r = ones(1,m)/m; % initial distribution for channel input
q = zeros(m,n);
error_tolerance = 1e-5/m;
r1 = [];
for i = 1:m
 p(i,:) = p(i,:)/sum(p(i,:));
end
for iter = 1:10000
 for j = 1:n
  q(:,j) = r'.*p(:,j);
  q(:,j) = q(:,j)/sum(q(:,j));
 end

 for i = 1:m
  r1(i) = prod(q(i,:).^p(i,:));
 end
 r1 = r1/sum(r1);

 if norm(r1 - r) < error_tolerance
  break
 else
  r = r1;
 end
end
C = 0;
for i = 1:m
 for j = 1:n
  if r(i) > 0 && q(i,j) > 0
   C = C+ r(i)*p(i,j)* log(q(i,j)/r(i));
 end
 end
end
C = C/log(2); % Capacity in bits

For example, we can compute the capacity of an erasure channel with erasure probability 0.1 as follows.

 

e = .1; p  = [1-e e 0; 0 e 1-e]
p =
 0.9000 0.1000 0
 0 0.1000 0.9000
>> [C r] = BlahutArimoto(p)
C =
0.9000
r =
0.5000 0.5000

 

The capacity is C=0.9. It is achieved by the channel input probability distribution r(1)=r(2)=0.5.

You can leave a comment, or trackback from your own site.

0 Comments

You can be the first to comment!

Leave a Reply

Your email address will not be published. Required fields are marked *

*

* Copy this password:

* Type or paste password here:

55,796 Spam Comments Blocked so far by Spam Free Wordpress

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>