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