[FFmpeg-devel] [PATCH] Common fixed-point ACELP routines (1/3) - math
Vladimir Voroshilov
voroshil
Wed Apr 23 19:58:08 CEST 2008
Michael Niedermayer wrote:
> On Tue, Apr 22, 2008 at 11:53:10PM +0700, Vladimir Voroshilov wrote:
> >
> > Michael Niedermayer wrote:
> > > On Tue, Apr 22, 2008 at 09:12:16AM +0700, Vladimir Voroshilov wrote:
> > > > On Tue, Apr 22, 2008 at 6:05 AM, Michael Niedermayer <michaelni at gmx.at> wrote:
[...]
> > > Sadly this is not bitexact.
> > > The question here is, is it possible at all to maintain bitexactness with
> > > common code?
> > > I assume the other acelp codecs use slightly different integer
> > > implementations?
> >
> > This is rhetorical question :) ?
> >
> > I prefer to keep bitexact code at least till commit into FFmpeg tree.
> > This makes developments much simple and safe in terms of breaking anything.
> > Of course i can keep bitexact routines in local tree only and commit more
> > precise math operations to tree.
> > Thus if new code gives acceptable results i'll not have anything against it.
>
> We could keep bitexact routines under #ifdef for each acelp variant if you
> want. That way there are both well written fast+precisse ones and the
> g729 bitexact ones. That surely could come in handy if we stumble across a
> g729 stream which decodes with artifacts.
I've put old code under #ifdef G729_BITAXACT
Don't sure that this is good solution, since it will lead to brokenness another
codecs sharing this code.
I'd prefer see working "+bitexact" option but i'm afraid such support will cause code
duplication (both bitexact and precise routines should be enabled).
[...]
> > int16_t ff_acelp_cos(uint16_t arg)
>
> ff_cos()
Fixed.
[...]
> > int ff_acelp_exp2(uint16_t power)
>
> ff_exp2()
Fixed.
[...]
> >
> > result= (result<<4) + ((result*exp2b[(power>>5)&31])>>16);
> > result= result + ((result*exp2c[ power &31])>>18);
>
> > return (result + 32) >> 6;
>
> IMHO the rounding should be done outside exp2(), this is just senslessly
> throwing bits away.
Fixed.
Also doxygen comments moved to header.
--
Regards,
Vladimir Voroshilov mailto:voroshil at gmail.com
Omsk State University
JID: voroshil at jabber.ru
ICQ: 95587719
-------------- next part --------------
diff --git a/libavcodec/acelp_math.c b/libavcodec/acelp_math.c
new file mode 100644
index 0000000..ab05f6b
--- /dev/null
+++ b/libavcodec/acelp_math.c
@@ -0,0 +1,208 @@
+/*
+ * Various fixed-point math operations
+ *
+ * Copyright (c) 2008 Vladimir Voroshilov
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include <inttypes.h>
+#include <limits.h>
+#include <assert.h>
+
+#include "avcodec.h"
+#include "acelp_math.h"
+
+#ifdef G729_BITEXACT
+/**
+ * Cosine table: base_cos[i] = (1<<15) * cos(i*PI/64)
+ */
+static const int16_t base_cos[64] =
+{
+ 32767, 32729, 32610, 32413, 32138, 31786, 31357, 30853,
+ 30274, 29622, 28899, 28106, 27246, 26320, 25330, 24279,
+ 23170, 22006, 20788, 19520, 18205, 16846, 15447, 14010,
+ 12540, 11039, 9512, 7962, 6393, 4808, 3212, 1608,
+ 0, -1608, -3212, -4808, -6393, -7962, -9512, -11039,
+ -12540, -14010, -15447, -16846, -18205, -19520, -20788, -22006,
+ -23170, -24279, -25330, -26320, -27246, -28106, -28899, -29622,
+ -30274, -30853, -31357, -31786, -32138, -32413, -32610, -32729
+};
+
+/**
+ * Slope used to compute cos(x)
+ *
+ * cos(ind*64+offset) = base_cos[ind]+offset*slope_cos[ind]
+ * values multiplied by 1<<19
+ */
+static const int16_t slope_cos[64] =
+{
+ -632, -1893, -3150, -4399, -5638, -6863, -8072, -9261,
+ -10428, -11570, -12684, -13767, -14817, -15832, -16808, -17744,
+ -18637, -19486, -20287, -21039, -21741, -22390, -22986, -23526,
+ -24009, -24435, -24801, -25108, -25354, -25540, -25664, -25726,
+ -25726, -25664, -25540, -25354, -25108, -24801, -24435, -24009,
+ -23526, -22986, -22390, -21741, -21039, -20287, -19486, -18637,
+ -17744, -16808, -15832, -14817, -13767, -12684, -11570, -10428,
+ -9261, -8072, -6863, -5638, -4399, -3150, -1893, -632
+};
+
+/**
+ * Table used to compute exp2(x)
+ *
+ * tab_exp2[i] = (1<<14) * exp2(i/32) = 2^(i/32) i=0..32
+ */
+static const uint16_t tab_exp2[33] =
+{
+ 16384, 16743, 17109, 17484, 17867, 18258, 18658, 19066, 19484, 19911,
+ 20347, 20792, 21247, 21713, 22188, 22674, 23170, 23678, 24196, 24726,
+ 25268, 25821, 26386, 26964, 27554, 28158, 28774, 29405, 30048, 30706,
+ 31379, 32066, 32767
+};
+
+int16_t ff_cos(uint16_t arg)
+{
+ uint8_t offset= arg;
+ uint8_t ind = arg >> 8;
+
+ assert(arg < 0x4000);
+
+ return FFMAX(base_cos[ind] + ((slope_cos[ind] * offset) >> 12), -0x8000);
+}
+
+int ff_exp2(uint16_t power)
+{
+ uint16_t frac_x0;
+ uint16_t frac_dx;
+ int result;
+
+ assert(power <= 0x7fff);
+
+ frac_x0 = power >> 10;
+ frac_dx = (power & 0x03ff) << 5;
+
+ result = tab_exp2[frac_x0] << 15;
+ result += frac_dx * (tab_exp2[frac_x0+1] - tab_exp2[frac_x0]);
+
+ return result >> 9;
+}
+
+#else // G729_BITEXACT
+
+/**
+ * Cosine table: base_cos[i] = (1<<15) * cos(i*PI/64)
+ */
+static const int16_t tab_cos[65] =
+{
+ 32767, 32738, 32617, 32421, 32145, 31793, 31364, 30860,
+ 30280, 29629, 28905, 28113, 27252, 26326, 25336, 24285,
+ 23176, 22011, 20793, 19525, 18210, 16851, 15451, 14014,
+ 12543, 11043, 9515, 7965, 6395, 4810, 3214, 1609,
+ 1, -1607, -3211, -4808, -6393, -7962, -9513, -11040,
+ -12541, -14012, -15449, -16848, -18207, -19523, -20791, -22009,
+ -23174, -24283, -25334, -26324, -27250, -28111, -28904, -29627,
+ -30279, -30858, -31363, -31792, -32144, -32419, -32616, -32736, -32768,
+};
+
+static const uint16_t exp2a[] =
+{
+ 0, 1435, 2902, 4400, 5932, 7496, 9096, 10730,
+ 12400, 14106, 15850, 17633, 19454, 21315, 23216, 25160,
+ 27146, 29175, 31249, 33369, 35534, 37747, 40009, 42320,
+ 44682, 47095, 49562, 52082, 54658, 57289, 59979, 62727,
+};
+
+static const uint16_t exp2b[]={
+ 0, 710, 1421, 2132, 2843, 3555, 4267, 4980,
+ 5694, 6408, 7122, 7837, 8552, 9268, 9984, 10701,
+ 11418, 12136, 12854, 13573, 14292, 15012, 15732, 16453,
+ 17174, 17896, 18618, 19340, 20063, 20787, 21511, 22236,
+};
+
+static const uint8_t exp2c[] =
+{
+ 0, 6, 11, 17, 22, 28, 33, 39,
+ 44, 50, 55, 61, 67, 72, 78, 83,
+ 89, 94, 100, 105, 111, 116, 122, 128,
+ 133, 139, 144, 150, 155, 161, 166, 172,
+};
+
+int16_t ff_cos(uint16_t arg)
+{
+ uint8_t offset= arg;
+ uint8_t ind = arg >> 8;
+
+ assert(arg <= 0x3fff);
+
+ return tab_cos[ind] + (offset * (tab_cos[ind+1] - tab_cos[ind]) >> 8);
+}
+
+int ff_exp2(uint16_t power)
+{
+ unsigned int result= exp2a[power>>10] + 0x10000;
+
+ assert(arg <= 0x7fff);
+
+ result= (result<<4) + ((result*exp2b[(power>>5)&31])>>16);
+ result= result + ((result*exp2c[ power &31])>>18);
+ return result;
+}
+
+#endif // else G729_BITEXACT
+
+/**
+ * Table used to compute log2(x)
+ *
+ * tab_log2[i] = (1<<15) * log2(1 + i/32), i=0..32
+ */
+static const uint16_t tab_log2[33] =
+{
+#ifdef G729_BITEXACT
+ 0, 1455, 2866, 4236, 5568, 6863, 8124, 9352,
+ 10549, 11716, 12855, 13967, 15054, 16117, 17156, 18172,
+ 19167, 20142, 21097, 22033, 22951, 23852, 24735, 25603,
+ 26455, 27291, 28113, 28922, 29716, 30497, 31266, 32023, 32767,
+#else
+ 4, 1459, 2870, 4240, 5572, 6867, 8127, 9355,
+ 10552, 11719, 12858, 13971, 15057, 16120, 17158, 18175,
+ 19170, 20145, 21100, 22036, 22954, 23854, 24738, 25605,
+ 26457, 27294, 28116, 28924, 29719, 30500, 31269, 32025, 32769,
+#endif
+};
+
+int ff_log2(int value)
+{
+ uint32_t result;
+ uint8_t power_int;
+ uint8_t frac_x0;
+ uint16_t frac_dx;
+
+ assert(value > 0);
+
+ // Stripping zeros from beginning
+ power_int = av_log2(value);
+ result = value << (31 - power_int);
+
+ // b31 is always non-zero now
+ frac_x0 = (result & 0x7c000000) >> 26; // b26-b31 and [32..63] -> [0..31]
+ frac_dx = (result & 0x03fff800) >> 11;
+
+ result = tab_log2[frac_x0] << 15;
+ result += frac_dx * (tab_log2[frac_x0+1] - tab_log2[frac_x0]);
+
+ return (power_int << 15) + (result >> 15);
+}
diff --git a/libavcodec/acelp_math.h b/libavcodec/acelp_math.h
new file mode 100644
index 0000000..7c6b36f
--- /dev/null
+++ b/libavcodec/acelp_math.h
@@ -0,0 +1,111 @@
+/*
+ * Various fixed-point math operations
+ *
+ * Copyright (c) 2008 Vladimir Voroshilov
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef FFMPEG_ACELP_MATH_H
+#define FFMPEG_ACELP_MATH_H
+
+/**
+ * \brief fixed-point implementation of cosine in [0; PI) domain
+ * \param arg fixed-point cosine argument, 0 <= arg < 0x4000
+ *
+ * \return value of (1<<15) * cos(arg * PI / (1<<14)), -0x8000 <= result <= 0x7fff
+ */
+int16_t ff_cos(uint16_t arg);
+
+/**
+ * \brief fixed-point implementation of exp2(x) in [0; 1] domain
+ * \param power argument to exp2, 0 <= power <= 0x7fff
+ *
+ * \return value of (1<<20) * exp2(power / (1<<15))
+ * 0x100000 <= result <= 0x1fffcd
+ */
+int ff_exp2(uint16_t power);
+
+/**
+ * \brief Calculates log2(x)
+ * \param value function argument, 0 < value <= 7fff ffff
+ *
+ * \return value of (1<<15) * log2(value)
+ */
+int ff_log2(int value);
+
+/**
+ * \brief multiplies 32-bit integer by another 16-bit and divides result by 2^15
+ * \param var_q24 32-bit integer
+ * \param var_15 16-bit integer
+ *
+ * \return result of (var_q24 * var_q15 >> 15) with clipping to [INT_MIN; INT_MAX] range
+ */
+static inline int mul_24_15(int var_q24, int16_t var_q15)
+{
+ int64_t tmp = (((int64_t)var_q24 * (int64_t)var_q15) >> 15);
+
+ if(tmp < INT_MIN)
+ return INT_MIN;
+ else if (tmp > INT_MAX)
+ return INT_MAX;
+ else
+ return tmp;
+}
+
+/**
+ * \brief Calculates sum of array elements multiplications
+ * \param speech array with input data
+ * \param cycles number elements to proceed
+ * \param offset offset for calculation sum of s[i]*s[i+offset]
+ * \param shift right shift by this value will be done before multiplication
+ *
+ * \return sum of multiplications
+ *
+ * \note array must be at least length+offset long!
+ */
+static int sum_of_squares(const int16_t* speech, int cycles, int offset, int shift)
+{
+ int n;
+ int sum = 0;
+
+ for(n=0; n<cycles; n++)
+ sum += (speech[n] >> shift) * (speech[n + offset] >> shift);
+
+ return av_clip(sum, -0x40000000, 0x3fffffff);
+}
+
+/**
+ * \brief Calculates sum of array elements absolute values
+ * \param speech array with input data
+ * \param cycles number elements to proceed
+ * \param shift right shift by this value will be done before addition
+ *
+ * \return sum of absolute values
+ */
+static int sum_of_absolute(const int16_t* speech, int cycles, int shift)
+{
+ int n;
+ int sum = 0;
+
+ for(n=0; n<cycles; n++)
+ sum += FFABS(speech[n] >> shift);
+
+ return sum;
+}
+
+#endif // FFMPEG_ACELP_MATH_H
More information about the ffmpeg-devel
mailing list