arith06.m

算术编码

      PutABit(0);      % indicate that M0==0
   else
      PutABit(1);      % indicate that M0~=0
      PutVLIC(M0);     % some bits for M0
      x=x-M0;          % translate x
   end
   M=max(x);
   PutVLIC(M);         % some bits for M
   % initialize model
   T=[ones(M+2,1);0];
   Tu=flipud((-1:(M+1))');   % (-1) since ESC never is used in Tu context
   % and code the symbols in the sequence x
   for l=1:L
      sc=T(1);
      m=x(l); 
      hc=T(m+1);lc=T(m+2);
      if hc==lc      % unused symbol, code ESC symbol first
         hc=T(M+2);lc=T(M+3);
         EncodeSymbol;      % code escape with T table
         sc=Tu(1);hc=Tu(m+1);lc=Tu(m+2);  % symbol with Tu table
         Tu(1:(m+1))=Tu(1:(m+1))-1;       % update Tu table
      end
      EncodeSymbol;  % code actual symbol with T table (or Tu table)
      % update T table, MUST be identical in EncodeVector and DecodeVector
      T(1:(m+1))=T(1:(m+1))+1; 
      if (rem(l,5000)==0)  
         dT=T(1:(M+2))-T(2:(M+3));
         dT=floor(dT*7/8+1/8);
         for m=(M+2):(-1):1; T(m)=T(m+1)+dT(m); end;
      end
   end
   EndPos=Byte*8-BitPos;     % used for counting bits
   bits=EndPos-StartPos;
end
return    % end of EncodeVector

function x = DecodeVector
global y Byte BitPos
global high low range ub hc lc sc K code

if GetABit
   x1=DecodeVector;
   x2=DecodeVector;
   L=length(x1)+length(x2);
   xm=median([x1;x2]);
   x=zeros(L,1);
   x(1)=x2(1);
   i1=0;i2=1;
   for i=2:L
      if (x(i-1) <= xm) 
         i1=i1+1; x(i)=x1(i1);
      else
         i2=i2+1; x(i)=x2(i2);
      end
   end
else
   L=GetVLIC;
   if (L>1)      % the normal (?) situation
      x=zeros(L,1);
      if GetABit
         M0=GetVLIC;   
      else
         M0=0;
      end
      M=GetVLIC;   
      % initialize model
      T=[ones(M+2,1);0];
      Tu=flipud((-1:(M+1))');   % (-1) since ESC never is used in Tu context
      % and decode the symbols in the sequence x
      for l=1:L
         sc=T(1);
         range=high-low+1;
         count=floor(( (code-low+1)*sc-1 )/range);
         m=2; while (T(m)>count); m=m+1; end; 
         hc=T(m-1);lc=T(m);m=m-2;
         RemoveSymbol;
         if (m>M)     % decoded ESC symbol, find symbol from Tu table
            sc=Tu(1);range=high-low+1;
            count=floor(( (code-low+1)*sc-1 )/range);
            m=2; while (Tu(m)>count); m=m+1; end; 
            hc=Tu(m-1);lc=Tu(m);m=m-2;
            RemoveSymbol;
            Tu(1:(m+1))=Tu(1:(m+1))-1;   % update Tu table
         end
         x(l)=m;
         % update T table, MUST be identical in EncodeVector and DecodeVector
         T(1:(m+1))=T(1:(m+1))+1; 
         if (rem(l,5000)==0)  
            dT=T(1:(M+2))-T(2:(M+3));
            dT=floor(dT*7/8+1/8);
            for m=(M+2):(-1):1; T(m)=T(m+1)+dT(m); end;
         end
      end
      if M0~=0
         x=x+M0;
      end
   elseif L==0
      if GetABit
         % length 0 was not confirmed
         La=GetVLIC;  % actual length
         x=GetVLIC*ones(La,1);
      else
         x=[];    % this was really a length 0 sequence
      end
   elseif L==1
      x=GetVLIC;
   else
      error('DecodeVector: illegal length of sequence.');
   end
end
return    % end of DecodeVector


% ------- Other subroutines ------------------------------------------------

% Functions to write and read a Variable Length Integer Code word
% This is a way of coding non-negative integers that uses fewer 
% bits for small integers than for large ones. The scheme is:
%   '00'   +  4 bit  - integers from 0 to 15
%   '01'   +  8 bit  - integers from 16 to 271
%   '10'   + 12 bit  - integers from 272 to 4367
%   '110'  + 16 bit  - integers from 4368 to 69903
%   '1110' + 20 bit  - integers from 69940 to 1118479
%   '1111' + 24 bit  - integers from 1118480 to 17895695
%   not supported  - integers >= 17895696 (=2^4+2^8+2^12+2^16+2^20+2^24)
function PutVLIC(N)
global y Byte BitPos
global high low range ub hc lc sc K code
if (N<0)
   error('Arith06-PutVLIC: Number is negative.'); 
elseif (N<16)
   PutABit(0);PutABit(0);
   for (i=4:-1:1); PutABit(bitget(N,i)); end;
elseif (N<272)
   PutABit(0);PutABit(1);
   N=N-16;
   for (i=8:-1:1); PutABit(bitget(N,i)); end;
elseif (N<4368)
   PutABit(1);PutABit(0);
   N=N-272;
   for (i=12:-1:1); PutABit(bitget(N,i)); end;
elseif (N<69940)
   PutABit(1);PutABit(1);PutABit(0);
   N=N-4368;
   for (i=16:-1:1); PutABit(bitget(N,i)); end;
elseif (N<1118480)
   PutABit(1);PutABit(1);PutABit(1);PutABit(0);
   N=N-69940;
   for (i=20:-1:1); PutABit(bitget(N,i)); end;
elseif (N<17895696)
   PutABit(1);PutABit(1);PutABit(1);PutABit(1);
   N=N-1118480;
   for (i=24:-1:1); PutABit(bitget(N,i)); end;
else
   error('Arith06-PutVLIC: Number is too large.'); 
end
return

function N=GetVLIC
global y Byte BitPos
global high low range ub hc lc sc K code
N=0;
if GetABit
   if GetABit
      if GetABit
         if GetABit
            for (i=1:24); N=N*2+GetABit; end;
            N=N+1118480;
         else
            for (i=1:20); N=N*2+GetABit; end;
            N=N+69940;
         end
      else
         for (i=1:16); N=N*2+GetABit; end;
         N=N+4368;
      end
   else
      for (i=1:12); N=N*2+GetABit; end;
      N=N+272;
   end
else
   if GetABit
      for (i=1:8); N=N*2+GetABit; end;
      N=N+16;
   else
      for (i=1:4); N=N*2+GetABit; end;
   end
end
return

% Aritmetic coding of a bit, probability is 0.5 for both 1 and 0
function PutABit(Bit)
global y Byte BitPos
global high low range ub hc lc sc K code
sc=2;
if Bit 
   hc=1;lc=0; 
else 
   hc=2;lc=1; 
end
EncodeSymbol;  % code the bit
return
   
function Bit=GetABit
global y Byte BitPos
global high low range ub hc lc sc K code
range=high-low+1;
sc=2;
count=floor(( (code-low+1)*sc-1 )/range);
if (1>count)
   Bit=1;hc=1;lc=0;
else
   Bit=0;hc=2;lc=1;
end
RemoveSymbol;
return;

% The EncodeSymbol function encode a symbol, (correspond to encode_symbol page 149)
function EncodeSymbol
global y Byte BitPos 
global high low range ub hc lc sc K code
range=high-low+1;
high=low+floor(((range*hc)/sc)-1);
low=low+floor((range*lc)/sc);
while 1          % for loop on page 149
   if bitget(high,K)==bitget(low,K)
      PutBit(bitget(high,K));
      while ub > 0
         PutBit(~bitget(high,K));
         ub=ub-1;
      end
   elseif (bitget(low,K-1) & (~bitget(high,K-1)))
      ub=ub+1;
      low=bitset(low,K-1,0);
      high=bitset(high,K-1,1);
   else
      break
   end
   low=bitset(low*2,K+1,0);
   high=bitset(high*2+1,K+1,0);
end
return

% The RemoveSymbol function removes (and fill in new) bits from
% file, y, to code
function RemoveSymbol 
global y Byte BitPos 
global high low range ub hc lc sc K code
range=high-low+1;
high=low+floor(((range*hc)/sc)-1);
low=low+floor((range*lc)/sc);
while 1
   if bitget(high,K)==bitget(low,K)
      % do nothing (shift bits out)
   elseif (bitget(low,K-1) & (~bitget(high,K-1)))
      code=bitset(code,K-1,~bitget(code,K-1));     % switch bit K-1
      low=bitset(low,K-1,0);
      high=bitset(high,K-1,1);
   else
      break
   end
   low=bitset(low*2,K+1,0);
   high=bitset(high*2+1,K+1,0);
   code=bitset(code*2+GetBit,K+1,0);
end
if (low > high); error('low > high'); end;
return

% Functions to write and read a Bit
function PutBit(Bit)
global y Byte BitPos
BitPos=BitPos-1;
if (~BitPos); Byte=Byte+1; BitPos=8; end; 
y(Byte) = bitset(y(Byte),BitPos,Bit);
return
   
function Bit=GetBit
global y Byte BitPos
BitPos=BitPos-1;
if (~BitPos); Byte=Byte+1; BitPos=8; end; 
Bit=bitget(y(Byte),BitPos);
return;
   
% this function is a variant of the standard hist function
function Hi=IntHist(W,i1,i2);
W=W(:);
L=length(W);
Hi=zeros(i2-i1+1,1);
if (i2-i1)>50
   for l=1:L
      i=W(l)-i1+1;
      Hi(i)=Hi(i)+1;
   end
else
   for i=i1:i2
      I=find(W==i);
      Hi(i-i1+1)=length(I);
   end
end
return;