huffman2.m

离散信源无失真信源编码（实现huffman编码） 可用于教学

%function [h,l]=huffman(p)                       
%功能：实现霍夫曼编码
%问题：1、霍夫曼编码的步骤？2、概率的顺序排列在程序中是怎么实现的？排序后返回的q、l值是什么含义？(程序的第27行)3、霍夫曼编码编出的码是不是唯一的(不是)？哪些因素可以
%使霍夫曼编码编出的码不同，但都是正确的？(概率由大到小排列；之后两个概率合并时，先标1后标0还是先0后1；合并的概率是否排在最上面）4、本程序是怎么得到的霍夫曼码（流程），确定可使霍夫曼编码不同的不确定因素是什么？5、如果我想编出的码
%正好和本程序的码相反（即0、1互换），则应该怎么修改程序？6、本程序信源的熵和编码效率是怎么求的？7、霍夫曼码不是唯一的，平均码长，编码效率是不是唯一的呢？（唯一）

                                                 %h是霍夫曼编完的码{111  10  0  1101  1100 }“该程序编码的方法：概率由大到小排列；之后两个概率合并时，先标1后标0（如果先标0后标1则应该怎么该？把36、37、42、44行语句中的0、1取反即可）；合并的概率尽量排在最上面；倒写”
                                                 % l是在第19行对概率升序排列之后，对应原来已知概
                                                 % 率的顺序，有了它，最后你编完的码与原来的信源好有
                                                 % 个对应关系啊，用它再来对应上（程序已经对应完）
                                                 %p是信源的概率[0.2,0.3,0.4,0.06,0.04];
p=[0.4,0.3,0.2,0.06,0.04];                       %信源的概率  改为这个非常好：p=[0.2,0.3,0.4,0.06,0.04]; 
N=4;  

if (length(find(p<0))~=0)                        %“~”表示Logical NOT；“~=”就是不等于                        
   error('Not a prob,negative component');       %p数组中概率有负的不行，拒绝编码（组成）
end  
if (abs(sum(p)-1)>10e-10)  
   error('Not a prob.vector,component do not add to 1') %如果不是一个完备集 所有概率之和与1差的绝对值大于10的负10次方的话就错误，不算了 probability vector  概率向量
end  

n=length(p);                                     %数组中元素的个数  n=5，后面循环操作用
q=p;                                             %信源的概率赋给q，以便后面排序、合并用，不用p是因为后面算平均码长时还得用原来的p这个概率
m=zeros(n-1,n);                                  %定义m为4行5列的零数组，以便后面记录排列顺序（小概率求和后在重新排列，这就存在排列顺序的问题）

for i=1:n-1                                      %i为1到4
   [q,l]=sort(q);                                %l是在这定义的  l是排列后的顺序 sort(q)是对q进行升序【sort（挑选） ascending（上升）】排列（就是实际中做题时概率递减），l为排完之后原来的顺序！ 第一次排完之后l为 5    4     1    2    3
   %l
   %q                                            %按照升序排列之后概率也要交换顺序，变后的为新的q
   m(i,:)=[l(1:n-i+1),zeros(1,i-1)];             %概率小的相加之后重新排列的顺序，所有的中间排列顺序都记录在m矩阵中
   q=[q(1)+q(2),q(3:n),1];                       %最后两个概率小的合并成一个 
   
end  
for i=1:n-1  
   c(i,:)=blanks(n*n);                           %c(i,:)   blanks(n)返回n个零或空格的字符串？blanks(n*n)是一个1行n*n列的数组，里面是空的
   %c0=c
   %blanks(n*n)                                  c(i,:)的每一行都是blanks(n*n)  i为1到4，4次循环，所以c为4*25的空矩阵
end  
c(n-1,n)='0';                                    %把c这个空数组的4行5列这个位置赋值0
c(n-1,2*n)='1';                                  %把c这个空数组的4行10列这个位置赋值1
%c2=c
for i=2:n-1  
   c(n-i,1:n-1)=c(n-i+1,n*(find(m(n-i+1,:)==1))...  %接续行  c(n-i,1:n-1)   find()查找字符串中的....函数
                -(n-2):n*(find(m(n-i+1,:)==1))); 
   c(n-i,n)='0';  
   c(n-i,n+1:2*n-1)=c(n-i,1:n-1);  
   c(n-i,2*n)='1';  
   for j=1:i-1  
      c(n-i,(j+1)*n+1:(j+2)*n)=c(n-i+1,...  
         n*(find(m(n-i+1,:)==j+1)-1)+1:n*find(m(n-i+1,:)==j+1));  
   end  
end   
 
   for i=1:n  
      h(i,1:n)=c(1,n*(find(m(1,:)==i)-1)+1:find(m(1,:)==i)*n);    %h是在这定义的，是霍夫曼编完的码字
      ll(i)=length(find(abs(h(i,:))~=32));                        %ll是在这定义的，含义为各码的码长
   end  
   h                                                              %怎么用图形的方式显示数组呢？？？？？？？？？
   la=sum(p.*ll)                                                  %la是平均码长
   r=2;
H=(sum(-p.*log2(p)))/(log2(r)+eps)                                %MATLAB定义的正的极小值=2.2204e-16,就跟PI定义为3.1415926   浮点精确？？？？？？？？！【直接这么写行不行？H=sum(-p.*log2(p))】
RW=H/la                                                           %编码效率