[FFmpeg-devel] [PATCH] Common fixed-point ACELP routines (1/3) - math

[Michael Niedermayer]

On Tue, Apr 22, 2008 at 09:12:16AM +0700, Vladimir Voroshilov wrote:
[...]
> +/**
> + * Table used to compute exp2(x)
> + *
> + * ff_acelp_tab_exp2[i] = (1<<14) * exp2(i/32) = 2^(i/32) i=0..32
> + */
> +static const uint16_t ff_acelp_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
> +};
> +

[...]
> +/**
> + * \brief fixed-point implementation of exp2(x) in [0; 1] domain
> + * \param power argument to exp2, 0 <= power <= 0x7fff
> + *
> + * \return value of (1<<14) * exp2(power / (1<<15)) rounded to nearest, 0x4000 <= result <= 0x7fff
> + */
> +int ff_acelp_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 = ff_acelp_tab_exp2[frac_x0] << 15;
> +    result += frac_dx * (ff_acelp_tab_exp2[frac_x0+1] - ff_acelp_tab_exp2[frac_x0]);
> +
> +    return (result + 0x4000) >> 15;
> +}

heres a more accurate exp2()

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,
};

int myexp2b(uint16_t power){
    unsigned int result= exp2a[power>>10] + 0x10000;

    result= (result<<4) + ((result*exp2b[(power>>5)&31])>>16);
    return result + ((result*exp2c[ power&31 ])>>18);
}

Aboves max error is 1/10 and the sum of squared errors is 1/100 of your
code.

Note, the output has more bits than your function thus it needs to be
scaled down.

[...]
-- 
Michael     GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB

I hate to see young programmers poisoned by the kind of thinking
Ulrich Drepper puts forward since it is simply too narrow -- Roman Shaposhnik
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Digital signature
URL: <http://lists.mplayerhq.hu/pipermail/ffmpeg-devel/attachments/20080422/edad0faa/attachment.pgp>