信息论matlab的实验用到的编码

Matlab源程序附录A信息熵%函数说明：%%H=entropy(P,r)为信息熵函数%%P为信源的概率矢量,r为进制数%%H为信息熵%%******************************%functionH=entropy(P,r)if(length(find(P=0))~=0)error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件endif(abs(sum(P)-1)10e-10)error('Notaprob.vector,componentdonotaddupto1');endH=(sum(-P.*log2(P)))/(log2(r)+eps);附录B离散无记忆信道容量的迭代计算%信道容量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(P0))~=0)error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件endif(abs(sum(P')-1)10e-10)error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1信息论基础及应用294end%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;whileabs(CC-C)=kn=n+1;%(1)先求PbPb=zeros(1,s);forj=1:sfori=1:rPb(j)=Pb(j)+Pa(i)*Pba(i,j);endend%(2)再求Pabsuma=zeros(1,s);forj=1:sfori=1:rPab(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%3)求信道容量CC=sum(suma);%4)求下一次Pa,即PaaL=zeros(1,r);sumaa=0;fori=1:rforj=1:sL(i)=L(i)+Pba(i,j)*log(Pab(j,i)+eps);enda(i)=exp(L(i));endsumaa=sum(a);fori=1:rPaa(i)=a(i)/(sumaa+eps);end%5)求下一次C,即CCCC=log2(sumaa);Pa=Paa;附录Matlab源程序295end%打印输出结果s0='很好！输入正确，迭代结果如下：';s1='最佳输入概率分布Pa:';s2='信道容量C:';s3='迭代次数n:';s4='输入符号数r:';s5='输出符号数s:';s6='迭代计算精度k:';fori=1:rB{i}=i;enddisp(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);附录CShannon编码%函数说明：%%[p,x]=array(P,X)为按降序排序的函数%%P为信源的概率矢量，X为概率元素的下标矢量%%p为排序后返回的信源的概率矢量%%x为排序后返回的概率元素的下标矢量%%*******************************************%function[p,x]=array(P,X)P=[P;X];[l,n]=size(P);fori=1:nmax=P(1,i);maxN=i;MAX=P(:,i);forj=i:nif(maxP(1,j))MAX=P(:,j);max=P(1,j);maxN=j;endendif(maxN1)信息论基础及应用296if(in)fork=(maxN-1):-1:iP(:,k+1)=P(:,k);endendendP(:,i)=MAX;endp=P(1,:);x=P(2,:);%shannon编码生成器%%函数说明：%%[W,L,q]=shannon(p)为shannon编码函数%%p为信源的概率矢量，W为编码返回的码字%%L为编码返回的平均码字长度，q为编码效率%%*****************************************%function[W,L,q]=shannon(p)%提示错误信息if(length(find(p=0))~=0)error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件endif(abs(sum(p)-1)10e-10)error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1end%1)排序n=length(p);x=1:n;[p,x]=array(p,x);%2)计算代码组长度ll=ceil(-log2(p));%3)计算累加概率PP(1)=0;n=length(p);fori=2:nP(i)=P(i-1)+p(i-1);end%4)求得二进制代码组W%a)将十进制数转为二进制数fori=1:nforj=1:l(i)temp(i,j)=floor(P(i)*2);P(i)=P(i)*2-temp(i,j);附录Matlab源程序297endend%b)给W赋ASCII码值，用于显示二进制代码组Wfori=1:nforj=1:l(i)if(temp(i,j)==0)W(i,j)=48;elseW(i,j)=49;endendendL=sum(p.*l);%计算平均码字长度H=entropy(p,2);%计算信源熵q=H/L;%计算编码效率%打印输出结果fori=1:nB{i}=i;end[n,m]=size(W);TEMP=32*ones(n,6);W=[W,TEMP];W=W';[n,m]=size(W);W=reshape(W,1,n*m);W=sprintf('%s',W);s0='很好！输入正确，编码结果如下：';s1='Shannon编码所得码字W:';s2='Shannon编码平均码字长度L:';s3='Shannon编码的编码效率q:';disp(s0);disp(s1),disp(B),disp(W);disp(s2),disp(L);disp(s3),disp(q);附录DFano编码%函数说明：%%[next_P,next_index,code_num]=compare(current_P,current_index)%%为比较函数，主要用于信源符号的分组%%current_P为当前分组的信源的概率矢量%%current_index为当前分组的信源的下标%%next_P为返回的下次分组的信源的概率矢量%信息论基础及应用298%next_index为返回的下次分组的信源的下标%%code_num为返回的ASCII值%%*********************************************************************%function[next_P,code_num,next_index]=compare(current_P,current_index);n=length(current_P);add(1)=current_P(1);%1)求概率的依次累加和fori=2:nadd(i)=0;add(i)=add(i-1)+current_P(i);end%2)求概率和最接近的两小组s=add(n);fori=1:ntemp(i)=abs(s-2*add(i));end[c,k]=min(temp);%3)对分组的信源赋ASCII值if(current_index=k)next_index=current_index;code_num=48;next_P=current_P(1:k);elsenext_index=current_index-k;code_num=49;next_P=current_P((k+1):n);end%fano编码生成器%%函数说明：%%[W,L,q]=fano(P)为fano编码函数%%P为信源的概率矢量，W为编码返回的码字%%L为编码返回的平均码字长度，q为编码效率%%*****************************************%function[W,L,q]=fano(P)%提示错误信息if(length(find(P=0))~=0)error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件endif(abs(sum(P)-1)10e-10)error('Notaprob.vector,componentdonotaddupto1')%判断是否符合概率和为1end%1)排序n=length(P);附录Matlab源程序299x=1:n;[P,x]=array(P,x);%2)将信源符号分组并得到对应的码字fori=1:ncurrent_index=i;j=1;current_P=P;while1[next_P,code_num,next_index]=compare(current_P,current_index);current_index=next_index;current_P=next_P;W(i,j)=code_num;j=j+1;if(length(current_P)==1)break;endendl(i)=length(find(abs(W(i,:))~=0));%得到各码字的长度endL=sum(P.*l);%计算平均码字长度H=entropy(P,2);%计算信源熵q=H/L;%计算编码效率%打印输出结果fori=1:nB{i}=i;end[n,m]=size(W);TEMP=32*ones(n,5);W=[W,TEMP];W=W';[n,m]=size(W);W=reshape(W,1,n*m);W=sprintf('%s',W);s0='很好！输入正确，编码结果如下：';s1='Fano编码所得码字W:';s2='Fano编码平均码字长度L:';s3='Fano编码的编码效率q:';disp(s0);disp(s1),disp(B),disp(W);disp(s2),disp(L);disp(s3),disp(q);信息论基础及应用300附录EHuffman编码Huffman编码（1）%huffman编码生成器%%函数说明：%%[W,L,q]=huffman(P)为huffman编码函数%%P为信源的概率矢量，W为编码返回的码字%%L为编码返回的平均码字长度，q为编码效率%%*****************************************%function[W,L,q]=huffman(P)if(length(find(P=0))~=0)error('Notaprob.vector,negativecomponent');%判断是否符合概率分布条件endif(abs(sum(P)-1)10e-10)er