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

Matlab源程序附录A信息炳%函数说明：％%H=entropy(P/)为信息炳函数％%P为信源的概率矢量,r为进制数%%11为信息炳%%******************************%fiinctionH=entiopy(Pj)if(length(fuid(P<=0)广=0)error(TSIotapiob.vector.negativecomponent1);%判断是否符合概率分布条件endif(abs(sum(P)-1)>10e-10)error(TSIotaptob.vector,componentdonotaddupto1');endH=(sum(-P.*log2(P)))/(log2(r)+eps);附录B离散无记忆信道容量的迭代计算TOC\o"1-5"\h\z%信道容量c的迭代算法％%函数说明：%%[CC,Paa]=ChaimelCap(P,k)为信道容量函数%%变量说明：％%P：输入的正向转移概率矩阵，k：迭代计算精度％%CC：最佳信道容量，Paa：最佳输入概率矩阵％%Pa：初始输入概率矩阵，Pba：正向转移概率矩阵%%Pb：输出概率矩阵，Pab：反向转移概率矩阵%%C：初始信道容量，1-：输入符号数，s：输出符号数%%**************************************************%fiinction[CC,Paa]=ChannelCap(P,k)%提示错误信息if(lengtli(fmd(P<0))~=0)error(TSIotapiob.vector.negativecomponent*);%判断是否符合概率分布条件endif(abs(sum(P,)-1)>10e~10)error(TSIotaprob.vector,componentdonotadduptoI1)%判断是否符合概率和为1end%1)初始化Pa[r,s]=size(P);Pa=(l/(r+eps))*ones(l,r);sumrow=zeros(1,r);Pba=P;%2)进行迭代计算n=0;C=0；CC=1;whileabs(CC-C)>=kn=n+l;%⑴先求PbPb=zeros(l,s);forj=l:sfori=l:rPb(j)=Pb(j)+Pa(i)*Pba(i,j);endend%(2)再求Pabsuma=zeros(l,s);forj=l:sfori=l:rPab(j,i)=Pa(i)*Pba(ij)/(Pb(j)+eps);suma(j)=suma(j)+Pa①*Pba(iJ)*log2((Pab(j,i)+eps)/(Pa(i)+eps));end%3)求信道容量CC=sum(sunia);%4)求下—*次Pa,即PaaL=zeios(l,r);sumaa=0;fbr1=1:rfbrj=l:sL(i)=L(i)+Pba(ij)*log(Pab(j,i)+eps);endam=exp(L(i));endsumaa=sum(a);fbr1=1:rPaa(i)=a(i)/(sumaa+eps);end%5)求下一次C,即CCCC=log2(sumaa);Pa=Paa;end%打印输出结果s0=，很好！输入正确，迭代结果如下：，;sl=，最佳输入概率分布Pa:';s2=，信道容量C:';S3。迭代次数H?；s4=，输入符号数1：'；s5=输出符号数s「；s6=，迭代计算精度k「；fbr1=1:rB{1}=1；enddisp(sO);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]=aiiay(P,X)为按降序排序的函数%P为信源的概率矢量，X为概率元素的下标矢量%%p为排序后返回的信源的概率矢量%%乂为排序后返回的概率元素的下标矢量%%*******************************************%fiinction[p.x]=array(P,X)P=[P；X]；[Ln]=size(P);fbri=l:nmax=P(lj);niaxN=i;MAX=P(:,i);fbrj=i:nif(niax<P(l,j))MAX=P(:J);max=P(l,j);maxN^j;endendif(maxN>l)

if(Kn)fbrk=(maxNT):T:iP(:,k+l)=p(:,k);endendendP(:,i)=MAX;endP=P(1,:);x=P(2,:);%%%%%%%%%%%[W,L,q]=shannon(p)为shaimon编码函数%p为信源的概率矢量，W为编码返回的码字%1为编码返回的平均码字长度，q为编码效率%*****************************************%function[W.L,q]=shaimon(p)%提示错误信息if(lengtli(fmd(p<=0))~=0)error(TSIotaprob.vector.negativecomponent1);%判断是否符合概率分布条件endif(abs(sum(p)-1)>1Oe-10)error(TSIotapiob.vector,componentdonotadduptoI1)%判断是否符合概率和为1end%1)排序n=length(p);x=l:n;[p,x]=array(p,x);%2)计算代码组长度1l=ceil(~log2(p));%3)计算累加概率PP⑴=0;n=length(p);fbri=2:nP(i>p(i-l)+p(i-l);end%4)求得二进制代码组W%a)将十进制数转为二进制数fori=l:ntemp(i,j)=floor(P(i)*2);P(i)=P(i)*2-temp(ij);endend%妙给可赋ASCII码值，用于显示二进制代码组Wfbri=l:nfbrj=l:l(i)if(temp(ij)==O)W(i,j)=48;elseW(i,j)=49;endendendL=sum(p.*l);%计算平均码字长度H=entropy(p,2);%计算信源炳q=H'L;%计算编码效率%打印输出结果fbri=l:nend[n,m]=size(W);TEMP=32*ones(n,6);W=[W.TEMP];W=W';[n,m]=size(W);W=reshape(W,1W=spnntR'%s',W);sO=，很好！输入正确，编码结果如下：,；sl='Shaimon编码所得码字W:';s2=Shannon编码平均码字长度L「；s3=,Shaiuion编码的编码效率q:'；disp(sO);disp(sl),disp(B),disp(W);disp(s2).disp(L);disp(s3)Jisp(q);附录DFano编码%函数说明：%%[nexcP.next_mdex,code_num]=compare(cuiTent_P.current_index)%%为比较函数，主要用于信源符号的分组%%current.P为当前分组的信源的概率矢量%%current_mdex为当前分组的信源的下标%%nexcP为返回的下次分组的信源的概率矢量%

%next.index为返回的下次分组的信源的下标%%code_num为返回的ASCII值%%*********************************************************************%function[next_P,code-nuni4iext_mdex]=compare(cuiTenCPxurrent_index);n=length(current_P);add(l)=current_P(1);%1)求概率的依次累加和fbri=2:nadd⑴=0;add(i)=add(iT)+cuirent_P(i);end%2)求概率和最接近的两小组s=add(n);fbri=l:ntemp(i)=abs(s-2*add(i));end[c,k]=niui(temp);%3)对分组的信源赋ASCII值if(current_indexv=k)nextmdex=currentindex;code_num=48;next_P=current_P(l:k);elsenextmdex=currentindex-k:code_num=49;next_P=current_P((k+l):n);end%%%%%%%%%%[W,L,q]=fhno(P)为fano编码函数%P为信源的概率矢量，w为编码返回的码字%1为编码返回的平均码字长度，q为编码效率%%*****************************************%function[W,L,q]=fhno(P)%提示错误信息if(length(fuid(P<=0))~=0)error(TSIotapiob.vector.negativecomponent*);%判断是否符合概率分布条件endif(abs(sum(P)-1)>10e-10)error(TSIotapiob.vector,componentdonotadduptoT)%判断是否符合概率和为1end%1)排序n=length(P);x=l:n;[P,x]=array(P5x);%2)将信源符号分组并得到对应的码字fori=l:ncunent_mdex=i;j=l;cunent_P=P;while1[next_P,code_num,next_index]=compare(cuiTent_P,current_nidex);currentindex=nextmdex;currentP=nextP;W(iJ)=code_num;J=J+1;if(length(cuiTent_P)==1)break;endendl(i)=lengtli(fiiid(abs(W(i,:))~=0));%得到各码字的长度endL=sum(P.*l);%计算平均码字长度H=entropy(P,2)；%计算信源炳q=H/L;%计算编码效率%打印输出结果fori=l:nend[njn]=size(W);TEMP=32*ones(n,5);W=[W.TEMP];W=W';[njn]=size(W);W=reshape(W,1,n*m);W=spnntR'%s',W);sO=，很好！输入正确，编码结果如下：\sl=Tano编码所得码字W「；s2=Tano编码平均码字长度L「；s3-Fano编码的编码效率q:1；disp(sO);disp(s1).disp(B),disp(W);disp(s2).disp(L);disp(s3):disp(q);附录EHuffman编码Huffman编码(1)TOC\o"1-5"\h\z%huffmaii编码生成器％%函数说明：％%[W,L,q]=huffman(P)为huffman编码函数％%P为信源的概率矢量，W为编码返回的码字%%1为编码返回的平均码字长度，q为编码效率%%*****************************************%fiinction[W,L,q]=hufIinan(P)if(length(fiiid(P<=0))~=0)error(TSIotapiob.vector,negativecomponent*);%判断是否符合概率分布条件endif(abs(sum(P)-1)>10e-10)error(TSIotaprob.vector,componentdonotadduptoI1)%判断是否符合概率和为1endn=length(P);%计算输入元素个数P=P；maik=zeros(n-l,n);%mark为nT行、n列矩阵，用来记录每行最小两概率叠加后概率排列次序%1)确定概率大小值的排列，得到mark矩阵。fori=l:n-l[p,num]=soit(p);%对输入元素排序并纪录mark(i,:)=[num(1:nT+1),zeros(1,i—1)];p=[p⑴+p(2),p(3:n),l];end%2)生成一个n-1行、nl(iiXn)列矩阵table,每行可看做11个段，%每段长为11,记录一个码字(每个码字的长度不会超过n)ofori=l:n-ltable。,:)=blanks(n*n);end%3)计算各个元素码字，循环n-2次，决定矩阵table%从倒数第二行开始到第-行的每段的码字值，到编码表格tabletable(nT,n)=T;%小值赋1table(nT,2*n)=。;%大值赋0fori=2:n-ltable(n-i,l:n-1)=table(n-i+1,n*(fhid(mark(n-i+1,:)==1))…-(n-2):n*(find(mark(n-i+l/)=!)));%按mark的记录依次赋值table(n-i,n)=,r；table(n-i,n+l:2*n-l)=table(n-i,l:nT);table(n-i,2*n)=O;fortable(n-i,(j+l)*n+l：(j+2)*n)=table(n-i+1,..・n*(find(mark(n-i+1,：)==j+1)"!)+!:n*fiiid(niark(n-i+1,：)=j+1))；%按mark的记录依次赋值end%4)得到编码后的码字fori=l:nW(i,1:n)=table(l4i*(fiiid(maik(1,:)==i)-1)+1:find(mark(1,:)=i)*n);l(i)=lengtli(fuid(abs(W(i,:))~=32));endL=sum(P.*l);%计算平均码字长度H=entropy(P,2)；%计算信源炳q=H/L;%计算编码效率%打印输出结果fori=l:nend[m,n]=size(W);TEMP=blaiiks(m);W=[W?TEMP\TEMP\TEMP*];[m,n]=size(W);W=reshape(W,1s0=，很好！输入正确，编码结果如下：,；sl^Huffinan编码所得码字W:1;s2=,Huffinan编码平均码字长度L「；s3=,Huffinan编码的编码效率q:\disp(sO);disp(sl).disp(B).disp(W);disp(s2).disp(L);disp(s3),disp(q);Huffman编码(2)%liuffmaii编码生成器%函数说明：%[W,L,V,q]=huffmaii_better(P)为liuffinan_better编码函数%P为信源的概率矢量，W为编码返回的码字，V为码字的方差％%1为编码返回的平均码字长度，q为编码效率％%*******************************************************%fiinction[W,L,V,q]=huffinaii_better(P)if(length(fiiid(P<=0))~=0)error(TSIotapiob.vector.negativecomponent*);%判断是否符合概率分布条件endif(abs(sum(P)-1)>10e-10)error(TSIotaprob.vector,componentdonotadduptof)%判断是否符合概率和为1endn=length(P);%计算输入元素个数P=P;maik=zeros(n-l,n);%mark为nT行、n列矩阵，用来记录每行最小两概率叠加后概率排列次序%1)确定概率大小值的排列，得到mark矩阵。L1；fbr[p,num]=soit(p);%对输入元素排序并纪录if(广=1)if(count~=0)k=max(a(t,:));fbrs=count：-l:lnum(k-s)=num(k-s+l);%若有与新项相等的项，贝V将新项尽可能排列在其右侧endnum(k)=l;endt=t+l;endmark(i,:)=[num(1:nT+l),zeros(1,i-1)];p=[p⑴+p(2),p(3:n),l];%前两项求和得新项count=0;%用于计数fbrj=2:n-iif(p(i)=po))%判断P中是否有与求和后的新项相等的项count=count+l;a(t,count)可；endendend%2)生成一个n-1行、nl(iiXn)列矩阵table,每行可看做n个段，%每段长为ii,记录一个码字(每个码字的长度不会超过以)。fbrtable(i,:)=blanks(n*n);end%3)计算各个元素码字，循环n-2次，决定矩阵table%从倒数第二行开始到第一行的每段的码字值，到编码表格tabletable(n-14i)~T;%小值赋1table(n-l,2*n)=,0,;%大值赋0fbri=2:nTtable(n-i,l:n-l)=table(nT+1,n*(fkid(mark(nT+1,：)==1))…-(n-2):n*(find(niark(ii-i+1,:)==1)));%按mark的i己录依次赋值table(n-i.n)-T;table(n-i,n+l:2*n-l)=table(n-i,l:n-l);table(n-i,2*n)=O;table(n-i,(j+1)*n+1:(j+2)*n)=table(n-i+1,…n*(find(mark(n-i+l5:)=^j+l)-l)+l:n*fiiid(mark(n-i+l,:)==j+l));%按mark的记录依次赋值endend%4)得到编码后的码字fori=l:nW(i,1:n)=table(l,n*(fiiid(maik(1,:)==i)T)+1:find(mark(1,:)=i)*n);l(i)=length(find(abs(W(i,:))〜=32));endL=sum(P.*l);%计算平均码字长度H=entropy(P,2)；%计算信源炳V=sum(P.*((l-L)."2));%计算码字的方差，以判断编码方法的优劣q=H，L;%计算编码效率%打印输出结果fori=l:nend[m,n]=size(W);TEMP=blaiiks(m);W=[W.TEMP\TEMP\TEMP,];[m,n]=size(W);W=reshape(W,1s0=，很好！输入正确，编码结果如下：，；sl=Huffinan编码所得码字W:1;s2=,Huffinan编码平均码字长度L「；s3=Huffinan编码所得码字W的方差V,s4=rHuffiiian编码的编码效率q:\disp(sO);disp(s1).disp(B),disp(W);disp(s2),disp(L);disp(s3),disp(V);disp(s4).disp(q)附录F信息率失真函数的迭代计算%信息率失真函数的迭代计算，迭代精度取为10-(-7)%%在信源的输入概率分布Pa和失真矩阵d己知的条件下求出信息率失真函数%%函数说明：%%[Pba,Riiuii,Dmax,Sniax]=RateDF(Pa,d,S)为信息率失真函数%%变量说明：%%Pa：信源的输入概率矩阵，d：失真矩阵，S：拉氏乘子%%Pba：最佳正向转移概率矩阵，Smax：最大拉氏乘子%%Riiiui：最小信息率，Dmax：允许的最大失真度%%Pb：信源的输出概率矩阵，D：允许的失真度，R：信息率%%r：输入信源数，s：输出信源数%%*********************************************************************%function[Pba,Rimn.Dmax,Smax]=RateDF(Pa.d.S)%提示错误信息[i\s]=size(d);if(lengtli(find(Pa<=0))~=0)erro^TJotapiob.vector,shoudbepositivecomponent!1);%判断是否符合概率分布条件endif(abs(sum(Pa)-1)>1Oe-10)errorfNotapiob.vector,componentdonotaddupto1「)％判断是否符合概率和为1endif(r~=length(Pa))euor(Theparametersdonotmatch!1);%判断参数是否-,致end%第一步pba=[];RS=[];%R(S)函数初始化DS=[];%D(S)函数初始化m=l;%m为S循环的次数while(1)%外层循环，对S的循环Pba(1:i\1:s,1)=1/s*ones(r,s);%求信道正向转移矩阵Pba%第二步fbrj=l:sPb(j,l)=O；fbri=l:iPb(j,l)=Pb(j,l)+Pa(i)*Pba(iJ,l);%求信源的输出概率矩阵，即Pb(j,l)endendfbr1=1:rtemp(i)=0;fbrj=l:stemp(i)=temp(i)+Pb(j,1)*exp(S(m)*d(iJ));%temp为临时项，求Pba(i,j,2)时表达式的分母endendfori=l:rfbrj=l:sPba(iJ,2)=(Pb(j,l)*exp(S(m)*d(ij)))/ten)p①；endendD⑴=0；fbri=l:rforj=l:sD(l)=D(l)+Pa(i)*Pba(ij,l)*d(iJ);%求D(l)endendR⑴=o；fbri=l:rfbrj=l:sif(Pba(i,j,l)~=0)R(l)-R(l)+Pa(i)*Pba(ij,l)*log2(Pba(ij,l)/Pb(j,l))；%求R(l)endendendn=2;%n为内层循环次数while⑴％内层循环，对精度的循环%第三步forj=l:sPb(j,n)=0；fbri=l:rPb(j,n)=Pb(j,ii)+Pa(i)*Pba(i,j,n);%求输出的信源概率分布endendfbri=l:itemp(i)=0;fbrj=l:stemp(i)=temp(i)+Pb(j,n)*exp(S(m)*d(i,j));%temp为临时项，求Pba(i,j,n+1)时表达式的分母endendfbri=l:ifbrj=l:sif(temp(i)~=0)Pba(iJ,n+l)=(Pb(j.n)*exp(S(m)*d(ij)))/temp①；%求Pba(i,j,n+1)endendend%第四步D(n)=O；fbri=l:ifbrj=l:sD(n)=D(n)+Pa(i)*Pba(iJ,n)*d(i,j);%求D(n)endendR(n)=0；fbri=l:iforj=l:sif(Pba(i,j,n)~=0)R(n)=R(n)+Pa(i)*Pba(i,j,H)*log2(Pba(i,j,n)/Pb(j,n));%求R(n)endendend%判断差别是否在允许的精度范围之内if(abs(R(n)-R(n-!))<=!0A(-7))%R(n)精度判断if(abs(D(n)-D(n-1))<=1(T(-7))%D(n)精度判断break;endendn=n+l;%内层循环次数加1end%内层循环结束%第五步S(m+1)=S(m)+0.5;%第六步if(abs(R(n)<10A(-7)))if(m<=10)disp(此时S的值为:)disp(S(m));errors初始拉氏乘子S取得大了，请取小些!)；%判断拉氏乘子S的初始值是否合适else[k,l,q]=size(pba);Pba=pba(:,:,q);Rmin=min(RS);Dmax=max(DS);Smax=S(mT);break;endendpba=[Pba(:,RS=[RSR(n)];DS=[DSD(n)];m=m+l;%外层循环次数加1end%外层循环结束%打印输出结果s0=，很好！输入正确，迭代结果如下二sl=，最佳转移概率分布Pba:1;s2=，最小信息率Rimn/;s3=，最大失真度Dmax:*;s4=，最大拉氏乘子Smax/;disp(sO);disp(sl),disp(Pba);disp(s2).disp(Riiiui);disp(s3).disp(Dniax);disp(s4).disp(Sniax);%画出信息率失真函数R(D)plot(DS,RS)xlabelf允许的失真度D)ylabel('信息率失真函数R(D)')title(，信息率失真函数R(D)的曲线图)1234567891011121314151617181920212223