clc;clear all;
N = input('輸入信源符號(hào)X的個(gè)數(shù)N='); 
M = input('輸出信源符號(hào)Y的個(gè)數(shù)M='); 
p_yx=zeros(N,M);   %程序設(shè)計(jì)需要信道矩陣初始化為零
fprintf('輸入信道矩陣概率\n')
for i=1:N
     for j=1:M
     p_yx(i,j)=input('p_yx=');%輸入信道矩陣概率
     if p_yx(i)<0
     error('不符合概率分布')
     end
     end
end
for i=1:N   %各行概率累加求和
     s(i)=0;
     for j=1:M
         s(i)=s(i)+p_yx(i,j);
     end
end
for i=1:N     %判斷是否符合概率分布
if (s(i)<=0.999999||s(i)>=1.000001)
   error('不符合概率分布')
end
end
b=input('輸入迭代精度：');%輸入迭代精度
for i=1:N
     p(i)=1.0/N;   %取初始概率為均勻分布
end
for j=1:M   %計(jì)算q(j)
     q(j)=0;
   for i=1:N
       q(j)=q(j)+p(i)*p_yx(i,j);
   end
end        
for i=1:N     %計(jì)算a(i)
     d(i)=0;
     for j=1:M
         if(p_yx(i,j)==0)
           d(i)=d(i)+0;
         else
         d(i)=d(i)+p_yx(i,j)*log(p_yx(i,j)/q(j));
       end
   end
   a(i)=exp(d(i));
end
u=0;
for i=1:N   %計(jì)算u
     u=u+p(i)*a(i);
end
IL=log2(u); %計(jì)算IL
IU=log2(max(a));%計(jì)算IU
n=1;
while((IU-IL)>=b) %迭代計(jì)算
     for i=1:N
         p(i)=p(i)*a(i)/u; %重新賦值p(i)
     end
     for j=1:M   %計(jì)算q(j)
         q(j)=0;
         for i=1:N
             q(j)=q(j)+p(i)*p_yx(i,j);
         end
     end
     for i=1:N   %計(jì)算a(i)
         d(i)=0;
         for j=1:M
             if(p_yx(i,j)==0)
               d(i)=d(i)+0;
             else
             d(i)=d(i)+p_yx(i,j)*log(p_yx(i,j)/q(j));
             end
       end
       a(i)=exp(d(i));
     end
     u=0;
     for i=1:N   %計(jì)算u
         u=u+p(i)*a(i);
     end
       IL=log2(u); %計(jì)算IL
       IU=log2(max(a));%計(jì)算IU
       n=n+1;
end
fprintf('信道矩陣為：\n');
disp(p_yx);
fprintf('迭代次數(shù)n=%d\n',n);
fprintf('信道容量C=%f比特/符號(hào)',IL);