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

[Michael Niedermayer]

On Thu, Apr 24, 2008 at 12:58:08AM +0700, Vladimir Voroshilov wrote:
> 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.
[...]
> +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);
> +}

int ff_log2(unsigned int value)
{
    uint8_t  power_int;
    uint8_t  frac_x0;
    uint16_t frac_dx;

    assert(value > 0);

    power_int = av_log2(value);
    value <<= 31 - power_int;

    frac_x0 = (value & 0x7c000000) >> 26
    frac_dx = (value & 0x03fff800) >> 11;

    value = tab_log2[frac_x0];
    value += (frac_dx * (tab_log2[frac_x0+1] - tab_log2[frac_x0])) >> 15;

    return (power_int << 15) + 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;
> +}

INT_MIN/MAX depends on how bits there are in int, this is not what g729
wants.


> +

> +/**
> + * \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;
> +}

The place where the shift is done is not optimal for precission.

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

The greatest way to live with honor in this world is to be what we pretend
to be. -- Socrates
-------------- 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/20080423/e4b2bbc0/attachment.pgp>