arith07.m

算术编码

MetS=[0,1,2,3];       % the different methods, direct, split, diff, diff+split
MetC=[9,3,3,1];       % and the counts (which gives the probabilities)
% first set how many bits needed to code the method
b0=log2(sum(MetC))-log2(MetC(1));  
b1=log2(sum(MetC))-log2(MetC(2));  
b2=log2(sum(MetC))-log2(MetC(3));  
b3=log2(sum(MetC))-log2(MetC(4));  

I=find(x(1:(L-1))==1);      % except last element in x
J=find(x(1:(L-1))==0);
L11=length(I);   
% perhaps is 'entropy' interesting
% N1=L11+x(L);
% N0=L-N1;
% b=N1*log2(L/N1)+N0*log2(L/N0);
% disp(['L=',int2str(L),',  N1=',int2str(N1),',  (e)bits=',num2str(b)]);
b0=b0+BitEst(L,L11+x(L));   % bits needed for the sequence without splitting
if Display
   disp(['EncodeBin: x of length ',int2str(L),' and ',int2str(L11+x(L)),...
         ' ones (p=',num2str((L11+x(L))/L,'%1.3f'),...
         ') can be coded by ',num2str(b0,'%6.0f'),' bits.']);
end
% diff with a binary sequence is to indicate wether or not a symbol is
% the same as preceeding symbol or not, x(0) is assumed to be '0' (zero).
% This is the DPCM coding scheme on a binary sequence
x3=abs(x-[0;x(1:(L-1))]);
I3=find(x3(1:(L-1))==1);      % except last element in x3
J3=find(x3(1:(L-1))==0);
L31=length(I3);   
b2=b2+BitEst(L,L31+x3(L));  % bits needed for the sequence without splitting
if Display
   disp(['EncodeBin: diff x, L=',int2str(L),' gives ',int2str(L31+x3(L)),...
         ' ones (p=',num2str((L31+x3(L))/L,'%1.3f'),...
         ') can be coded by ',num2str(b2,'%6.0f'),' bits.']);
end
%
if (L>40)   
   % only now we try to split the sequences x and x3
   if (L11>3) & ((L-L11)>3) 
      % try to split x into x1 and x2, depending on previous symbol
      if L11<(L/2)
         x1=x(I+1);
         x2=[x(1);x(J+1)];
         L1=L11;L2=L-L11;
      else
         x1=[x(1);x(I+1)];
         x2=x(J+1);
         L1=L11+1;L2=L-L11-1;
      end
      b11=BitEst(L1,length(find(x1)));  % bits needed for x1
      b12=BitEst(L2,length(find(x2)));  % bits needed for x2
      % b1 is bits to code: Method, L11, x1 and x2
      b1=b1+log2(L)+b11+b12;  
      if Display
         disp(['EncodeBin, x -> x1+x2:      lengths are ',int2str(L1),'+',int2str(L2),...
               '  bits are ',num2str(b11,'%6.0f'),'+',num2str(b12,'%6.0f'),...
               '.  Total is ',num2str(b1,'%6.0f'),' bits.']);
      end
   else
      b1=b0+1;      % just to make this larger
   end
   if (L31>3) & ((L-L31)>3) 
      % try to split x3 into x4 and x5, depending on previous symbol
      if L31<(L/2)
         x4=x3(I3+1);
         x5=[x3(1);x3(J3+1)];
         L4=L31;L5=L-L31;
      else
         x4=[x3(1);x3(I3+1)];
         x5=x3(J3+1);
         L4=L31+1;L5=L-L31-1;
      end
      b31=BitEst(L4,length(find(x4)));  % bits needed for x4
      b32=BitEst(L5,length(find(x5)));  % bits needed for x5
      % b3 is bits to code: Method, L31, x4 and x5
      b3=b3+log2(L)+b31+b32;  
      if Display
         disp(['EncodeBin, diff x -> x4+x5: lengths are ',int2str(L4),'+',int2str(L5),...
               '  bits are ',num2str(b31,'%6.0f'),'+',num2str(b32,'%6.0f'),...
               '.  Total is ',num2str(b3,'%6.0f'),' bits.']);
      end
   else
      b3=b2+1;      % just to make this larger
   end
else
   b1=b0+1;      % just to make this larger
   b3=b2+1;      % just to make this larger
end
% now code x by the best method of those investigated
[b,MetI]=min([b0,b1,b2,b3]);
MetI=MetI-1;
PutS(MetI,MetS,MetC);   % code which method to use
%
if MetI==0
   % code the sequence x
   N1=L11+x(L);
   N0=L-N1;
   PutN(N1,L);            % code N1, 0<=N1<=L
   for n=1:L
      if ~(N0*N1); break; end;
      PutS(x(n),[0,1],[N0,N1]);   % code x(n)
      N0=N0-1+x(n);N1=N1-x(n);    % update model (of rest of the sequence)
   end
elseif MetI==1
   % code x1 and x2
   clear x4 x5 x3 x I3 J3 I J
   PutN(L11,L-1);         % 0<=L11<=(L-1) 
   EncodeBin(x1,L1);
   EncodeBin(x2,L2);
elseif MetI==2
   % code the sequence x3
   N1=L31+x3(L);
   N0=L-N1;
   PutN(N1,L);            % code N1, 0<=N1<=L
   for n=1:L
      if ~(N0*N1); break; end;
      PutS(x3(n),[0,1],[N0,N1]);   % code x3(n)
      N0=N0-1+x3(n);N1=N1-x3(n);    % update model (of rest of the sequence)
   end
elseif MetI==3
   % code x4 and x5
   clear x1 x2 x3 x I3 J3 I J
   PutN(L31,L-1);         % 0<=L31<=(L-1) 
   EncodeBin(x4,L4);
   EncodeBin(x5,L5);
end
return    % end of EncodeBin

function x = DecodeBin(L)
global y Byte BitPos
global high low range ub hc lc sc K code

% these must be as in EncodeBin
MetS=[0,1,2,3];       % the different methods, direct, split, diff, diff+split
MetC=[9,3,3,1];       % and the counts (which gives the probabilities)
MetI=GetS(MetS,MetC);    % encode which method to use

if (MetI==1) | (MetI==3)   % a split was done
   L11=GetN(L-1);         % 0<=L11<=(L-1) 
   if L11<(L/2)
      L1=L11;L2=L-L11;
   else
      L1=L11+1;L2=L-L11-1;
   end
   x1=DecodeBin(L1);
   x2=DecodeBin(L2);
   % build sequence x from x1 and x2
   x=zeros(L,1);
   if L11<(L/2)
      x(1)=x2(1);
      n1=0;n2=1;   % index for the last in x1 and x2
   else
      x(1)=x1(1);
      n1=1;n2=0;   % index for the last in x1 and x2
   end
   for n=2:L
      if (x(n-1))
         n1=n1+1;
         x(n)=x1(n1);
      else
         n2=n2+1;
         x(n)=x2(n2);
      end
   end
else                      % no split
   N1=GetN(L);
   N0=L-N1;
   x=zeros(L,1);
   for n=1:L
      if (N0==0); x(n:L)=1; break; end;
      if (N1==0); break; end;
      x(n)=GetS([0,1],[N0,N1]);   % decode x(n)
      N0=N0-1+x(n);N1=N1-x(n);    % update model (of rest of the sequence)
   end
end

if (MetI==2) | (MetI==3)   % x is diff coded
   for n=2:L
      x(n)=x(n-1)+x(n);
   end
   x=rem(x,2);
end

return    % end of DecodeBin


% ------- 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('Arith07-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('Arith07-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 number or symbol x, the different symbols (numbers)
% are given in the array S, and the counts (which gives the probabilities)
% are given in C, an array of same length as S.
% x must be a number in S.
% example with symbols 0 and 1, where probabilites are P{0}=0.8, P{1}=0.2
% PutS(x,[0,1],[4,1]) and x=GetS([0,1],[4,1])
% An idea: perhaps the array S should be only in the calling function, and
% that it do not need to be passed as an argument at all.
% Hint: it may be best to to the most likely symbols (highest counts) in
% the beginning of the tables S and C.
function PutS(x,S,C)
global y Byte BitPos
global high low range ub hc lc sc K code
N=length(S);     % also length(C) 
m=find(S==x);    % m is a single value, index in S, 1 <= m <= N
sc=sum(C);
lc=sc-sum(C(1:m));
hc=lc+C(m);
% disp(['PutS: lc=',int2str(lc),' hc=',int2str(hc),' sc=',int2str(sc),' m=',int2str(m)]);
EncodeSymbol;  % code the bit
return

function x=GetS(S,C)
global y Byte BitPos
global high low range ub hc lc sc K code
range=high-low+1;
sc=sum(C);
count=floor(( (code-low+1)*sc-1 )/range);
m=1;
lc=sc-C(1);
while (lc>count); m=m+1; lc=lc-C(m); end; 
hc=lc+C(m);
x=S(m);
% disp(['GetS: lc=',int2str(lc),' hc=',int2str(hc),' sc=',int2str(sc),' m=',int2str(m)]);
RemoveSymbol;
return

% Aritmetic coding of a number x, 0<=x<=N, P{0}=P{1}=...=P{N}=1/(N+1)
function PutN(x,N)     % 0<=x<=N 
global y Byte BitPos
global high low range ub hc lc sc K code
sc=N+1;
lc=x;
hc=x+1;
EncodeSymbol;  % code the bit
return

function x=GetN(N)
global y Byte BitPos
global high low range ub hc lc sc K code
range=high-low+1;
sc=N+1;
x=floor(( (code-low+1)*sc-1 )/range);
hc=x+1;lc=x;
RemoveSymbol;
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
   
function b=BitEst(N,N1);
if (N1>(N/2)); N1=N-N1; end;
N0=N-N1;
if (N>1000)
   b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-1.3256;
elseif (N1>20)
   b=(N+3/2)*log2(N)-(N0+1/2)*log2(N0)-(N1+1/2)*log2(N1)-0.020984*log2(log2(N))-1.25708;
else
   b=log2(N+1)+sum(log2(N-(0:(N1-1))))-sum(log2(N1-(0:(N1-1))));  
end
return