信道容量迭代算法

啦啦啦啦，开始菜鸟的第一篇博客。
这次的实验对我来说真是太艰难了，一部分因为对信道容量的计算方式本身就不熟，一部分原因是不太看得懂算法流程的意思，另一部分原因就是对MATLAB的矩阵运用不熟练。网上找的代码很多，但是没有什么解析，而且有的太繁琐，有的错误比较多。
信道容量和输入概率分布与信道转移概率有关，本次实验转移概率固定，是要求最佳输入概率分布。求信道容量就是求平均互信息量的最大值。在循环过程中，输入概率分布和信道容量不断迭代。其实至于数学表达式如何转化为算法步骤的，我看了半天还不是很清楚，可能复习完一遍线性代数概率论以后，回过头来看，我能明白点吧。
那就先尝试把复杂公式用程序实现，直接根据算法步骤转化MATLAB代码：
一开始看到这些奇怪的式子我真的是头大，$\sum$，i,j 这个求和那个求和，MATLAB这列乘以那列乱七八糟都理不出头绪。
下面是我的代码和一些详细的注释。

clc;clear;
%检验矩阵
%[0.98 0.02;0.05 0.95][0.6 0.4;0.05 0.95][0.8 0.15 0.05;0.05 0.15 0.8]
%[0.99 0.01 0;0.005 0.99 0.005;0 0.01 0.99]
Pij=input('请输入信道转移矩阵P')
[r,s]=size(Pij);%r为信源数 s为信道数
disp('原始信源分布:');
Pi=repmat(1/r,1,r);%产生信源Pi
disp(Pi)
offset=0.00000001;%误差
C(1)=-100000;%相当于无穷小
for k=2:10000
    %计算fai
    %无论j是几，转移矩阵的每一行都要×对应的输入概率，故将输入概率扩展为ij型矩阵后进行点乘
   Pi_copy=repmat(Pi',1,s);
   faizi=Pij.*Pi_copy;
   faimu=sum(faizi,1);%求和为i 则对结果的一列求和
   faimu_copy=repmat(faimu,r,1);%因为每一列除的值都是一样的
   fai=faizi./(faimu_copy);
   
   %计算新信源Pi
   Pizi=exp(sum(Pij.*log(fai+eps),2));%对j求和 行相加 转移概率为0 会出现NaN的数据,要进行特殊处理
   Pimu=sum(Pizi,1);%Pimu经过了两次求和 变成了标量数字
   Pi=Pizi/(Pimu);
   Pi=Pi'%方便下次迭代
   
   %计算信道容量
   C(k)=log(Pimu)/log(2);%注意信道容量以2为底
   if (C(k)-C(k-1))/C(k)<offset
       break;
   end;
end
disp('迭代次数：k='),disp(k-1)
disp('最大信道容量时的信源分布:p='),disp(Pi)
disp('最大信道容量:C='),disp(C(k))

这是运行结果：

遇到的问题：前三组数据通过检验，最后一组出错，并且出现NaN。第四组转移概率出现零。 首先想到的是怎么对零进行特殊处理，可是整体的算法思路看下来，并没有哪里可以更改，不管是不是零都是一样的。于是将0变成非常非常小的数。分母加上eps避免除以零，无效。改变循环次数，发现在第一次循环中迭代Pi时log出错，在工作区判断出是出现了log(0)，故加eps。
总结：
1.repmat()函数可以扩展矩阵，方便实现点乘，避免了矩阵乘法的繁琐。
2.sum()作用于方便实现i,j的求和。
3.ij分不清怎么乘可以尝试固定一个变量看。