Implement Fixed-Point Log2 Using Lookup Table

This example shows how to implement fixed-point log2 using a lookup table. Lookup tables generate efficient code for embedded devices.

Setup

To assure that this example does not change your preferences or settings, this code stores the original state, and you will restore it at the end.

originalFormat = get(0, 'format'); format long g
originalWarningState = warning('off','fixed:fi:underflow');
originalFiprefState = get(fipref); reset(fipref)

Log2 Implementation

The log2 algorithm is summarized here.

  1. Declare the number of bits in a byte, B, as a constant. In this example, B=8.

  2. Use function fi_normalize_unsigned_8_bit_byte() described in example Normalize Data for Lookup Tables to normalize the input u>0 such that u = x * 2^n and 1 <= x < 2.

  3. Extract the upper B-bits of x. Let x_B denote the upper B-bits of x.

  4. Generate lookup table, LOG2LUT, such that the integer i = x_B - 2^(B-1) + 1 is used as an index to LOG2LUT so that log2(x_B) can be evaluated by looking up the index log2(x_B) = LOG2LUT(i).

  5. Use the remainder, r = x - x_B, interpreted as a fraction, to linearly interpolate between LOG2LUT(i) and the next value in the table LOG2LUT(i+1). The remainder, r, is created by extracting the lower w - B bits of x, where w denotes the word length of x. It is interpreted as a fraction by using function reinterpretcast().

  6. Finally, compute the output using the lookup table and linear interpolation:

log2( u ) = log2( x * 2^n )
          = n + log2( x )
          = n + LOG2LUT( i ) + r * ( LOG2LUT( i+1 ) - LOG2LUT( i ) )
function y = fi_log2lookup_8_bit_byte(u) %#codegen
    % Load the lookup table
    LOG2LUT = log2_lookup_table();
    % Remove fimath from the input to insulate this function from math
    % settings declared outside this function.
    u = removefimath(u);
    % Declare the output
    y = coder.nullcopy(fi(zeros(size(u)), numerictype(LOG2LUT), fimath(LOG2LUT)));
    B = 8; % Number of bits in a byte
    w = u.WordLength;
    for k = 1:numel(u)
        assert(u(k)>0,'Input must be positive.');
        % Normalize the input such that u = x * 2^n and 1 <= x < 2
        [x,n] = fi_normalize_unsigned_8_bit_byte(u(k));
        % Extract the high byte of x
        high_byte = storedInteger(bitsliceget(x, w, w - B + 1));
        % Convert the high byte into an index for LOG2LUT
        i = high_byte - 2^(B-1) + 1;
        % Interpolate between points.
        % The upper byte was used for the index into LOG2LUT
        % The remaining bits make up the fraction between points.
        T_unsigned_fraction = numerictype(0, w-B, w-B);
        r = reinterpretcast(bitsliceget(x,w-B,1), T_unsigned_fraction);
        y(k) = n + LOG2LUT(i) + ...
               r*(LOG2LUT(i+1) - LOG2LUT(i)) ;
    end
    % Remove fimath from the output to insulate the caller from math settings
    % declared inside this function.
    y = removefimath(y);
end

Log2 Lookup Table

Function log2_lookup_table loads the lookup table of log2 values. You can create the table by running:

B = 8;
log2_table = log2((2^(B-1) : 2^(B)) / 2^(B - 1))
function LOG2LUT = log2_lookup_table()
    B = 8;  % Number of bits in a byte
    % log2_table = log2((2^(B-1) : 2^(B)) / 2^(B - 1))
    log2_table = [0.000000000000000   0.011227255423254   0.022367813028454   0.033423001537450 ...
                  0.044394119358453   0.055282435501190   0.066089190457773   0.076815597050831 ...
                  0.087462841250339   0.098032082960527   0.108524456778169   0.118941072723507 ...
                  0.129283016944966   0.139551352398794   0.149747119504682   0.159871336778389 ...
                  0.169925001442312   0.179909090014934   0.189824558880017   0.199672344836364 ...
                  0.209453365628950   0.219168520462162   0.228818690495881   0.238404739325079 ...
                  0.247927513443586   0.257387842692652   0.266786540694901   0.276124405274238 ...
                  0.285402218862248   0.294620748891627   0.303780748177103   0.312882955284355 ...
                  0.321928094887362   0.330916878114617   0.339850002884625   0.348728154231078 ...
                  0.357552004618084   0.366322214245816   0.375039431346925   0.383704292474052 ...
                  0.392317422778760   0.400879436282184   0.409390936137702   0.417852514885898 ...
                  0.426264754702098   0.434628227636725   0.442943495848728   0.451211111832329 ...
                  0.459431618637297   0.467605550082997   0.475733430966398   0.483815777264256 ...
                  0.491853096329675   0.499845887083205   0.507794640198696   0.515699838284042 ...
                  0.523561956057013   0.531381460516312   0.539158811108031   0.546894459887637 ...
                  0.554588851677637   0.562242424221073   0.569855608330948   0.577428828035749 ...
                  0.584962500721156   0.592457037268080   0.599912842187128   0.607330313749611 ...
                  0.614709844115208   0.622051819456376   0.629356620079610   0.636624620543649 ...
                  0.643856189774725   0.651051691178929   0.658211482751795   0.665335917185176 ...
                  0.672425341971496   0.679480099505446   0.686500527183218   0.693486957499325 ...
                  0.700439718141092   0.707359132080883   0.714245517666123   0.721099188707185 ...
                  0.727920454563199   0.734709620225838   0.741466986401147   0.748192849589460 ...
                  0.754887502163469   0.761551232444479   0.768184324776926   0.774787059601173 ...
                  0.781359713524660   0.787902559391432   0.794415866350106   0.800899899920305 ...
                  0.807354922057604   0.813781191217037   0.820178962415188   0.826548487290915 ...
                  0.832890014164742   0.839203788096944   0.845490050944375   0.851749041416058 ...
                  0.857980995127572   0.864186144654280   0.870364719583405   0.876516946565000 ...
                  0.882643049361841   0.888743248898259   0.894817763307943   0.900866807980749 ...
                  0.906890595608518   0.912889336229962   0.918863237274595   0.924812503605781 ...
                  0.930737337562886   0.936637939002571   0.942514505339240   0.948367231584678 ...
                  0.954196310386875   0.960001932068081   0.965784284662087   0.971543553950772 ...
                  0.977279923499916   0.982993574694310   0.988684686772166   0.994353436858858 ...
                  1.000000000000000];
     % Cast to fixed point with the most accurate rounding method
     WL = 4*B;  % Word length
     FL = 2*B;  % Fraction length
     LOG2LUT = fi(log2_table,1,WL,FL,'RoundingMethod','Nearest');
     % Set fimath for the most efficient math operations
     F = fimath('OverflowAction','Wrap',...
                'RoundingMethod','Floor',...
                'SumMode','SpecifyPrecision',...
                'SumWordLength',WL,...
                'SumFractionLength',FL,...
                'ProductMode','SpecifyPrecision',...
                'ProductWordLength',WL,...
                'ProductFractionLength',2*FL);
     LOG2LUT = setfimath(LOG2LUT,F);
 end

Example

u = fi(linspace(0.001,20,100));

y = fi_log2lookup_8_bit_byte(u);

y_expected = log2(double(u));
%%3
clf
subplot(211)
plot(u,y,u,y_expected)
legend('Output','Expected output','Location','Best')

subplot(212)
plot(u,double(y)-y_expected,'r')
legend('Error')
figure(gcf)

Cleanup

Restore original state.

set(0, 'format', originalFormat);
warning(originalWarningState);
fipref(originalFiprefState);