算术编码的原理和MATLAB实现

转自http://blog.sina.com.cn/s/blog_4a8f0cbc01000b6e.html

算术编码与译码原理： 
1、 编码过程
算术编码方法是将被编码的一则消息或符号串（序列）表示成0和1之间的一个间隔（Interval），即对一串符号直接编码成[0,1]区间上的一个浮点小数。符号序列越长，编码表示它的间隔越小，表示这一间隔所需的位数就越多。信源中的符号序列仍然要根据某种模式生成概率的大小来减少间隔。可能出现的符号概率要比不太可能出现的符号减少范围小，因此，只正加较少的比特位。
在传输任何符号串之前，0符号串的完整范围设为[0,1]。当一个符号被处理时，这一范围就依据分配给这一符号的那一范围变窄。算术编码的过程，实际上就是依据信源符号的发生概率对码区间分割的过程。
举例说明如下：
假设一则消息“static_tree”具有如下的概率分布：
             字符                   概率
     ---------------------------------------------------------------
            ＿（space）             0.1
             a                      0.1
             e                      0.3
             r                      0.1
             s                      0.1
             t                      0.3
下面用算术编码方法给该消息编码。
一旦字符的概率已知，就沿着“概率线”为每一个单独的符号设定一个范围，哪一个被设定到哪一段范围并不重要，只要编码和解码都以同样方式进行就可以，这里所用的6个字符被分配的范围（range）如下：
     字符             概率             范围
   ＿(space)          0.1            0≤r<0.1
      a               0.1            0.1≤r<0.2
      e               0.3            0.2≤r<0.5
      r               0.1            0.5≤r<0.6
      s               0.1            0.6≤r<0.7
      t               0.3            0.7≤r<1.0
     ----------------------------------------------------------------
          对“state_tree”的算术编码过程为：
（1） 初始化时，被分割的范围range=high-low=[0,1）,下一个范围的低、高端分别由下式计算：
             Low=low+range×range low
             High=low+range×range high
其中等号右边的low为上一个被编码字符的范围低；range low和range high分别为被编码符号已给定的字符出现概率范围的low和high。
（2） 对消息第一字符s编码：s的range low=0.6，         s的range high=0.7因此，下一个区间的low和high为：
                Low=low+range×range low=0+1×0.6=0.6
                High=low+range×range high=0+1×0.7=0.7
                  Range=high-low=0.7-0.6=0.1
                  S将区间[0,1）=>[0.6,0.7）
      (3)对第二个字符t编码，使用的新生范围为[0.6,0.7），因为t的range low=0.7,range high=1.0,因此下一个low,high分别为
                   Low=0.6+0.1×0.7＝0.67
                   High=0.6+0.1×1.0＝0.70
                   Range=0.7-0.67=0.03
                   t将[0.6,0.7）=>[0.67,0.70）
       (4)对第三个字符a编码，在新生成的[0.67,0.70）中进行分割，因为a的range low=0.10,range high=0.2,因此下一个low,high分别为
                   Low=0.67+0.03×0.1＝0.673
                   High=0.67+0.03×0.2＝0.676
                   Range=0.676-0.673=0.003
                   a将[0.67,0.70）=>[0.673,0.676）
        (5)对第四个字符t编码，在新生成的[0.673,0.676）上进行分割。因为t的range low=0.70,range high=1.0,则下一个low,high分别为
                   Low=0.673+0.003×0.7＝0.6751
                   High=0.673+0.003×1.0＝0.676
                   Range=0.0009
                   t将[0.673,0.676）=>[0.6751,0.676）
         同理得到下面各字符e,＿,s,t,r,e,e编码所得到的范围分别为[0.67528,0.67555），[0.67528,0.675307），[0.675 298 9,0.675 307），[0.675 302 95,0.675 303 76），[0.675 303 112,0.675 303 355），[0.675 303 160 6,0.675 303 233 5）
      将编码后的区间范围综合如下：
          
    上述编码过程可以用下面的区间分割过程表示：


     
2、 解码过程
    解码是编码的逆过程，了解了编码过程后，理解解码过程的操作就相对容易了。通过编码，最后的下界值0.675 303 160 6就是消息“state_tree”的唯一编码。然后解码时，通过判断哪一个符号能拥有我们已编码的消息落在的空间来找出消息中的第一个符号。由于0.675 303 160 6落在[0.6,0.7）之间，因此可解得第一个符号是s。
在解出s后，由于我们知道s的范围的上界和下界，利用编码的逆作用，首先去掉s的下界值0.6，得0.075 303 160 6，然后用s的范围range=0.1去除所得到的0.075 303 160 6，即得到0.753 031 606，接着找出它所在的区间，就是t的原来范围[0.7,1.0）。去掉t的下界值0.7，得0.053 031 606，然后用t的range=0.3除该数得出0.176 772 02，该值所属范围就是字符a……如此操作下去便得到消息的准确译码。
综述，可以得到解码公式为：
            (Number-range low)/range=>number
其中，number为字符串的编码。

matlab 程序:

%I=imread('001.bmp')
%imshow(I);
clear
I=[3 3 1 1 3 3 1 2;2 3 3 1 3 2 3 2;1 2 3 3 3 3 1 2];
%I=[1 1 1 1 0 0 1 0 1 1 1 0];
[m,n]=size(I);
% 第一列为灰度值,第二列为个数,第三列为概率百分数，应该也可以用imhist
table = tabulate(I();  % 注意的是，tabulate要求I的元素必须为非负整数
% 否则，以采用如下方法求解
% 如[1 2 3；1 2 2]，则统计出结果1是2个，2是3个，3是1个
% sortM=sort(M();                                                              
% uniqueM=([diff(sortM);1]>0);                                                   
% count = [sortM(uniqueM) diff(find([1;uniqueM]))]    

% 即color,p如下所示
color = table(:,1)';
p = table(:,3)'/100;

% 计算上下限
csump = cumsum(table(:,3)');
allLow =[0,csump(1:end-1)/100];
allHigh = csump/100;

numberlow = 0;
numberhigh = 1;

for k = 1:m
   for kk = 1:n
       data = I(k,kk);
       low = allLow(data==color);
       high = allHigh(data==color);
       range = numberhigh-numberlow;
       tmp = numberlow;
       numberlow = tmp+range*low;
       numberhigh = tmp+range*high;
   end
end
% the result range
fprintf('算术编码范围下限为%16.15f\n\n',numberlow);
fprintf('算术编码范围上限为%16.15f\n\n',numberhigh);


% decoding
Mat = zeros(m,n);
for k = 1:m
   for kk = 1:n
       indx = numberlow<allLow;
       indx = [indx 1];
       ind = diff(indx);
       ind = logical(ind);
       Mat(k,kk) = color(ind);
       low = allLow(ind);
       high = allHigh(ind);
       range = high - low;
       numberlow = numberlow-low;
       numberlow = numberlow/range;
   end
end

fprintf('原矩阵为:\n')
disp(I);
fprintf('\n');
fprintf('解码矩阵:\n');
disp(Mat);