chcap1.m

％ 信道容量C的迭代算法 ％ ％ 函数说明： ％ ％ [CC,Paa]=ChannelCap(P,k) 为信道容量函数 ％ ％ 变量说明： ％ ％ P：输入的正向转移概率矩阵

% 信道容量C的迭代算法                              %
% 函数说明：                                       %
% [CC,Paa]=ChannelCap(P,k) 为信道容量函数        	%
% 变量说明：                                       %
% P：输入的正向转移概率矩阵，k：迭代计算精度       %
% CC：最佳信道容量，Paa：最佳输入概率矩阵          %
% Pa：初始输入概率矩阵，Pba：正向转移概率矩阵    	%
% Pb：输出概率矩阵                                	%
% C：初始信道容量， r：输入符号数，s：输出符号数   %
%************************************************  %
function[CC,Paa]=ChannelCap(P,K)
%初始化Pa
P=input('请输入转移概率矩阵P:[直接回车表示P=[1/3,1/3,1/6,1/6;1/6,1/3,1/6,1/3]\n');
if isempty(P)
P=[1/3,1/3,1/6,1/6;1/6,1/3,1/6,1/3];
end
[r,s]=size(P);
Pa=(1/(r+eps))*ones(1,r);
Pba=P;
%进行迭代计算
n=1;
C=0;
CC=1;
k=input('请输入迭代计算精度k:[直接回车表示k=0.0001]\n');
if isempty(k)
k=0.0001;
end
while abs(CC-C)>=k
    n=n+1;
%求Pb
Pb=zeros(1,s);
for j=1:s
  for i=1:r
        Pb(j)=Pb(j)+Pa(i)*Pba(i,j);
    end
end
 %求Paa
a=zeros(1,r);
b=zeros(1,r);
  for i=1:r
     for j=1:s
        a(i)=a(i)+Pba(i,j)*log2((Pba(i,j)+eps)/(Pb(j)+eps));
    end
b(i)=exp(a(i));
end
sum1=0;
for i=1:r
sum1=sum1+Pa(i)*(b(i)+eps);
end
sum2=0;
for i=1:r
    if sum2<b(i)
        sum2=b(i);
    end
end
for i=1:r
    Paa(i)=Pa(i)*b(i)/(sum1+eps);
end
%求信道容量
C=log(sum1);
CC=log(sum2);
Pa=Paa;
end
%打印输出结果
s1='最佳输入概率分步Paa';
s2='信道容量C';
s3='迭代次数n';
s4='输入符号r';
s5='输出符号s';
s6='迭代精度k';
for i=1:r
    B{i}=i;
end
disp(s1),disp(B),disp(Paa);
disp(s4),disp(r);
disp(s5),disp(s);
disp(s2),disp(CC);
disp(s6),disp(k);
disp(s3),disp(n);