arith06.m

算术编码

function varargout = Arith06(xC)
% Arith06     Arithmetic encoder or decoder 
% Vectors of integers are arithmetic encoded, 
% these vectors are collected in a cell array, xC.
% If first argument is a cell array the function do encoding,
% else decoding is done.
%
% [y, Res] = Arith06(xC);                    % encoding
% y = Arith06(xC);                           % encoding
% xC = Arith06(y);                           % decoding
% ------------------------------------------------------------------
% Arguments:
%  y        a column vector of non-negative integers (bytes) representing 
%           the code, 0 <= y(i) <= 255. 
%  Res      a matrix that sum up the results, size is (NumOfX+1)x4
%           one line for each of the input sequences, the columns are
%           Res(:,1) - number of elements in the sequence
%           Res(:,2) - unused (=0) 
%           Res(:,3) - bits needed to code the sequence
%           Res(:,4) - bit rate for the sequence, Res(:,3)/Res(:,1)
%           Then the last line is total (which include bits needed to store NumOfX)
%  xC       a cell array of column vectors of integers representing the
%           symbol sequences. (should not be to large integers)
%           If only one sequence is to be coded, we must make the cell array
%           like: xC=cell(2,1); xC{1}=x; % where x is the sequence
% ------------------------------------------------------------------
% Note: this routine is extremely slow since it is all Matlab code
% This function do recursive encoding like Huff06.
% An alternative (a perhaps better) aritmethic coder is Arith07,
% which is a more "pure" arithmetic coder

% SOME NOTES ON THE FUNCTION
% The descrition of the encoding algorithm is in 
% chapter 5 of "The Data Compression Book" by Mark Nelson. 
% The actual coding algorithm is practical identical, it is a translation
% from C code to MatLab code, but some differences have been made.
% The system model, T, keep record of the symbols that have been encoded. 
% Based on this table the probabiltity of each symbol is estimated. Probability 
% for symbol m is: (T(m+1)-T(m+2))/T(1)
% The symbols are 0,1,...,M and Escape (M+1), Escape is used to indicate an 
% unused symbol, which is then coded by another table, the Tu table.
% POSSIBLE IMPROVEMENTS
% - better decision wether to split a sequence or not 
% - for long sequences, update frequency table T=floor(T*a)  (ex: 0.2 < a < 0.9)
%   and do this for every La samples (ex: 100 < La < 5000)
%   We must not set any non-zero probabilities to zero during this adaption!!
% - Display some information (so users know something is happening)

%----------------------------------------------------------------------
% Copyright (c) 1999-2001.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  14.04.1999  KS: Function made
% Ver. 2.0  10.04.2001  KS: made function more like Huff06
%----------------------------------------------------------------------

% these global variables are used to read from or write to the compressed sequence
global y Byte BitPos      
% and these are used by the subfunctions for arithmetic coding
global high low range ub hc lc sc K code

Mfile='Arith06';
K=24;  % number of bits to use in integers  (4 <= K <= 24)
Display=1;      % display progress and/or results

% check input and output arguments, and assign values to arguments
if (nargin < 1); 
   error([Mfile,': function must have input arguments, see help.']); 
end
if (nargout < 1); 
   error([Mfile,': function must have output arguments, see help.']); 
end

if (~iscell(xC))
   Encode=0;Decode=1;
   y=xC(:);            % first argument is y
   y=[y;0;0;0;0];      % add some zeros to always have bits available
else
   Encode=1;Decode=0;
   NumOfX = length(xC);
end

Byte=0;BitPos=1;         % ready to read/write into first position
low=0;high=2^K-1;ub=0;   % initialize the coder
code=0;

if Encode
   Res=zeros(NumOfX,4);
   % initalize the global variables
   y=zeros(10,1);    % put some zeros into y initially
   % start encoding, first write VLIC to give number of sequences
   PutVLIC(NumOfX);
   % now encode each sequence continuously
   Ltot=0;
   for num=1:NumOfX
      x=xC{num};
      x=full(x(:));        % make sure x is a non-sparse column vector
      L=length(x);Ltot=Ltot+L;
      y=[y(1:Byte);zeros(50+2*L,1)];  % make more space available in y
      % now find some info about x to better code it
      maxx=max(x);
      minx=min(x);
      rangex=maxx-minx+1;
      if (minx<0) 
         Negative=1; 
         maxx=max(abs(x));
         minx=min(abs(x));
         rangex=maxx-minx+1;
      else 
         Negative=0; 
      end
      if ( (((rangex*4)>L) | (rangex>1023)) & (L>1) & (maxx>minx))  
         LogCode=1;    % this could be 0 if LogCode is not wanted
      else
         LogCode=0;
      end
      PutABit(LogCode);
      PutABit(Negative);
      I=find(x);                      % non-zero entries in x
      Sg=(sign(x(I))+1)/2;            % the signs may be needed later, 0/1
      x=abs(x);   
      if LogCode
         xa=x;                        % additional bits
         x(I)=floor(log2(x(I)));
         xa(I)=xa(I)-2.^x(I);
         x(I)=x(I)+1;
      end
      bits=EncodeVector(x);           % store the (abs and/or log) values
      if Negative                     % store the signs
         for i=1:length(Sg); PutABit(Sg(i)); end;   
         bits=bits+length(Sg);
      end
      if LogCode                     % store the additional bits
         for i=1:L
            for ii=(x(i)-1):(-1):1
               PutABit(bitget(xa(i),ii));
            end
         end
         bits=bits+sum(x)-length(I);
      end
      if L>0; Res(num,1)=L; else Res(num,1)=1; end;
      Res(num,2)=0;
      Res(num,3)=bits;
      if Display
         disp([Mfile,': Sequence ',int2str(num),' of ',int2str(L),' symbols ',...
               'encoded using ',int2str(bits),' bits.']);
      end
   end
   % flush the arithmetic coder
   PutBit(bitget(low,K-1));     
   ub=ub+1;
   while ub>0
      PutBit(~bitget(low,K-1));
      ub=ub-1;
   end
   % flush is finished
   y=y(1:Byte);   
   varargout(1) = {y};
   if (nargout >= 2) 
      % now calculate results for the total
      if Ltot<1; Ltot=1; end;   % we do not want Ltot to be zero
      Res(NumOfX+1,3)=Byte*8;
      Res(NumOfX+1,1)=Ltot;
      Res(NumOfX+1,2)=0;
      Res(:,4)=Res(:,3)./Res(:,1);
      varargout(2) = {Res}; 
   end
end

if Decode
   for k=1:K
      code=code*2;
      code=code+GetBit;   % read bits into code
   end
   NumOfX=GetVLIC;   % first read number of sequences
   xC=cell(NumOfX,1);
   for num=1:NumOfX
      LogCode=GetABit;
      Negative=GetABit;
      x=DecodeVector;   % get the (abs and/or log) values
      L=length(x);
      I=find(x);
      if Negative
         Sg=zeros(size(I));
         for i=1:length(I); Sg(i)=GetABit; end;   % and the signs   (0/1)
         Sg=Sg*2-1;                               % (-1/1)
      else
         Sg=ones(size(I));                    
      end
      if LogCode          % read additional bits too
         xa=zeros(L,1);
         for i=1:L
            for ii=2:x(i)
               xa(i)=2*xa(i)+GetABit;
            end
         end
         x(I)=2.^(x(I)-1);
         x=x+xa;
      end
      x(I)=x(I).*Sg;
      xC{num}=x;
   end
   varargout(1) = {xC}; 
end

return     % end of main function, Arith06
%----------------------------------------------------------------------
%----------------------------------------------------------------------

% ------- The main functions: EncodeVector and DecodeVector ---------------
% the EncodeVector and DecodeVector functions are the ones
% where actual coding is going on.
% These function may call themselves recursively
function bits = EncodeVector(x)
global y Byte BitPos
global high low range ub hc lc sc K code

StartPos=Byte*8-BitPos;     % used for counting bits
L=length(x);
DoSplit=0;
if L>50
   % try to find a good way to decide if sequence should be splitted
   Hi=IntHist(x,min(x),max(x));  % find the histogram
   Hinz=nonzeros(Hi);
   temp=length(Hinz);
   ent=log2(L)-sum(Hinz.*log2(Hinz))/L;  % find entropy
   % find x1 and x2  (that is do the split)
   xm=median(x);       % median in MatLab is slow
   x1=zeros(L,1);x2=zeros(L,1);
   x2(1)=x(1);i1=0;i2=1;
   for i=2:L
      if (x(i-1) <= xm) 
         i1=i1+1; x1(i1)=x(i);
      else
         i2=i2+1; x2(i2)=x(i);
      end
   end
   x1=x1(1:i1);x2=x2(1:i2);
   %
   L1=length(x1);
   Hi=IntHist(x1,min(x1),max(x1));  % find the histogram
   Hinz=nonzeros(Hi);
   ent1=log2(L1)-sum(Hinz.*log2(Hinz))/L1;  % find entropy
   L2=length(x2);
   Hi=IntHist(x2,min(x2),max(x2));  % find the histogram
   Hinz=nonzeros(Hi);
   ent2=log2(L2)-sum(Hinz.*log2(Hinz))/L2;  % find entropy
   % display results
   if 0
   disp(['Arith06-EncodeVector: sequence x is length ',int2str(L),...
         ' entropy ',num2str(ent),' giving bits ',int2str(ceil(L*ent))]);
   disp(['L1=',int2str(L1),' ent1=',num2str(ent1),' bits1=',int2str(ceil(L1*ent1)),...
      ',  L2=',int2str(L2),' ent2=',num2str(ent2),' bits2=',int2str(ceil(L2*ent2))]);
   disp(['Difference is ',int2str(ceil(L*ent-L1*ent1-L2*ent2)),...
         ' and number of symbols is ',int2str(temp)]);
   end
   % decision
   if (L*ent-L1*ent1-L2*ent2)>(temp*5)   
      DoSplit=1;
   end
end

% Handle some special possible exceptions,
if L==0
   PutABit(0);      % indicate that a sequence is coded
   PutVLIC(0);      % with length 0 (0 is 6 bits)
   PutABit(0);      % 'confirm' this by a '0'
   bits=8;
   return    % end of EncodeVector
end
if L==1
   PutABit(0);      % indicate that a sequence is coded
   PutVLIC(1);      % with length 1 (6 bits) 
   PutVLIC(x(1));   % containing this integer    
   EndPos=Byte*8-BitPos;     % used for counting bits
   bits=EndPos-StartPos;
   return    % end of EncodeVector
end
if max(x)==min(x)
   PutABit(0);      % indicate that a sequence is coded
   PutVLIC(0);      % with length 0 (0 is 6 bits)
   PutABit(1);      % 'deny' this by a '1'
   PutVLIC(L);      % actual length is L  
   PutVLIC(x(1));   % all entries containing this integer    
   EndPos=Byte*8-BitPos;     % used for counting bits
   bits=EndPos-StartPos;
   return    % end of EncodeVector
end
% end of specail cases

if DoSplit       % we split by using the median
   PutABit(1);       % indicate this sequence is splitted into two
   bits1=EncodeVector(x1);
   bits2=EncodeVector(x2);
   bits=bits1+bits2+1;
else
   PutABit(0);       % indicate that a sequence is coded
   PutVLIC(L);       % of this length
   M0=min(x);
   if (M0==0)