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)=kn=n+1;% (1)先求 Pb Pb=zeros(1,s); for j=1:sfor i=1:rPb(j)=Pb(j)+Pa(i)*Pba(i,j); endend% (2)再求 Pabsuma=zeros(1,s);for j=1:sfor i=1:rPab(j,i)=Pa(i)*Pba(i,j)/(Pb(j)+eps);suma(j)=suma(j)+Pa(i)*Pba(i,j)*log2(Pab(j,
4、i)+eps)/(Pa(i)+eps);end% 3)求信道容量 CC=sum(suma);% 4)求下一次 Pa,即 PaaL=zeros(1,r);sumaa=0;for i=1:rfor j=1:sL(i)=L(i)+Pba(i,j)*log(Pab(j,i)+eps); enda(i)=exp( L(i); endsumaa=sum(a);for i=1:rPaa(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:rBi=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:nmax=P(1,i);maxN=i;MAX=P(:,i);for j=i:nif (max1)if (i10e-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);f
7、or i=2:nP(i)=P(i-1)+p(i-1);end% 4)求得二进制代码组 W% a)将十进制数转为二进制数for i=1:nfor 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:nfor j=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; % 计算编码效率% 打印输
8、出结果for i=1:nBi=i;endn,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);附录 D Fano 编码% 函数说明： % next_P,next_in
9、dex,code_num=compare(current_P,current_index) % % 为比较函数，主要用于信源符号的分组 % 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)=c
10、urrent_P(1);% 1)求概率的依次累加和for i=2:nadd(i)=0;add(i)=add(i-1)+current_P(i);end% 2)求概率和最接近的两小组s=add(n);for i=1:ntemp(i)=abs(s-2*add(i);endc,k=min(temp);% 3)对分组的信源赋 ASCII 值if (current_index10e-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
11、,x);% 2)将信源符号分组并得到对应的码字for i=1:ncurrent_index=i;j=1;current_P=P;while 1next_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); % 计算平均码字长度
12、H=entropy(P,2); % 计算信源熵q=H/L; % 计算编码效率% 打印输出结果for i=1:nBi=i;endn,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); 附
13、录 E Huffman 编码Huffman 编码（ 1）% huffman 编码生成器 % 函数说明： % W,L,q=huffman(P) 为 huffman 编码函数 % % 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=zero
14、s(n-1,n); % mark 为 n-1 行、n 列矩阵，用来记录每行最小两概率叠加后概率排列次序% 1) 确定概率大小值的排列，得到 mark 矩阵。for i=1:n-1 p,num=sort(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(n n)列矩阵 table，每行可看做 n 个段，% 每段长为 n，记录一个码字（每个码字的长度不会超过 n）。 for i=1:n-1 table(i,:)=blanks(n*n); end % 3)
15、 计算各个元素码字，循环 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); % 按 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;
16、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 % 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); %
17、计算信源熵q=H/L; % 计算编码效率% 打印输出结果for i=1:nBi=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=Huffman 编码平均码字长度 L:;s3=Huffman 编码的编码效率 q:;disp(s0);disp(s1),disp(B),disp(W);disp(s2),disp(L);disp(s3),disp(q); Huffman 编码（ 2）% huffman
18、编码生成器 % 函数说明： % W,L,V,q=huffman_better(P) 为 huffman_better 编码函数 % % P 为信源的概率矢量，W 为编码返回的码字， V 为码字的方差 % L 为编码返回的平均码字长度，q 为编码效率 %*%function 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
19、-1,n); % mark 为 n-1 行、n 列矩阵，用来记录每行最小两概率叠加后概率排列次序% 1) 确定概率大小值的排列，得到 mark 矩阵。t=1;for i=1:n-1 p,num=sort(p); % 对输入元素排序并纪录if (i=1)if (count=0) k=max(a(t,:);for s=count:-1:1num(k-s)= num(k-s+1); % 若 有 与 新 项 相 等 的 项 ， 则 将 新 项 尽 可 能 排 列 在 其 右 侧endnum(k)=1;endt=t+1;endmark(i,:)=num(1:n-i+1),zeros(1,i-1); p=
20、p(1)+p(2),p(3:n),1; % 前两项求和得新项count=0; % 用于计数for j=2:n-iif (p(1)=p(j) % 判断 p 中是否有与求和后的新项相等的项count=count+1;a(t,count)=j; endendend % 2) 生成一个 n-1 行、n1(n n)列矩阵 table，每行可看做 n 个段，% 每段长为 n，记录一个码字（每个码字的长度不会超过 n）。 for i=1:n-1 table(i,:)=blanks(n*n); end % 3) 计算各个元素码字，循环 n-2 次，决定矩阵 table% 从倒数第二行开始到第一行的每段的码字值
21、，到编码表格 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); % 按 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,.
