function [ps C] = capacity_blahut(Q)
% Input: Q = channel transition probability matrix
% Output: C = channel capacity
% ps = row vector containing pmf that achieves capacity
tl = 1e-8; % tolerance (for the stopping condition)
n = 1000; % max number of iterations (in case the stopping condition
% is "never" reached")
nx = size(Q,1); pT = ones(1,nx)/nx; % First, guess uniform X.
for k = 1:n
qT = pT*Q;
% Eliminate the case with 0
% Column-division by qT
temp = Q.*(ones(nx,1)*(1./qT));
%Eliminate the case of 0/0
l2 = log2(temp);
l2(find(isnan(l2) | (l2==-inf) | (l2==inf)))=0;
logc = (sum(Q.*(l2),2))';
CT = 2.^(logc);
A = log2(sum(pT.*CT)); B = log2(max(CT));
if((B-A)