huffman.m

现代通信系统（Proakis）第二版

function [h,l]=huffman(p);
%HUFFMAN 	Huffman code generator
%	 	[h,l]=huffman(p), Huffman code generator
%	 	returns h the Huffman code matrix, and l the
%        	average codeword length for a source with 
%	 	probability vector p. 

%[h,l]=huffman(p),输入为一维行矩阵p，p为各符号的概率分布，概率和为1，各元素值为
%正，输出H矩阵为对应每个符号概率的码字，L为输出码字的平均码长。Huffman .m运用典
%型的IF和FOR控制流循环语句，该程序包括两个IF 控制流和5个FOR 循环结构。
%第一个IF 语句判断输入P矩阵各元素是否全为大于零的有效概率值；第二个IF 语句判断
%输入矩阵的概率和是否为合理值1。  
%N取输入行向量P的长度，即需要编码元素个数。  
%M为N-1行、N列矩阵，用来记录每行最小两概率叠加后概率排列次序。
%第一个FOR 循环确定概率大小值的排列，得到M矩阵。
%第二个FOR循环生成一个N-1行、N2（N×N）列矩阵C，每行可看作N个段，每段长为N，记
%录一个码字（每个码字的长度不会超过N）。
%给C矩阵的N-1行的第一个段赋值0，第二个段赋值1，这两个码字对应编码中最后相加为一
%的两个概率。
%第三个循环是本程序的主要部分，循环N-2次，决定矩阵C从倒数第二行开始到第一行的每
%段的码字值。每一行值都从下一行值得到，找到在下一行码字中相加本行最小两个概率得
%到的概率的对应码字，本行两个最小概率对应码字分别为此码字最后加“0”，加“1”。
%嵌套的第四个FOR循环找到其余的本行在下一行对应的码字，该码字保持不变。循环结束
%后，C矩阵第一行的N段对应输入N个概率所对应符号的码字。该码字按码字长短排列。
%第五个FOR循环根据M矩阵第一行记录的概率排序位置分配给每个概率对应符号的码字。  
%FOR EXAMPLE：  
%P=[1/6,1/4,5/12,1/6];  
%N=4;  
%M矩阵：  
%     m =[ 1 4 2 3；2 1 3 0；2 1 0 0]  



%N矩阵：  

%     n =[1110     1110          110       0；  
%           10       11            0        ；  
%            0        1                      ]  



if length(find(p<0))~=0,
  error('Not a prob. vector, negative component(s)')
end
if abs(sum(p)-1)>10e-10,
  error('Not a prob. vector, components do not add up to 1')
end
n=length(p);
q=p;
m=zeros(n-1,n);
for i=1:n-1
  [q,l]=sort(q);
  m(i,:)=[l(1:n-i+1),zeros(1,i-1)];
  q=[q(1)+q(2),q(3:n),1];
end
for i=1:n-1
  c(i,:)=blanks(n*n);
end
c(n-1,n)='0';
c(n-1,2*n)='1';
for i=2:n-1
  c(n-i,1:n-1)=c(n-i+1,n*(find(m(n-i+1,:)==1))...
  -(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);
  l1(i)=length(find(abs(h(i,:))~=32));
end
l=sum(p.*l1);