channelcapacity.m

信道容量的计算

clc;clear all;
r = input('the number of input r='); 
s = input('the number of output s='); 
p_matrix=zeros(r,s);  
fprintf('input the value of p_matrix:\n')
for i=1:r
    for j=1:s
    p_matrix(i,j)=input('p_matrix=');
    if p_matrix(i,j)<0
    error('Do not accord with the distribution')
    end
    end
end
for i=1:r   
    sum_line(i)=0;
    for j=1:s
        sum_line(i)=sum_line(i)+p_matrix(i,j);
    end
end
for i=1:r    
if (sum_line(i)<=0.999999||sum_line(i)>=1.000001)
   error('Do not accord with the distribution')
end
end
e=input('input the iterative precision:');
for i=1:r
    p(i)=1.0/r;  % the initial probability accord with uniform distribution
end
for j=1:s   %the sum of q(j)
    q(j)=0;
   for i=1:r
       q(j)=q(j)+p(i)*p_matrix(i,j);
   end
end        
for i=1:r    %the sum of a(i)
    d(i)=0;
    for j=1:s
        if(p_matrix(i,j)==0)
           d(i)=d(i)+0;
        else
         d(i)=d(i)+p_matrix(i,j)*log(p_matrix(i,j)/q(j));
       end
   end
   a(i)=exp(d(i));
end
u=0;
for i=1:r  %calculate the value of u
    u=u+p(i)*a(i);
end
IL=log2(u); %calculate the value of IL
IU=log2(max(a));%calculate the value of IU
n=1;
while((IU-IL)>=e) %the iterative calculation
    for i=1:r
        p(i)=p(i)*a(i)/u; %refresh the value of p(i)
    end
    for j=1:s   %calculate the value of q(j)
        q(j)=0;
        for i=1:r
            q(j)=q(j)+p(i)*p_matrix(i,j);
        end
    end
    for i=1:r   %calculate the value of a(i)
        d(i)=0;
        for j=1:s
            if(p_matrix(i,j)==0)
              d(i)=d(i)+0;
            else
             d(i)=d(i)+p_matrix(i,j)*log(p_matrix(i,j)/q(j));
            end
       end
       a(i)=exp(d(i));
    end
    u=0;
    for i=1:r  %calculate the value of u
        u=u+p(i)*a(i);
    end
      IL=log2(u); %calculate the value of IL
      IU=log2(max(a));%calculate the value of IU
      n=n+1;
end
fprintf('the channel matrix is:\n');
disp(p_matrix);
fprintf('the iterative number n=%d\n',n);
fprintf('the channel capacity C=%fbit/symbol',IL);