function y = perform_arith_fixed(x, h, n)% perform_arith_fixed - arithmetic coding% % Coding:%   y = perform_arith_fixed(x, h, n);% Decoding:%   x = perform_arith_fixed(y, h);%%   Copyright (c) 2008 Gabriel Peyreif nargin==2    dir=1;elseif nargin==3    dir=-1;else    error('Wrong number of arumgnets');endh = round(h*100000); if dir==1    y = arithenco(x,h);else    y = arithdeco(x, h, n);endfunction code = arithenco(seq, counts)%ARITHENCO Encode a sequence of symbols using arithmetic coding.%   CODE = ARITHENCO(SEQ, COUNTS) generates binary arithmetic code %   corresponding to the sequence of symbols specified in the vector SEQ. %   The vector COUNTS contains the symbol counts (the number of times each%   symbol of the source's alphabet occurs in a test data set) and represents%   the source's statistics.%%   Example: %     Consider a source whose alphabet is {x, y, z}. A 177-symbol test data %     set from the source contains 29 x's, 48 y's and 100 z's. To encode the %     sequence yzxzz, use these commands:%%       seq = [2 3 1 3 3];%       counts = [29 48 100];%       code = arithenco(seq, counts) %                   %   See also ARITHDECO.%   Copyright 1996-2002 The MathWorks, Inc.%   $Revision: 1.3 $ $Date: 2002/06/17 12:22:10 $%   References:%         [1] Sayood, K., Introduction to Data Compression, %         Morgan Kaufmann, 2000, Chapter 4, Section 4.4.3.         % Check the incoming orientation and adjust if necessary[row_s, col_s] = size(seq);if (row_s > 1),    seq = seq.';end[row_c, col_c] = size(counts);if (row_c > 1),    counts = counts.';end% Compute the cumulative counts vector from the counts cum_counts = [0, cumsum(counts)];% Compute the Word Length required.total_count = cum_counts(end);N = ceil(log2(total_count)) + 2;% Initialize the lower and upper bounds.dec_low = 0;dec_up = 2^N-1;E3_count = 0;% Obtain an over estimate for the length of CODE and initialize CODEcode_len = length(seq) * ( ceil(log2(length(counts))) + 2 ) + N;code = zeros(1, code_len);code_index = 1;% Loop for each symbol in SEQfor k = 1:length(seq)    symbol = seq(k);    % Compute the new lower bound    dec_low_new = dec_low + floor( (dec_up-dec_low+1)*cum_counts(symbol+1-1)/total_count );    % Compute the new upper bound    dec_up = dec_low + floor( (dec_up-dec_low+1)*cum_counts(symbol+1)/total_count )-1;    % Update the lower bound    dec_low = dec_low_new;        % Check for E1, E2 or E3 conditions and keep looping as long as they occur.    while( isequal(bitget(dec_low, N), bitget(dec_up, N)) || ...        (isequal(bitget(dec_low, N-1), 1) && isequal(bitget(dec_up, N-1), 0) ) ),                % If it is an E1 or E2 condition,        if isequal(bitget(dec_low, N), bitget(dec_up, N)),            % Get the MSB            b = bitget(dec_low, N);            code(code_index) = b;            code_index = code_index + 1;                    % Left shifts            dec_low = bitshift(dec_low, 1) + 0;            dec_up = bitshift(dec_up, 1) + 1;                        % Check if E3_count is non-zero and transmit appropriate bits            if (E3_count > 0),                % Have to transmit complement of b, E3_count times.                code(code_index:code_index+E3_count-1) = bitcmp(b, 1).*ones(1, E3_count);                code_index = code_index + E3_count;                E3_count = 0;            end            % Reduce to N for next loop            dec_low = bitset(dec_low, N+1, 0);            dec_up  = bitset(dec_up, N+1, 0);                    % Else if it is an E3 condition            elseif ( (isequal(bitget(dec_low, N-1), 1) && ...            isequal(bitget(dec_up, N-1), 0) ) ),                        % Left shifts            dec_low = bitshift(dec_low, 1) + 0;            dec_up  = bitshift(dec_up, 1) + 1;            % Reduce to N for next loop            dec_low = bitset(dec_low, N+1, 0);            dec_up  = bitset(dec_up, N+1, 0);                        % Complement the new MSB of dec_low and dec_up            dec_low = bitxor(dec_low, 2^(N-1) );            dec_up  = bitxor(dec_up, 2^(N-1) );                        % Increment E3_count to keep track of number of times E3 condition is hit.            E3_count = E3_count+1;        end    endend % Terminate encodingbin_low = de2bi(dec_low, N, 'left-msb');if E3_count==0,    % Just transmit the final value of the lower bound bin_low           code(code_index:code_index + N - 1) = bin_low;    code_index = code_index + N;else   % Transmit the MSB of bin_low.    b = bin_low(1);   code(code_index) = b;   code_index = code_index + 1;      % Then transmit complement of b (MSB of bin_low), E3_count times.    code(code_index:code_index+E3_count-1) = bitcmp(b, 1).*ones(1, E3_count);   code_index = code_index + E3_count;   % Then transmit the remaining bits of bin_low   code(code_index:code_index+N-2) = bin_low(2:N);   code_index = code_index + N - 1;end          % Output only the filled valuescode = code(1:code_index-1);% Set the same output orientation as seqif (row_s > 1)    code = code.';end%----------------------------------------------------------function eStr = errorchk(seq, counts)% Function for validating the input parameters. eStr.ecode = 0;eStr.emsg = '';% Check to make sure a vector has been entered as input and not a matrixif (length(find(size(seq)==1)) ~= 1)    eStr.emsg = ['The symbol sequence parameter must be a vector of positive ',...                 'finite integers.'];    eStr.ecode = 1; return;    end% Check to make sure a character array is not specified for SEQif ischar(seq)    eStr.emsg = ['The symbol sequence parameter must be a vector of positive ',...                 'finite integers.'];    eStr.ecode = 1; return;    end% Check to make sure that finite positive integer values (non-complex) are % entered for SEQif ~all(seq > 0) || ~all(isfinite(seq)) || ~isequal(seq, round(seq)) || ...    ~isreal(seq)    eStr.emsg = ['The symbol sequence parameter must be a vector of positive ',...                 'finite integers.'];    eStr.ecode = 1; return;    endif length(find(size(counts)==1)) ~= 1    eStr.emsg = ['The symbol counts parameter must be a vector of positive ',...                 'finite integers.'];    eStr.ecode = 1; return;    end% Check to make sure that finite positive integer values (non-complex) are% entered for COUNTSif ~all(counts > 0) || ~all(isfinite(counts)) || ~isequal(counts, round(counts)) || ...     ~isreal(counts)    eStr.emsg = ['The symbol counts parameter must be a vector of positive ',...                 'finite integers.'];    eStr.ecode = 1; return;    end% Check to ensure that the maximum value in the SEQ vector is less than or equal% to the length of the counts vector COUNTS.if max(seq) > length(counts)    eStr.emsg = ['The symbol sequence parameter can take values only between',...                 ' 1 and the length of the symbol counts parameter.'];    eStr.ecode = 1; return;    end% [EOF]function dseq = arithdeco(code, counts, len)%ARITHDECO Decode binary code using arithmetic decoding.%   DSEQ = ARITHDECO(CODE, COUNTS, LEN) decodes the binary arithmetic code%   in the vector CODE (generated using ARITHENCO) to the corresponding%   sequence of symbols. The vector COUNTS contains the symbol counts (the%   number of times each symbol of the source's alphabet occurs in a test%   data set) and represents the source's statistics. LEN is the number of%   symbols to be decoded. %   %   Example: %     Consider a source whose alphabet is {x, y, z}. A 177-symbol test data %     set from the source contains 29 x's, 48 y's and 100 z's. To encode the %     sequence yzxzz, use these commands:%%       seq = [2 3 1 3 3];%       counts = [29 48 100];%       code = arithenco(seq, counts)   %            %     To decode this code (and recover the sequence of  %     symbols it represents) use this command:%            %       dseq = arithdeco(code, counts, 5)%            %   See also ARITHENCO.%   Copyright 1996-2002 The MathWorks, Inc.%   $Revision: 1.2 $ $Date: 2002/04/14 20:12:32 $%   References:%         [1] Sayood, K., Introduction to Data Compression, %         Morgan Kaufmann, 2000, Chapter 4, Section 4.4.3.% Check the incoming orientation and adjust if necessary[row_cd, col_cd] = size(code);if (row_cd > 1),    code = code.';end[row_c, col_c] = size(counts);if (row_c > 1),    counts = counts.';end% Compute the cumulative counts vector from the counts vectorcum_counts = [0, cumsum(counts)];% Compute the Word Length (N) required.total_count = cum_counts(end);N = ceil(log2(total_count)) + 2;% Initialize the lower and upper bounds.dec_low = 0;dec_up = 2^N-1;% Read the first N number of bits into a temporary tag bin_tagbin_tag = code(1:N);dec_tag = bi2de(bin_tag, 'left-msb');% Initialize DSEQdseq = zeros(1,len);dseq_index = 1;k=N;ptr = 0;% This loop runs untill all the symbols are decoded into DSEQwhile (dseq_index <= len)        % Compute dec_tag_new    dec_tag_new =floor( ((dec_tag-dec_low+1)*total_count-1)/(dec_up-dec_low+1) );        % Decode a symbol based on dec_tag_new    ptr = pick(cum_counts, dec_tag_new);        % Update DSEQ by adding the decoded symbol    dseq(dseq_index) = ptr;    dseq_index = dseq_index + 1;        % Compute the new lower bound    dec_low_new = dec_low + floor( (dec_up-dec_low+1)*cum_counts(ptr-1+1)/total_count );        % Compute the new upper bound    dec_up = dec_low + floor( (dec_up-dec_low+1)*cum_counts(ptr+1)/total_count )-1;        % Update the lower bound    dec_low = dec_low_new;        % Check for E1, E2 or E3 conditions and keep looping as long as they occur.     while ( isequal(bitget(dec_low, N), bitget(dec_up, N)) | ...        ( isequal(bitget(dec_low, N-1), 1) & isequal(bitget(dec_up, N-1), 0) ) ),                % Break out if we have finished working with all the bits in CODE        if ( k==length(code) ), break, end;        k = k + 1;        % If it is an E1 or E2 condition, do        if isequal(bitget(dec_low, N), bitget(dec_up, N)),            % Left shifts and update            dec_low = bitshift(dec_low, 1) + 0;            dec_up  = bitshift(dec_up,  1) + 1;            % Left shift and read in code            dec_tag = bitshift(dec_tag, 1) + code(k);            % Reduce to N for next loop            dec_low = bitset(dec_low, N+1, 0);            dec_up  = bitset(dec_up,  N+1, 0);            dec_tag = bitset(dec_tag, N+1, 0);                % Else if it is an E3 condition                elseif ( isequal(bitget(dec_low, N-1), 1) & ...            isequal(bitget(dec_up, N-1), 0) ),            % Left shifts and update            dec_low = bitshift(dec_low, 1) + 0;            dec_up  = bitshift(dec_up,  1) + 1;            % Left shift and read in code            dec_tag = bitshift(dec_tag, 1) + code(k);                        % Reduce to N for next loop            dec_low = bitset(dec_low, N+1, 0);            dec_up  = bitset(dec_up,  N+1, 0);            dec_tag = bitset(dec_tag, N+1, 0);                        % Complement the new MSB of dec_low, dec_up and dec_tag            dec_low = bitxor(dec_low, 2^(N-1) );            dec_up  = bitxor(dec_up,  2^(N-1) );            dec_tag = bitxor(dec_tag, 2^(N-1) );        end    end % end whileend % end while length(dseq)% Set the same output orientation as codeif (row_cd > 1)    dseq = dseq.';end%-------------------------------------------------------------function [ptr] = pick(cum_counts, value);% This internal function is used to find where value is positioned% Check for this case and quickly exitif value == cum_counts(end)    ptr = length(cum_counts)-1;    returnendc = find(cum_counts <= value);ptr = c(end);% EOF