function Code = huffman(p)
%HUFFMAN Builds a variable-lenth Huffman code for a symbol source.
% CODE = HUFFMAN(P) returns a Huffman code as a binary strings in 
% cell array CODE for input symbol probability vector P. Each word 
% CODE corresponds to a symbol whose probability is at the corresponding
% index of P
%
%   See also MAT2HUFF

% Check the input arguments for reasonableness
error(nargchk(1,1, nargin));
if (ndims(p)~=2) | (min(size(p))>1) | ~isreal(p) | ~isnumeric(p)
    error('P must be a real positive vector.');
end

% Global variable surviving all recursions of function 'makecode'
global Code
Code = cell(length(p),1);   % Initialize the global cell array

if length(p)>1                % When more than one symbol ...    
    p = p/sum(p);             % Confirm that p conatins probabilities ...
    s = reduce(p);            % Do Huffman source symbol reductions
    makecode(s,[]);           % Recursiveley generate the code
else
    Code={'1'};               % Else, trivial one symbol case!
end

%----------------------------------------------------------------------
function s = reduce(p);
% Create a Huffamn source reduction tree in a Matlab cell structure
% by performing source symbol reductions until there re only two reduced
% symbols remaining

s = cell(length(p),1);

% Generate a starting tree with symbol nodes 1,2,3, ... to reference 
% the symbol probabilities
for i=1:length(p)
    s{i} = i;
end

while numel(s)>2
    [p, i] = sort(p);       % Sort the symbol probabilities
    p(2) = p(1) + p(2);     % Merge the 2 lowest probabilities 
    p(1) = [];              % prune the lowest one
    
    s = s(i);               % Reorder tree for new probabilities 
    s{2} = {s{1}, s{2}};     % and merge & prune its nodes
    s(1) = [];              % to match the probabilities
end

%------------------------------------------------------------------------
function makecode(sc, codeword)
% Scan the codes of a Huffman source reduction tree recursively to generate
% the indicated variable length code words

% Global variable surviving all recursive calls
global Code

if isa(sc, 'cell')                      % For cell array nodes,
    makecode(sc{1}, [codeword 0]);      % add a 0 if the first element 
    makecode(sc{2}, [codeword 1]);      % or 1 if the second element
else                                    % For leaf (numeric) nodes,
    Code{sc} = char('0' + codeword);
end