function varargout = Huff06(xC, ArgLevel, ArgSpeed)
% Huff06      Huffman encoder/decoder with (or without) recursive splitting
% Vectors of integers are Huffman 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] = Huff06(xC, Level, Speed);                    % encoding
% y = Huff06(xC);                                         % encoding
% xC = Huff06(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) - zero-order entropy of the sequence
%           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
%  Level    How many levels of splitting that is allowed, legal values 1-8
%           If Level=1, no further splitting of the sequences will be done
%           and there will be no recursive splitting.
%  Speed    For complete coding set Speed to 0. Set Speed to 1 to cheat 
%           during encoding, y will then be a sequence of zeros only,
%           but it will be of correct length and the other output 
%           arguments will be correct. 
% ------------------------------------------------------------------

% SOME NOTES ON THE FUNCTION
% huff06 depends on other functions for Huffman code, and the functions in this file
%  HuffLen     - find length of codewords (HL)
%  HuffTabLen  - find bits needed to store Huffman table information (HL)
%  HuffCode    - find huffman codewords
%  HuffTree    - find huffman tree

%----------------------------------------------------------------------
% Copyright (c) 1999-2000.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@tn.his.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  13.06.2000  KS: Function made based on huff04
% Ver. 1.1  20.06.2000  KS: Handle some more exceptions
% Ver. 1.2  21.06.2000  KS: Handle also negative values
% Ver. 1.3  23.06.2000  KS: Use logarithms for some sequences (line 114)
% Ver. 1.4  31.07.2000  KS: If a sequence has many zeros, Run + Value coding
%   is done. (from line 255 and some more)
% Ver. 1.5  02.08.2000  KS: May have larger integers in PutVLIC and GetVLIC
% Ver. 1.6  18.01.2001  KS: MaxL in line 218 was reduced from 2^16 to 50000.
%   For some sequences we may have length of code word larger than 16, even
%   if probability was larger than 2^(-16). Ex: Hi=[12798,14241,7126,7159,3520,...
%   3512,1857,1799,1089,1092,681,680,424,431,320,304,201,204,115,118,77,83,45,...
%   40,24,26,18,14,4,12,3,3,4,2,2,0,1]', sum(Hi)=58029
% Ver. 1.7  21.08.2001  KS: MaxL in line 218 and 420 must be the same
%   We may now have long code words (also see HuffTabLen.m)
%----------------------------------------------------------------------

global y Byte BitPos Speed Level

Mfile='Huff06';
Debug=0;    % note Debug is defined in EncodeVector and DecodeVector too

% 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
else
   Encode=1;Decode=0;
   if (nargin < 3); Speed=0; else Speed=ArgSpeed; end;
   if (nargin < 2); Level=8; else Level=ArgLevel; end;
   if ((length(Speed(:))~=1)); 
      error([Mfile,': Speed argument is not scalar, see help.']); 
   end
   if Speed; Speed=1; end;
   if ((length(Level(:))~=1)); 
      error([Mfile,': Level argument is not scalar, see help.']); 
   end
   Level=floor(Level);
   if (Level < 1); Level=1; end;
   if (Level > 8); Level=8; end;
   NumOfX = length(xC);
end
   
if Encode
   Res=zeros(NumOfX,4);
   % initalize the global variables
   y=zeros(10,1);    % put some zeros into y initially
   Byte=0;BitPos=1;  % ready to write into first position
   % start encoding, first write VLIC to give number of sequences
   PutVLIC(NumOfX);
   if Debug
      disp([Mfile,' (Encode): Level=',int2str(Level),'  Speed=',int2str(Speed),...
            '  NumOfX=',int2str(NumOfX)]);
   end
   % 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);
      if (minx<0) 
         Negative=1; 
      else 
         Negative=0; 
      end
      if ( (((maxx*4)>L) | (maxx>1023)) & (L>1) & (maxx>minx))  
         % the test for LogCode could be better, I think, (ver. 1.3)
         LogCode=1;    % this could be 0 if LogCode is not wanted
      else
         LogCode=0;
      end
      PutBit(LogCode);
      PutBit(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, ent]=EncodeVector(x);   % store the (abs and/or log) values
      if Negative                    % store the signs
         for i=1:length(Sg); PutBit(Sg(i)); end;   
         bits=bits+length(Sg);
         ent=ent+length(Sg)/L;
      end
      if LogCode                     % store the additional bits
         for i=1:L
            for ii=(x(i)-1):(-1):1
               PutBit(bitget(xa(i),ii));
            end
         end
         bits=bits+sum(x)-length(I);
         ent=ent+(sum(x)-length(I))/L;
      end
      if L>0; Res(num,1)=L; else Res(num,1)=1; end;
      Res(num,2)=ent;
      Res(num,3)=bits;
   end
   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)=sum(Res(1:NumOfX,1).*Res(1:NumOfX,2))/Ltot;
      Res(:,4)=Res(:,3)./Res(:,1);
      varargout(2) = {Res}; 
   end
end

if Decode
   % initalize the global variables, y is set earlier
   Byte=0;BitPos=1;  % ready to read from first position
   NumOfX=GetVLIC;   % first read number of sequences
   if Debug
      disp([Mfile,'(Decode):  NumOfX=',int2str(NumOfX),'  length(y)=',int2str(length(y))]);
   end
   xC=cell(NumOfX,1);
   for num=1:NumOfX
      LogCode=GetBit;
      Negative=GetBit;
      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)=GetBit; 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)+GetBit;
            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, huff06

% the EncodeVector and DecodeVector functions are the ones
% where actual coding is going on.
% This function calls itself recursively
function [bits, ent] = EncodeVector(x, bits, HL, Maxx, Meanx)
global y Byte BitPos Speed Level 
Debug=0;
Level = Level - 1;
MaxL=50000;         % longer sequences is split in the middle
L=length(x);
% first handle some special possible exceptions,
if L==0
   PutBit(0);       % indicate that a sequence is coded
   PutVLIC(L);      % with length 0 (0 is 6 bits)
   PutBit(0);       % 'confirm' this by a '0', Run + Value is indicated by a '1'
   bits=2+6;
   ent=0;
   Level = Level + 1;
   return    % end of EncodeVector
end
if L==1
   PutBit(0);       % indicate that a sequence is coded
   PutVLIC(L);      % with length 1 (6 bits) 
   PutVLIC(x(1));   % containing this integer    
   bits=1+2*6;
   if (x(1)>=16); bits=bits+4; end;
   if (x(1)>=272); bits=bits+4; end;
   if (x(1)>=4368); bits=bits+5; end;
   if (x(1)>=69904); bits=bits+5; end;
   if (x(1)>=1118480); bits=bits+4; end;
   ent=0;
   Level = Level + 1;
   return    % end of EncodeVector
end
if max(x)==min(x)
   PutBit(0);       % indicate that a sequence is coded
   PutVLIC(L);      % with length L
   for i=1:7; PutBit(1); end;   % write end of Huffman Table
   PutVLIC(x(1));   % containing this integer    
   bits=1+6+7+6;      
   if (x(1)>=16); bits=bits+4; end;
   if (x(1)>=272); bits=bits+4; end;
   if (x(1)>=4368); bits=bits+5; end;
   if (x(1)>=69904); bits=bits+5; end;
   if (x(1)>=1118480); bits=bits+4; end;
   if (L>=16); bits=bits+4; end;
   if (L>=272); bits=bits+4; end;
   if (L>=4368); bits=bits+5; end;
   if (L>=69904); bits=bits+5; end;
   if (L>=1118480); bits=bits+4; end;
   ent=0;
   Level = Level + 1;
   return    % end of EncodeVector
end
% here we test if Run + Value coding should be done
I=find(x);   % the non-zero indices of x
if (L/2-length(I))>50
   Maxx=max(x);
   Hi=IntHist(x,0,Maxx);  % find the histogram
   Hinz=nonzeros(Hi);
   ent=log2(L)-sum(Hinz.*log2(Hinz))/L;  % find entropy
   % there are few non-zero indices => Run+Value coding of x
   x2=x(I);  % the values  
   I=[I(:);L+1];   % include length of x
   for i=length(I):(-1):2; I(i)=I(i)-I(i-1); end;
   x1=I-1;   % the runs  
   % code this as an unconditional split (like if L is large)
   if Speed
      Byte=Byte+1;    % since we add 8 bits
   else
      PutBit(0);       % this is idicated like when a sequence 
      PutVLIC(0);      % of length 0 is coded, but we add one extra bit
      PutBit(1);       % Run + Value is indicated by a '1'
   end;
   [bits1, temp] = EncodeVector(x1);
   [bits2, temp] = EncodeVector(x2);
   bits=bits1+bits2+8;
   Level = Level + 1;
   return    % end of EncodeVector
end

if (nargin==1)
   Maxx=max(x);
   Meanx=mean(x);
   Hi=IntHist(x,0,Maxx);  % find the histogram
   Hinz=nonzeros(Hi);
   ent=log2(L)-sum(Hinz.*log2(Hinz))/L;  % find entropy
   HL=HuffLen(Hi);
   HLlen=HuffTabLen(HL);
   % find number of bits to use, store L, HL and x
   bits=6+HLlen+sum(HL.*Hi);
   if (L>=16); bits=bits+4; end;
   if (L>=272); bits=bits+4; end;
   if (L>=4368); bits=bits+5; end;
   if (L>=69904); bits=bits+5; end;
   if (L>=1118480); bits=bits+4; end;
   if Debug
      disp(['bits=',int2str(bits),'  HLlen=',int2str(HLlen),...
         '   HClen=',int2str(sum(HL.*Hi))]);
   end
else                % arguments are given, do not need to be calculated
   ent=0;
end
%
% Here we have: x, bits, L, HL, Maxx, Meanx, ent
if (L>MaxL)   % we split sequence anyway (and the easy way; in the middle)
   L1=ceil(L/2);L2=L-L1;
   x1=x(1:L1);x2=x((L1+1):L);
elseif ((Level > 0) & (L>10))      
   xm=median(x); % median in MatLab is slow, could be calulated faster by using the histogram
   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);
   % find bits1 and bits2 for x1 and x2
   L1=length(x1);L2=length(x2);
   Maxx1=max(x1);Maxx2=max(x2);
   Meanx1=mean(x1);Meanx2=mean(x2);
   Hi1=IntHist(x1,0,Maxx1);  % find the histogram
   Hi2=IntHist(x2,0,Maxx2);  % find the histogram
   HL1=hufflen(Hi1);HL2=hufflen(Hi2);
   HLlen1=HuffTabLen(HL1);
   HLlen2=HuffTabLen(HL2);
   bits1=6+HLlen1+sum(HL1.*Hi1);
   bits2=6+HLlen2+sum(HL2.*Hi2);
   if (L1>=16); bits1=bits1+4; end;
   if (L1>=272); bits1=bits1+4; end;
   if (L1>=4368); bits1=bits1+5; end;
   if (L1>=69904); bits1=bits1+5; end;
   if (L1>=1118480); bits1=bits1+4; end;
   if (L2>=16); bits2=bits2+4; end;
   if (L2>=272); bits2=bits2+4; end;
   if (L2>=4368); bits2=bits2+5; end;
   if (L2>=69904); bits2=bits2+5; end;
   if (L2>=1118480); bits2=bits2+4; end;
else
   bits1=bits;bits2=bits;
end
% Here we may have: x1, bits1, L1, HL1, Maxx1, Meanx1
% and               x2, bits2, L2, HL2, Maxx2, Meanx2
% but at least we have bits1 and bits2  (and bits)
if Debug
   disp(['Level=',int2str(Level),'  bits=',int2str(bits),'  bits1=',int2str(bits1),...
         '  bits2=',int2str(bits2),'  sum=',int2str(bits1+bits2)]);
end

if (L>MaxL)
   if Speed
      BitPos=BitPos-1;
      if (~BitPos); Byte=Byte+1; BitPos=8; end; 
   else
      PutBit(1);       % indicate sequence is splitted into two
   end;
   [bits1, temp] = EncodeVector(x1);
   [bits2, temp] = EncodeVector(x2);
   bits=bits1+bits2+1;
elseif ((bits1+bits2) < bits) 
   if Speed
      BitPos=BitPos-1;
      if (~BitPos); Byte=Byte+1; BitPos=8; end; 
   else
      PutBit(1);       % indicate sequence is splitted into two
   end;
   [bits1, temp] = EncodeVector(x1, bits1, HL1, Maxx1, Meanx1);
   [bits2, temp] = EncodeVector(x2, bits2, HL2, Maxx2, Meanx2);
   bits=bits1+bits2+1;
else
   bits=bits+1;      % this is how many bits we are going to write
   if Debug
      disp(['EncodeVector: Level=',int2str(Level),'  ',int2str(L),...
            ' sybols stored in ',int2str(bits),' bits.']);
   end
   if Speed
      % advance Byte and BitPos without writing to y
      Byte=Byte+floor(bits/8);
      BitPos=BitPos-mod(bits,8);
      if (BitPos<=0); BitPos=BitPos+8; Byte=Byte+1; end;
   else
      % put the bits into y
      StartPos=Byte*8-BitPos;     % control variable
      PutBit(0);       % indicate that a sequence is coded
      PutVLIC(L);       
      PutHuffTab(HL);
      HK=huffcode(HL);
      for i=1:L;
         n=x(i)+1;    % symbol number (value 0 is first symbol, symbol 1)
         for k=1:HL(n)
            PutBit(HK(n,k));
         end
      end
      % check if one has used as many bits as calculated
      BitsUsed=Byte*8-BitPos-StartPos;
      if (BitsUsed~=bits)
         disp(['L=',int2str(L),'  max(x)=',int2str(max(x)),'  min(x)=',int2str(min(x))]);
         disp(['BitsUsed=',int2str(BitsUsed),'  bits=',int2str(bits)]);
         error(['Huff06-EncodeVector: Logical error, (BitsUsed~=bits).']); 
      end
   end
end
Level = Level + 1;
return    % end of EncodeVector

function x = DecodeVector
global y Byte BitPos
MaxL=50000;      % as in the EncodeVector function (line 216)
if GetBit
   x1=DecodeVector;
   x2=DecodeVector;
   L=length(x1)+length(x2);
   if (L>MaxL)
      x=[x1(:);x2(:)];
   else
      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
   end
else
   L=GetVLIC;
   if (L>1)
      x=zeros(L,1);
      HL=GetHuffTab;
      if length(HL)
         Htree=HuffTree(HL);
         root=1;pos=root;     
         l=0;  % number of symbols decoded so far
         while l<L
            if GetBit
               pos=Htree(pos,3);
            else
               pos=Htree(pos,2);
            end
            if Htree(pos,1)           % we have arrived at a leaf
               l=l+1;
               x(l)=Htree(pos,2)-1;   % value is one less than symbol number
               pos=root;              % start at root again
            end
         end
      else     % HL has length 0, that is empty Huffman table
         x=x+GetVLIC;
      end
   elseif L==0
      if GetBit
         % this is a Run + Value coded sequence
         x1=DecodeVector;
         x2=DecodeVector;
         % now build the actual sequence
         I=x1;      % runs  
         I=I+1;
         L=length(I);  % one more than the number of values in x
         for i=2:L;I(i)=I(i-1)+I(i); end;
         x=zeros(I(L)-1,1);
         x(I(1:(L-1)))=x2;  % values
      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

% Functions to write and read the Huffman Table Information
% The format is defined in HuffTabLen, we repeat it here
% Function assume that the table information is stored in the following format
%        previous symbol is set to the initial value 2, Prev=2
%        Then we have for each symbol a code word to tell its length 
%         '0'             - same length as previous symbol
%         '10'            - increase length by 1, and 17->1
%         '1100'          - decrease length by 1, and 0->16
%         '11010'         - increase length by 2, and 17->1, 18->2
%         '11011'         - One zero, unused symbol (twice for two zeros)
%         '111xxxx'       - set code length to CL=Prev+x (where 3 <= x <= 14)
%                           and if CL>16; CL=CL-16
%        we have 4 unused 7 bit code words, which we give the meaning
%         '1110000'+4bits - 3-18 zeros
%         '1110001'+8bits - 19-274 zeros, zeros do not change previous value
%         '1110010'+4bits - for CL=17,18,...,32, do not change previous value
%         '1111111'       - End Of Table

function PutHuffTab(HL)
global y Byte BitPos

HL=HL(:);
% if (max(HL) > 32) 
%    disp(['PutHuffTab: To large value in HL, max(HL)=',int2str(max(HL))]); 
% end
% if (min(HL) < 0)
%    disp(['PutHuffTab: To small value in HL, min(HL)=',int2str(min(HL))]); 
% end
Prev=2;
ZeroCount=0;
L=length(HL);

for l=1:L
   if HL(l)==0
      ZeroCount=ZeroCount+1;
   else
      while (ZeroCount > 0)
         if ZeroCount<3
            for i=1:ZeroCount
               PutBit(1);PutBit(1);PutBit(0);PutBit(1);PutBit(1);
            end
            ZeroCount=0; 
         elseif ZeroCount<19 
            PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(0);PutBit(0);
            for (i=4:-1:1); PutBit(bitget(ZeroCount-3,i)); end;
            ZeroCount=0; 
         elseif ZeroCount<275
            PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(0);PutBit(1);
            for (i=8:-1:1); PutBit(bitget(ZeroCount-19,i)); end;
            ZeroCount=0; 
         else 
            PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(0);PutBit(1);
            for (i=8:-1:1); PutBit(1); end;
            ZeroCount=ZeroCount-274; 
         end
      end
      if HL(l)>16
         PutBit(1);PutBit(1);PutBit(1);PutBit(0);PutBit(0);PutBit(1);PutBit(0);
         for (i=4:-1:1); PutBit(bitget(HL(l)-17,i)); end;
      else
         Inc=HL(l)-Prev;
         if Inc<0; Inc=Inc+16; end;
         if (Inc==0) 
            PutBit(0);
         elseif (Inc==1) 
            PutBit(1);PutBit(0);
         elseif (Inc==2) 
            PutBit(1);PutBit(1);PutBit(0);PutBit(1);PutBit(0);
         elseif (Inc==15) 
            PutBit(1);PutBit(1);PutBit(0);PutBit(0);
         else 
            PutBit(1);PutBit(1);PutBit(1);
            for (i=4:-1:1); PutBit(bitget(Inc,i)); end;
         end
         Prev=HL(l);
      end
   end
end
for (i=7:-1:1); PutBit(1); end;       % the EOT codeword

return;  % end of PutHuffTab

function HL=GetHuffTab
global y Byte BitPos

Debug=0;
Prev=2;
ZeroCount=0;
HL=zeros(10000,1);
HLi=0;
EndOfTable=0;

while ~EndOfTable
   if GetBit
      if GetBit
         if GetBit
            Inc=0;
            for (i=1:4); Inc=Inc*2+GetBit; end;
            if Inc==0
               ZeroCount=0;
               for (i=1:4); ZeroCount=ZeroCount*2+GetBit; end;
               HLi=HLi+ZeroCount+3;
            elseif Inc==1
               ZeroCount=0;
               for (i=1:8); ZeroCount=ZeroCount*2+GetBit; end;
               HLi=HLi+ZeroCount+19;
            elseif Inc==2           % HL(l) is large, >16
               HLi=HLi+1;
               HL(HLi)=0;
               for (i=1:4); HL(HLi)=HL(HLi)*2+GetBit; end;
               HL(HLi)=HL(HLi)+17;
            elseif Inc==15
               EndOfTable=1;
            else
               Prev=Prev+Inc;
               if Prev>16; Prev=Prev-16; end;
               HLi=HLi+1;HL(HLi)=Prev;
            end
         else
            if GetBit
               if GetBit
                  HLi=HLi+1;
               else
                  Prev=Prev+2;
                  if Prev>16; Prev=Prev-16; end;
                  HLi=HLi+1;HL(HLi)=Prev;
               end
            else
               Prev=Prev-1;
               if Prev<1; Prev=16; end;
               HLi=HLi+1;HL(HLi)=Prev;
            end
         end
      else
         Prev=Prev+1;
         if Prev>16; Prev=1; end;
         HLi=HLi+1;HL(HLi)=Prev;
      end
   else
      HLi=HLi+1;HL(HLi)=Prev;
   end
end
if HLi>0
   HL=HL(1:HLi);
else
   HL=[];
end

if Debug
   % check if this is a valid Huffman table
   temp=sum(2.^(-nonzeros(HL)));
   if temp ~=1
      error(['GetHuffTab: HL table is no good, temp=',num2str(temp)]);
   end
end

return;  % end of GetHuffTab

% 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
if (N<0)
   error('Huff06-PutVLIC: Number is negative.'); 
elseif (N<16)
   PutBit(0);PutBit(0);
   for (i=4:-1:1); PutBit(bitget(N,i)); end;
elseif (N<272)
   PutBit(0);PutBit(1);
   N=N-16;
   for (i=8:-1:1); PutBit(bitget(N,i)); end;
elseif (N<4368)
   PutBit(1);PutBit(0);
   N=N-272;
   for (i=12:-1:1); PutBit(bitget(N,i)); end;
elseif (N<69940)
   PutBit(1);PutBit(1);PutBit(0);
   N=N-4368;
   for (i=16:-1:1); PutBit(bitget(N,i)); end;
elseif (N<1118480)
   PutBit(1);PutBit(1);PutBit(1);PutBit(0);
   N=N-69940;
   for (i=20:-1:1); PutBit(bitget(N,i)); end;
elseif (N<17895696)
   PutBit(1);PutBit(1);PutBit(1);PutBit(1);
   N=N-1118480;
   for (i=24:-1:1); PutBit(bitget(N,i)); end;
else
   error('Huff06-PutVLIC: Number is too large.'); 
end
return

function N=GetVLIC
global y Byte BitPos
N=0;
if GetBit
   if GetBit
      if GetBit
         if GetBit
            for (i=1:24); N=N*2+GetBit; end;
            N=N+1118480;
         else
            for (i=1:20); N=N*2+GetBit; end;
            N=N+69940;
         end
      else
         for (i=1:16); N=N*2+GetBit; end;
         N=N+4368;
      end
   else
      for (i=1:12); N=N*2+GetBit; end;
      N=N+272;
   end
else
   if GetBit
      for (i=1:8); N=N*2+GetBit; end;
      N=N+16;
   else
      for (i=1:4); N=N*2+GetBit; end;
   end
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(:);
%if (rem(i1,1) | rem(i2,1));   error('Non integers'); end;
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;


