libavcodec/celp_math.c

   1 /*
   2  * Various fixed-point math operations
   3  *
   4  * Copyright (c) 2008 Vladimir Voroshilov
   5  *
   6  * This file is part of Libav.
   7  *
   8  * Libav is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU Lesser General Public
  10  * License as published by the Free Software Foundation; either
  11  * version 2.1 of the License, or (at your option) any later version.
  12  *
  13  * Libav is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * Lesser General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU Lesser General Public
  19  * License along with Libav; if not, write to the Free Software
  20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  21  */
  22
  23 #include <inttypes.h>
  24 #include <limits.h>
  25 #include <assert.h>
  26
  27 #include "avcodec.h"
  28 #include "celp_math.h"
  29 #include "mathops.h"
  30
  31 #include "libavutil/common.h"
  32
  33 static const uint16_t exp2a[]=
  34 {
  35      0,  1435,  2901,  4400,  5931,  7496,  9096, 10730,
  36  12400, 14106, 15850, 17632, 19454, 21315, 23216, 25160,
  37  27146, 29175, 31249, 33368, 35534, 37747, 40009, 42320,
  38  44682, 47095, 49562, 52082, 54657, 57289, 59979, 62727,
  39 };
  40
  41 static const uint16_t exp2b[]=
  42 {
  43      3,   712,  1424,  2134,  2845,  3557,  4270,  4982,
  44   5696,  6409,  7124,  7839,  8554,  9270,  9986, 10704,
  45  11421, 12138, 12857, 13576, 14295, 15014, 15734, 16455,
  46  17176, 17898, 18620, 19343, 20066, 20790, 21514, 22238,
  47 };
  48
  49 int ff_exp2(uint16_t power)
  50 {
  51     unsigned int result= exp2a[power>>10] + 0x10000;
  52
  53     assert(power <= 0x7fff);
  54
  55     result= (result<<3) + ((result*exp2b[(power>>5)&31])>>17);
  56     return result + ((result*(power&31)*89)>>22);
  57 }
  58
  59 /**
  60  * Table used to compute log2(x)
  61  *
  62  * tab_log2[i] = (1<<15) * log2(1 + i/32), i=0..32
  63  */
  64 static const uint16_t tab_log2[33] =
  65 {
  66       4,   1459,   2870,   4240,   5572,   6867,   8127,   9355,
  67   10552,  11719,  12858,  13971,  15057,  16120,  17158,  18175,
  68   19170,  20145,  21100,  22036,  22954,  23854,  24738,  25605,
  69   26457,  27294,  28116,  28924,  29719,  30500,  31269,  32025,  32769,
  70 };
  71
  72 int ff_log2_q15(uint32_t value)
  73 {
  74     uint8_t  power_int;
  75     uint8_t  frac_x0;
  76     uint16_t frac_dx;
  77
  78     // Stripping zeros from beginning
  79     power_int = av_log2(value);
  80     value <<= (31 - power_int);
  81
  82     // b31 is always non-zero now
  83     frac_x0 = (value & 0x7c000000) >> 26; // b26-b31 and [32..63] -> [0..31]
  84     frac_dx = (value & 0x03fff800) >> 11;
  85
  86     value = tab_log2[frac_x0];
  87     value += (frac_dx * (tab_log2[frac_x0+1] - tab_log2[frac_x0])) >> 15;
  88
  89     return (power_int << 15) + value;
  90 }
  91
  92 int64_t ff_dot_product(const int16_t *a, const int16_t *b, int length)
  93 {
  94     int i;
  95     int64_t sum = 0;
  96
  97     for (i = 0; i < length; i++)
  98         sum += MUL16(a[i], b[i]);
  99
 100     return sum;
 101 }