《信息论matlab的实验用到的编码(共17页).doc》由会员分享，可在线阅读，更多相关《信息论matlab的实验用到的编码(共17页).doc（17页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上 Matlab源程序附录A 信息熵% 函数说明： % H=entropy(P,r) 为信息熵函数 % P为信源的概率矢量, r为进制数 % H为信息熵 %* %function H=entropy(P,r)if (length(find(P10e-10) error(Not a prob.vector,component do not add up to 1); end H=(sum(-P.*log2(P)/(log2(r)+eps);附录B 离散无记忆信道容量的迭代计算% 信道容量C的迭代算法 % 函数说明： % CC,Paa=ChannelCap(P,k) 为信道
2、容量函数 % 变量说明： % P：输入的正向转移概率矩阵，k：迭代计算精度 % CC：最佳信道容量，Paa：最佳输入概率矩阵 % Pa：初始输入概率矩阵，Pba：正向转移概率矩阵 % Pb：输出概率矩阵，Pab：反向转移概率矩阵 % C：初始信道容量， r：输入符号数，s：输出符号数 %*%function CC,Paa=ChannelCap(P,k)% 提示错误信息 if (length(find(P10e-10) error(Not a prob.vector,component do not add up to 1) % 判断是否符合概率和为1end % 1)初始化Par,s=size(
3、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); endend% (2)再求Pabsuma=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
4、(j,i)+eps)/(Pa(i)+eps);end% 3)求信道容量CC=sum(suma);% 4)求下一次Pa,即PaaL=zeros(1,r);sumaa=0;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); endsumaa=sum(a);for i=1:r Paa(i)=a(i)/(sumaa+eps); end% 5)求下一次C,即CCCC=log2(sumaa);Pa=Paa;end% 打印输出结果s0=很好！输入正确，迭代结果如下：;s1=最佳输入概率分布Pa:;s2=信
5、道容量C:;s3=迭代次数n:;s4=输入符号数r:;s5=输出符号数s:;s6=迭代计算精度k:;for i=1:r Bi=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);附录C Shannon编码% 函数说明： % p,x=array(P,X) 为按降序排序的函数 % % P为信源的概率矢量，X为概率元素的下标矢量 % p为排序后返回的信源的概率矢量 % x为排序后返回的概率元素的下
6、标矢量 %*%function p,x=array(P,X)P=P;X;l,n=size(P);for i=1:n max=P(1,i); maxN=i; MAX=P(:,i); for j=i:n if (max1) if (in) for k=(maxN-1):-1:i P(:, k+1)=P(:,k); end end end P(:,i)=MAX;endp=P(1,:);x=P(2,:);% shannon编码生成器 % 函数说明： % W,L,q=shannon(p) 为shannon编码函数 % % p为信源的概率矢量，W为编码返回的码字% L为编码返回的平均码字长度，q为编码效率
7、%*%function W,L,q=shannon(p)% 提示错误信息 if (length(find(p10e-10) error(Not a prob.vector,component do not add up to 1) % 判断是否符合概率和为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);for i=2:n P(i)=P(i-1)+p(i-1);end% 4)求得二进制代码组W% a)将十进制数转为二进制数for i=1:
8、n for j=1:l(i) temp(i,j)=floor(P(i)*2); P(i)=P(i)*2-temp(i,j); endend% b)给W赋ASCII码值，用于显示二进制代码组Wfor i=1:n for j=1:l(i) if (temp(i,j)=0) W(i,j)=48; else W(i,j)=49; end endendL=sum(p.*l); % 计算平均码字长度H=entropy(p,2); % 计算信源熵q=H/L; % 计算编码效率% 打印输出结果for i=1:n Bi=i;endn,m=size(W);TEMP=32*ones(n,6);W=W,TEMP;W=
9、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);附录D Fano编码% 函数说明： % next_P,next_index,code_num=compare(current_P,current_index) % %为比较函数，主要用于信源符号的分组
10、 % current_P为当前分组的信源的概率矢量 % current_index为当前分组的信源的下标 % next_P为返回的下次分组的信源的概率矢量 % 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)求概率的依次累加和for i=2:n add(i)=0; add(i)=add(i-1)+current_P(i);
11、end% 2)求概率和最接近的两小组s=add(n);for i=1:n temp(i)=abs(s-2*add(i);endc,k=min(temp);% 3)对分组的信源赋ASCII值if (current_index=k) next_index=current_index; code_num=48; next_P=current_P(1:k);else next_index=current_index-k; code_num=49; next_P=current_P(k+1):n);end% fano编码生成器 % 函数说明： % W,L,q=fano(P) 为fano编码函数 % %
12、P为信源的概率矢量，W为编码返回的码字 % L为编码返回的平均码字长度，q为编码效率 %*%function W,L,q=fano(P)% 提示错误信息 if (length(find(P10e-10) error(Not a prob.vector,component do not add up to 1) % 判断是否符合概率和为1end % 1)排序n=length(P); x=1:n;P,x=array(P,x);% 2)将信源符号分组并得到对应的码字for i=1:n current_index=i; j=1; current_P=P; while 1 next_P,code_num
13、,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; end end l(i)=length(find(abs(W(i,:) =0); % 得到各码字的长度endL=sum(P.*l); % 计算平均码字长度H=entropy(P,2); % 计算信源熵q=H/L; % 计算编码效率% 打印输出结果for i=1:n Bi=i;endn,m=size(W);T
14、EMP=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); 附录E Huffman编码Huffman编码（1）% huffman编码生成器 % 函数说明： % W,L,q=huffman(P) 为huffman编码函数
15、% % P为信源的概率矢量，W为编码返回的码字 % L为编码返回的平均码字长度，q为编码效率%*%function W,L,q=huffman(P) if (length(find(P10e-10) error(Not a prob.vector,component do not add up to 1) % 判断是否符合概率和为1end n=length(P); % 计算输入元素个数p=P; mark=zeros(n-1,n); % mark为n-1行、n列矩阵，用来记录每行最小两概率叠加后概率排列次序% 1) 确定概率大小值的排列，得到mark矩阵。for i=1:n-1 p,num=so
16、rt(p); % 对输入元素排序并纪录 mark(i,:)=num(1:n-i+1),zeros(1,i-1); p=p(1)+p(2),p(3:n),1; end % 2) 生成一个n-1行、n1(nn)列矩阵table，每行可看做n个段，% 每段长为n，记录一个码字（每个码字的长度不会超过n）。 for i=1:n-1 table(i,:)=blanks(n*n); end % 3) 计算各个元素码字，循环n-2次，决定矩阵table% 从倒数第二行开始到第一行的每段的码字值，到编码表格tabletable(n-1,n)=1; % 小值赋1table(n-1,2*n)=0; % 大值赋0
17、for i=2:n-1 table(n-i,1:n-1)=table(n-i+1,n*(find(mark(n-i+1,:)=1). -(n-2):n*(find(mark(n-i+1,:)=1); % 按mark的记录依次赋值 table(n-i,n)=1; table(n-i,n+1:2*n-1)=table(n-i,1:n-1); table(n-i,2*n)=0; for j=1:i-1 table(n-i,(j+1)*n+1:(j+2)*n)=table(n-i+1,. n*(find(mark(n-i+1,:)=j+1)-1)+1:n*find(mark(n-i+1,:)=j+1)
18、; % 按mark的记录依次赋值end % 4)得到编码后的码字for i=1:n W(i,1:n)=table(1,n*(find(mark(1,:)=i)-1)+1:find(mark(1,:)=i)*n); l(i)=length(find(abs(W(i,:) =32); end L=sum(P.*l); % 计算平均码字长度H=entropy(P,2); % 计算信源熵q=H/L; % 计算编码效率% 打印输出结果for i=1:n Bi=i;endm,n=size(W);TEMP=blanks(m);W=W,TEMP,TEMP,TEMP;m,n=size(W);W=reshape(
19、W,1,m*n);s0=很好！输入正确，编码结果如下：;s1=Huffman编码所得码字W:;s2=Huffman编码平均码字长度L:;s3=Huffman编码的编码效率q:;disp(s0);disp(s1),disp(B),disp(W);disp(s2),disp(L);disp(s3),disp(q); Huffman编码（2）% huffman编码生成器 % 函数说明： % W,L,V,q=huffman_better(P) 为huffman_better编码函数 % % P为信源的概率矢量，W为编码返回的码字，V为码字的方差% L为编码返回的平均码字长度，q为编码效率 %*%fun
20、ction W,L,V,q=huffman_better(P) if (length(find(P10e-10) error(Not a prob.vector,component do not add up to 1) % 判断是否符合概率和为1end n=length(P); % 计算输入元素个数p=P; mark=zeros(n-1,n); % mark为n-1行、n列矩阵，用来记录每行最小两概率叠加后概率排列次序% 1) 确定概率大小值的排列，得到mark矩阵。t=1;for i=1:n-1 p,num=sort(p); % 对输入元素排序并纪录 if (i=1) if (count=
21、0) k=max(a(t,:); for s=count:-1:1 num(k-s)= num(k-s+1); % 若有与新项相等的项，则将新项尽可能排列在其右侧 end num(k)=1; end t=t+1; end mark(i,:)=num(1:n-i+1),zeros(1,i-1); p=p(1)+p(2),p(3:n),1; % 前两项求和得新项 count=0; % 用于计数 for j=2:n-i if (p(1)=p(j) % 判断p中是否有与求和后的新项相等的项 count=count+1; a(t,count)=j; end endend % 2) 生成一个n-1行、n1
22、(nn)列矩阵table，每行可看做n个段，% 每段长为n，记录一个码字（每个码字的长度不会超过n）。 for i=1:n-1 table(i,:)=blanks(n*n); end % 3) 计算各个元素码字，循环n-2次，决定矩阵table% 从倒数第二行开始到第一行的每段的码字值，到编码表格tabletable(n-1,n)=1; % 小值赋1table(n-1,2*n)=0; % 大值赋0 for i=2:n-1 table(n-i,1:n-1)=table(n-i+1,n*(find(mark(n-i+1,:)=1). -(n-2):n*(find(mark(n-i+1,:)=1);
23、 % 按mark的记录依次赋值 table(n-i,n)=1; table(n-i,n+1:2*n-1)=table(n-i,1:n-1); table(n-i,2*n)=0; for j=1:i-1 table(n-i,(j+1)*n+1:(j+2)*n)=table(n-i+1,. n*(find(mark(n-i+1,:)=j+1)-1)+1:n*find(mark(n-i+1,:)=j+1); % 按mark的记录依次赋值 end end % 4)得到编码后的码字for i=1:n W(i,1:n)=table(1,n*(find(mark(1,:)=i)-1)+1:find(mark
24、(1,:)=i)*n); l(i)=length(find(abs(W(i,:)=32); end L=sum(P.*l); % 计算平均码字长度H=entropy(P,2); % 计算信源熵V=sum(P.*(l-L).2); % 计算码字的方差，以判断编码方法的优劣q=H/L; % 计算编码效率% 打印输出结果for i=1:n Bi=i;endm,n=size(W);TEMP=blanks(m);W=W,TEMP,TEMP,TEMP;m,n=size(W);W=reshape(W,1,m*n);s0=很好！输入正确，编码结果如下：;s1=Huffman编码所得码字W:;s2=Huffma
25、n编码平均码字长度L:;s3=Huffman编码所得码字W的方差V:;s4=Huffman编码的编码效率q:;disp(s0);disp(s1),disp(B),disp(W);disp(s2),disp(L);disp(s3),disp(V); disp(s4),disp(q)附录F 信息率失真函数的迭代计算% 信息率失真函数的迭代计算，迭代精度取为10(-7) % 在信源的输入概率分布Pa和失真矩阵d已知的条件下求出信息率失真函数 % 函数说明： % Pba,Rmin,Dmax,Smax=RateDF(Pa,d,S) 为信息率失真函数 % 变量说明： % Pa：信源的输入概率矩阵，d：失真
26、矩阵，S：拉氏乘子 % % Pba：最佳正向转移概率矩阵， Smax：最大拉氏乘子 % % Rmin：最小信息率，Dmax：允许的最大失真度 % % Pb：信源的输出概率矩阵，D：允许的失真度，R：信息率% r：输入信源数，s：输出信源数 %*%function Pba,Rmin,Dmax,Smax=RateDF(Pa,d,S)% 提示错误信息 r,s=size(d);if (length(find(Pa10e-10) error(Not a prob.vector,component do not add up to 1!) % 判断是否符合概率和为1end if (r=length(Pa)
27、 error(The parameters do not match!); % 判断参数是否一致end% 第一步pba=;RS=; % R(S)函数初始化DS=; % D(S)函数初始化m=1; % m为S循环的次数while (1) % 外层循环，对 S 的循环 Pba(1:r,1:s,1)=1/s*ones(r,s); % 求信道正向转移矩阵Pba % 第二步 for j=1:s Pb(j,1)=0; for i=1:r Pb(j,1)=Pb(j,1)+Pa(i)*Pba(i,j,1);% 求信源的输出概率矩阵,即Pb(j,1) end end for i=1:r temp(i)=0; f
28、or j=1:s temp(i)=temp(i)+Pb(j,1)*exp(S(m)*d(i,j);% temp为临时项,求Pba(i,j,2)时表达式的分母 endend for i=1:r for j=1:s Pba(i,j,2)=(Pb(j,1)*exp(S(m)*d(i,j)/temp(i); end end D(1)=0; for i=1:r for j=1:s D(1)=D(1)+Pa(i)*Pba(i,j,1)*d(i,j); % 求D(1) end end R(1)=0; for i=1:r for j=1:s if (Pba(i,j,1) =0) R(1)=R(1)+Pa(i)
29、*Pba(i,j,1)*log2(Pba(i,j,1)/Pb(j,1); % 求R(1) end end end n=2; % n为内层循环次数 while (1) % 内层循环，对精度的循环% 第三步 for j=1:s Pb(j,n)=0; for i=1:r Pb(j,n)=Pb(j,n)+Pa(i)*Pba(i,j,n); % 求输出的信源概率分布 end end for i=1:r temp(i)=0; for j=1:s temp(i)=temp(i)+Pb(j,n)*exp(S(m)*d(i,j);% temp为临时项,求Pba(i,j,n+1)时表达式的分母 end end f
30、or i=1:r for j=1:s if (temp(i) =0) Pba(i,j,n+1)=(Pb(j,n)*exp(S(m)*d(i,j)/temp(i); % 求Pba(i,j,n+1) end end end% 第四步 D(n)=0; for i=1:r for j=1:s D(n)=D(n)+Pa(i)*Pba(i,j,n)*d(i,j); % 求D(n) end end R(n)=0; for i=1:r for j=1:s if (Pba(i,j,n) =0) R(n)=R(n)+Pa(i)*Pba(i,j,n)*log2(Pba(i,j,n)/Pb(j,n);% 求R(n)
31、end end end% 判断差别是否在允许的精度范围之内 if (abs(R(n)-R(n-1)=10(-7) % R(n)精度判断 if (abs(D(n)-D(n-1)=10(-7) % D(n)精度判断 break; end end n=n+1; % 内层循环次数加1 end % 内层循环结束% 第五步 S(m+1)=S(m)+0.5;% 第六步 if (abs(R(n)10(-7) if (m=10) disp(此时S的值为:),disp(S(m); error(初始拉氏乘子S取得大了，请取小些!); % 判断拉氏乘子S的初始值是否合适 else k,l,q=size(pba); Pba=pba(:,:,q); Rmin=min(RS); Dmax=max(DS); Smax=S(m-1); break; end end pba=Pba(:,:,:);