channelcap.m

用matlab实现的一种求信道容量的改进方法

% 信道容量C的迭代算法                            %
% 函数说明：                                     %
% [CC,Paa]=ChannelCap(P,k) 为信道容量函数        %
% 变量说明：                                     %
% P：输入的正向转移概率矩阵，k：迭代计算精度     % 
% CC：最佳信道容量，Paa：最佳输入概率矩阵        % 
% Pa：初始输入概率矩阵，Pba：正向转移概率矩阵    %          
% Pb：输出概率矩阵，Pab：反向转移概率矩阵        %
% C：初始信道容量， r：输入符号数，s：输出符号数 %
%************************************************%

function [CC,Paa]=ChannelCap(P,k)

% 提示错误信息 
if (length(find(P<0))~=0)  
    error('Not a prob.vector,negative component');  % 判断是否符合概率分布条件
end  

if (abs(sum(P')-1)>10e-10)  
    error('Not a prob.vector,components do not add up to 1')  % 判断是否符合概率和为1
end  

% 1)初始化Pa
[r,s]=size(P);
Pa=(1/(r+eps))*ones(1,r);
sumrow=zeros(1,r);
             
Pba=P;
     
% 2)进行迭代计算
n=0;
C=0; 
CC=1;

while abs(CC-C)>=k  % 迭代开始
      n=n+1;
      
% (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

% (2)再求Pab
suma=zeros(1,s);

for j=1:s
    for i=1:r
          Pab(j,i)=Pa(i)*Pba(i,j)/(Pb(j)+eps);
          suma(j)=suma(j)+Pa(i)*Pba(i,j)*log2((Pab(j,i)+eps)/(Pa(i)+eps));
    end
end

% 3)求信道容量C
C=sum(suma);

% 4)求下一次Pa,即Paa
L=zeros(1,r);

for i=1:r
    for j=1:s
          L(i)=L(i)+Pba(i,j)*log(Pab(j,i)+eps);   
    end
    a(i)=exp( L(i));  
end
    
Paa=a./(sum(a)+eps); 

% 5)求下一次C,即CC
CC=log2(sum(a)+eps);
Pa=Paa;

end  % 迭代结束

% 打印输出结果
 s0='很好！输入正确，迭代结果如下：';
 s1='最佳输入概率分布Pa:';
 s2='信道容量C:';
 s3='迭代次数n:';
 s4='输入符号数r:';
 s5='输出符号数s:';
 s6='迭代计算精度k:';
 
 for i=1:r
       B{i}=i;    
 end
 
disp(s0);
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);