%
% From http://electronicdesign.com/tabid/57/default.aspx?topic=algorithm_delivers_lossless_compression_to_adc_samples&catpath=&fltrTitle=&fltrSummary=&fltrPublication.aspx?nl=1
%
function [cr, m, e, d] = simpleSamplify(x)
%
%   function [cr, m, e, d] = simpleSamplify(x)            Al Wegener          August 2009
%
%   This function applies a simple lossless compression algorithm to a
%   sampled integer waveform and returns the mantissas and exponents and the
%   resulting compression ratio.
%
%   Input:      x       integer input array (array length should be a multiple of 4)
%
%   Output:     cr      compression ratio
%               m       output mantissas
%               e       1 exponent per 4 mantissas 
%               d       derivative (0, 1, or 2) that uses the fewest bits
%               
 
%   Ensure that the x array length is a multiple of 4, and make sure
%   we have a row vector.
 
N = 4 * floor(length(x)/4);
x = x(1:N);
x = x(:)';
 
%   Determine the number of bits in the input array, which is the
%   product of the array length and the number of bits per sample.
 
bitsPerSample = 1 + ceil(log(max(abs(x))) / log(2));
bits_in = N * bitsPerSample;
 
%   Generate the first and second derivatives of the input signal
 
x1 = [x(1) diff(x)];
x2 = [x(1) x1(2) diff(x1(2:end))];
 
%   Determine the length of each bit-packed array.
 
[e0, N0] = calcExps(x);
[e1, N1] = calcExps(x1);
[e2, N2] = calcExps(x2);
 
[y, i]   = min([N0 N1 N2]);         %  which derivative array requires the fewest bits ?
d = i - 1;                          %  "winner" derivative = i - 1;
 
%   Return the compressed array, i.e. the derivative that requires the fewest bits.
 
switch i
    
    case 1,   e = e0;  m = x;   N = N0;
    case 2,   e = e1;  m = x1;  N = N1;
    case 3,   e = e2;  m = x2;  N = N2;
        
    otherwise, fprintf(1,'ERROR\n');
        
end
 
cr = bits_in / N;
 
return
 
function [exps, Nbits] = calcExps(x)
 
x1    = reshape(x, 4, length(x)/4);
xmax  = min(max(abs(x1)), 32768);           %  limit input to 16-bit  samples
exps  = 1 + ceil(log(xmax+0.1) / log(2));   %  calculate exponents for each 4-sample mantissa group
exps  = (exps > 0.9) .* exps;               %  if xmax was zero, exps should be zero too!
 
mant_bits = 4 * sum(exps);                  %  Calculate the total number of mantissa bits
 
%   Calculate the number of bits used for the Huffman-encoded exponent differences
 
eDiffs    = [exps(1) diff(exps)];           %  generate the exponent differences
h_exps    = huffman(eDiffs);                %  determine how many bits are in the Huffman-encoded exponent array
exp_bits  = sum(h_exps);                    %  calculate the total number of Huffman-encoded exponent bits
 
%   Combine the mantissa bit count, the exponent bit count, and the 2-bit derivative flag.
 
Nbits = mant_bits + exp_bits + 2;           %  combine the mantissa and exponent bits, plus
                                            %  2 bits to hold the derivative in the packed array
return
 
function eBits = huffman(eDiff)
%
%   Huffman encoding of exponent differences.  This function assumes the following mapping of
%   exponent differences to Huffman tokens, and also shows the corresponding Huffman token length.
%
%       Exp difference              Huffman token           Token length (# bits)
%       --------------              -------------           ---------------------
%             0                            0                        1
%          -1, +1                   100,       101                  3
%          -2, +2                   1100,      1101                 4
%          -3, +3                   1110 0000, 1111 0000            8
%          -4, +4                   1110 0001, 1111 0001            8
%             :                              :                      :
%          -16, +16                 1110 1101, 1111 1101            8
%          
 
huffTable = 8 * ones(1,17);     %  most Huffman table entries take 8 bits, but
huffTable(1:3) = [1 3 4];       %    the first three Huffman table entries use 1, 3, and 4 bits
                                %    Exp diffs of 0, +/-1, and +/-2 usually account for more than
                                %      80% of all exp diffs, so their token length should be short!
 
j = 1 + abs(eDiff);             %   (exp diff range is {0, 16}, but huffTable indices are {1, 17})
eBits = huffTable(j);           %  calculate the Huffman token length for each exponent difference
 
return


