哈夫曼编码的MATLAB实现

在手动计算时，对哈夫曼编码的流程还是非常熟悉的，所以难点就是代码实现。在实现以前，我考虑了以下一些问题：一个符号的哈夫曼编码值可以用一个向量进行存储，符号合并以后的结果用什么表示，对其赋值0,1如何作用在被合并的原始符号上面。由于哈夫曼编码是从树根开始，存在继承上一节点编码的现象，故要与手动计算的方向相反，先根据概率及其合并得出各层次之间的关系，然后从底层向上开始编码。
编码中用到的一个核心函数是sort函数，用法如下：
[B,I]=sort(A)
如果 A=[5 4 7 0 2 ]
B=[0 2 4 5 7]
I=[4 5 2 1 3]是排序以后的数在原序列中的位置
以输入概率P=[0.2 0.15 0.4 0.25]为例对编码思路进行分析。
每经过一次合并，排序的概率少一个，最终由n个变成2个。为此可开辟一个n-1行 n*n列的编码矩阵c，第一行是对n个符号进行编码，直到n-1行是对2个符号进行编码。上一层编码是在下一层的基础上增加0或1 但是是在下一层哪个编码的基础上增加呢，这就要知道概率和是由原来的哪两个概率就和得到的。//为解决这个问题，需要用到sort()函数返回I的用法（上文提到的）。构建n-1行n列的合并关系矩阵m，与编码矩阵相对应。m矩阵对应位置的值就是指编码矩阵中该位置在排序前是在哪一位。
由于哈夫曼编码的不唯一性，规定概率小编码0，大的编码1。
生成m矩阵
红框中3的含义是：四个符号按概率从小到大进行排序，将概率较小的两个合并，变成三个符号。这三个符号进行新的排序以后，该位置上的符号是原始序列的第三个。

下面是代码及其详细解释:

p=[0.2 0.15 0.4 0.25]
%判断信源的合理性
if length(find(p<0))~=0
error('错误！出现小于0的概率')
end

if abs(sum(p)-1)>10e-10
error('错误！信源不满足完备性')
end
n=length(p);
q=p;
%构建概率合并关系矩阵
m=zeros(n-1,n);
for i=1:n-1%合并n-1次
[q,e]=sort(q);%E是个行向量 值代表q中该位置的概率在原序列中的位置
m(i,:)=[e(1:n-i+1),zeros(1,i-1)];
%m的特点：第一行n个非零元素，逐行递减，n-1行2个非零元素
q=[q(1)+q(2),q(3:n),1];
end
%构建编码矩阵
for i=1:n-1
c(i,1:n*n)=blanks(n*n); %c矩阵特点 与m对应  每一个n块对应一个概率的符号 第一行有n个 n-1行有2个 其余为空
end
%开始编码 因为哈夫曼编码不唯一 自规定概率小的赋值0 概率大的赋值1
c(n-1,n)='0';
c(n-1,2*n)='1';
for i=2:n-1
%被合并的前两个需要分配0和1
index=find(m(n-i+1,:)==1);%下一行对应位置为1 则为该行较小的两个概率合并以后的总概率
c(n-i,1:n-1)=c(n-i+1,n*(index)-(n-2):n*index);%前缀是对应的合并后概率的编码值 注意：1：n-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';

%未合并的 从每一行第三个开始编码 无需新分配0 1 直接继承上一个节点的编码
for j=1:i-1
index=find(m(n-i+1,:)==j+1);%从index=2开始到index=i结束 i=3 即n-3行 就是index从2到3
c(n-i,(j+1)*n+1:(j+2)*n)=c(n-i+1,n*(index-1)+1:n*index);
end
end
for i=1:n
index=find(m(1,:)==i);
h(i,1:n)=c(1,n*(index-1)+1:index*n); %将与输入概率对应的编码放入h中
is0=abs(h(i,:));%abs()将字符变为ASCII码 空为32
len(i)=length(find(is0~=32)); %非空符号的个数为码长
end
disp('编码结果');h
disp('平均码长');
e=sum(p.*len) %计算平均码长

总结：最困难的应该是理解概率合并关系矩阵m和编码矩阵c的关系，其实m就像一张地图，告诉c在自下向上编码时，继承的编码应该是哪个节点的。当然，里面涉及到提取与赋值矩阵的部分元素，用到i,j等，要非常清晰边界是什么，要注意赋值时矩阵类型要一致，否则就很容易在细节处出错。