vbrquantize.c

音频编码

/*
 *	MP3 quantization
 *
 *	Copyright (c) 1999-2000 Mark Taylor
 *	Copyright (c) 2000-2005 Robert Hegemann
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/* $Id: vbrquantize.c,v 1.103.2.1 2005/11/20 14:08:25 bouvigne Exp $ */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

#include <assert.h>
#include "util.h"
#include "l3side.h"
#include "quantize_pvt.h"
#include "vbrquantize.h"

#ifdef WITH_DMALLOC
#  include <dmalloc.h>
#endif



struct algo_s;
typedef struct algo_s algo_t;

typedef void (*quantize_f) (const algo_t *);

typedef int (*find_f) (const FLOAT *, const FLOAT *, FLOAT, int, int);

typedef int (*alloc_sf_f) (const algo_t *, int *, const int *, int);

struct algo_s {
    find_f  find;
    quantize_f quantize;
    alloc_sf_f alloc;
    const FLOAT *xr34orig;
    lame_internal_flags *gfc;
    gr_info *cod_info;
    int mingain_l;
    int mingain_s[3];
};



/*  Remarks on optimizing compilers:
 *
 *  the MSVC compiler may get into aliasing problems when accessing
 *  memory through the fi_union. declaring it volatile does the trick here
 *
 *  the calc_sfb_noise_* functions are not inlined because the intel compiler
 *  optimized executeables won't work as expected anymore
 */

#ifdef _MSC_VER
#  define VOLATILE volatile
#else
#  define VOLATILE
#endif

typedef VOLATILE union {
    float   f;
    int     i;
} fi_union;



#define DOUBLEX double

#define MAGIC_FLOAT_def (65536*(128))
#define MAGIC_INT_def    0x4b000000

#ifdef TAKEHIRO_IEEE754_HACK
#  define ROUNDFAC_def -0.0946f
#else
/*********************************************************************
 * XRPOW_FTOI is a macro to convert floats to ints.
 * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x]
 *                                         ROUNDFAC= -0.0946
 *
 * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]
 *                                   ROUNDFAC=0.4054
 *********************************************************************/
#  define QUANTFAC(rx)  adj43[rx]
#  define ROUNDFAC_def 0.4054f
#  define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
#endif

static int const MAGIC_INT = MAGIC_INT_def;
static DOUBLEX const ROUNDFAC = ROUNDFAC_def;
static DOUBLEX const MAGIC_FLOAT = (65536 * (128));
static DOUBLEX const ROUNDFAC_plus_MAGIC_FLOAT = ROUNDFAC_def + MAGIC_FLOAT_def;



static int
valid_sf(int sf)
{
    return (sf >= 0 ? (sf <= 255 ? sf : 255) : 0);
}



static  FLOAT
max_x34(const FLOAT * xr34, unsigned int bw)
{
    FLOAT   xfsf = 0;
    int     j = bw >> 1;
    int     remaining = j % 2;
    assert(bw >= 0);
    for (j >>= 1; j > 0; --j) {
        if (xfsf < xr34[0]) {
            xfsf = xr34[0];
        }
        if (xfsf < xr34[1]) {
            xfsf = xr34[1];
        }
        if (xfsf < xr34[2]) {
            xfsf = xr34[2];
        }
        if (xfsf < xr34[3]) {
            xfsf = xr34[3];
        }
        xr34 += 4;
    }
    if (remaining) {
        if (xfsf < xr34[0]) {
            xfsf = xr34[0];
        }
        if (xfsf < xr34[1]) {
            xfsf = xr34[1];
        }
    }
    return xfsf;
}



static int
find_lowest_scalefac(const FLOAT xr34)
{
    FLOAT   xfsf;
    int     sf = 128, sf_ok = 10000, delsf = 128, i;
    for (i = 0; i < 8; ++i) {
        delsf >>= 1;
        xfsf = IPOW20(sf) * xr34;
        if (xfsf <= IXMAX_VAL) {
            sf_ok = sf;
            sf -= delsf;
        }
        else {
            sf += delsf;
        }
    }
    if (sf_ok < 255) {
        sf = sf_ok;
    }
    return sf;
}



static void
k_34_4(DOUBLEX x[4], int l3[4])
{
#ifdef TAKEHIRO_IEEE754_HACK
    fi_union fi[4];

    assert(x[0] <= IXMAX_VAL && x[1] <= IXMAX_VAL && x[2] <= IXMAX_VAL && x[3] <= IXMAX_VAL);
    x[0] += MAGIC_FLOAT;
    fi[0].f = x[0];
    x[1] += MAGIC_FLOAT;
    fi[1].f = x[1];
    x[2] += MAGIC_FLOAT;
    fi[2].f = x[2];
    x[3] += MAGIC_FLOAT;
    fi[3].f = x[3];
    fi[0].f = x[0] + adj43asm[fi[0].i - MAGIC_INT];
    fi[1].f = x[1] + adj43asm[fi[1].i - MAGIC_INT];
    fi[2].f = x[2] + adj43asm[fi[2].i - MAGIC_INT];
    fi[3].f = x[3] + adj43asm[fi[3].i - MAGIC_INT];
    l3[0] = fi[0].i - MAGIC_INT;
    l3[1] = fi[1].i - MAGIC_INT;
    l3[2] = fi[2].i - MAGIC_INT;
    l3[3] = fi[3].i - MAGIC_INT;
#else
    assert(x[0] <= IXMAX_VAL && x[1] <= IXMAX_VAL && x[2] <= IXMAX_VAL && x[3] <= IXMAX_VAL);
    XRPOW_FTOI(x[0], l3[0]);
    XRPOW_FTOI(x[1], l3[1]);
    XRPOW_FTOI(x[2], l3[2]);
    XRPOW_FTOI(x[3], l3[3]);
    x[0] += QUANTFAC(l3[0]);
    x[1] += QUANTFAC(l3[1]);
    x[2] += QUANTFAC(l3[2]);
    x[3] += QUANTFAC(l3[3]);
    XRPOW_FTOI(x[0], l3[0]);
    XRPOW_FTOI(x[1], l3[1]);
    XRPOW_FTOI(x[2], l3[2]);
    XRPOW_FTOI(x[3], l3[3]);
#endif
}



static void
k_34_2(DOUBLEX x[2], int l3[2])
{
#ifdef TAKEHIRO_IEEE754_HACK
    fi_union fi[2];

    assert(x[0] <= IXMAX_VAL && x[1] <= IXMAX_VAL);
    x[0] += MAGIC_FLOAT;
    fi[0].f = x[0];
    x[1] += MAGIC_FLOAT;
    fi[1].f = x[1];
    fi[0].f = x[0] + adj43asm[fi[0].i - MAGIC_INT];
    fi[1].f = x[1] + adj43asm[fi[1].i - MAGIC_INT];
    l3[0] = fi[0].i - MAGIC_INT;
    l3[1] = fi[1].i - MAGIC_INT;
#else
    assert(x[0] <= IXMAX_VAL && x[1] <= IXMAX_VAL);
    XRPOW_FTOI(x[0], l3[0]);
    XRPOW_FTOI(x[1], l3[1]);
    x[0] += QUANTFAC(l3[0]);
    x[1] += QUANTFAC(l3[1]);
    XRPOW_FTOI(x[0], l3[0]);
    XRPOW_FTOI(x[1], l3[1]);
#endif
}



static void
k_iso_4(DOUBLEX x[4], int l3[4])
{
#ifdef TAKEHIRO_IEEE754_HACK
    fi_union fi[4];

    assert(x[0] <= IXMAX_VAL && x[1] <= IXMAX_VAL && x[2] <= IXMAX_VAL && x[3] <= IXMAX_VAL);
    x[0] += ROUNDFAC_plus_MAGIC_FLOAT;
    fi[0].f = x[0];
    x[1] += ROUNDFAC_plus_MAGIC_FLOAT;
    fi[1].f = x[1];
    x[2] += ROUNDFAC_plus_MAGIC_FLOAT;
    fi[2].f = x[2];
    x[3] += ROUNDFAC_plus_MAGIC_FLOAT;
    fi[3].f = x[3];
    l3[0] = fi[0].i - MAGIC_INT;
    l3[1] = fi[1].i - MAGIC_INT;
    l3[2] = fi[2].i - MAGIC_INT;
    l3[3] = fi[3].i - MAGIC_INT;
#else
    l3[0] = x[0] + ROUNDFAC;
    l3[1] = x[1] + ROUNDFAC;
    l3[2] = x[2] + ROUNDFAC;
    l3[3] = x[3] + ROUNDFAC;
#endif
}



static void
k_iso_2(DOUBLEX x[2], int l3[2])
{
#ifdef TAKEHIRO_IEEE754_HACK
    fi_union fi[2];

    assert(x[0] <= IXMAX_VAL && x[1] <= IXMAX_VAL);
    x[0] += ROUNDFAC_plus_MAGIC_FLOAT;
    fi[0].f = x[0];
    x[1] += ROUNDFAC_plus_MAGIC_FLOAT;
    fi[1].f = x[1];
    l3[0] = fi[0].i - MAGIC_INT;
    l3[1] = fi[1].i - MAGIC_INT;
#else
    l3[0] = x[0] + ROUNDFAC;
    l3[1] = x[1] + ROUNDFAC;
#endif
}



/*  do call the calc_sfb_noise_* functions only with sf values
 *  for which holds: sfpow34*xr34 <= IXMAX_VAL
 */

static  FLOAT
calc_sfb_noise_x34(const FLOAT * xr, const FLOAT * xr34, unsigned int bw, int sf)
{
    DOUBLEX x[4];
    int     l3[4];
    const int SF = valid_sf(sf);
    const FLOAT sfpow = POW20(SF); /*pow(2.0,sf/4.0); */
    const FLOAT sfpow34 = IPOW20(SF); /*pow(sfpow,-3.0/4.0); */

    FLOAT   xfsf = 0;
    int     j = bw >> 1;
    int     remaining = j % 2;
    assert(bw >= 0);
    for (j >>= 1; j > 0; --j) {
        x[0] = sfpow34 * xr34[0];
        x[1] = sfpow34 * xr34[1];
        x[2] = sfpow34 * xr34[2];
        x[3] = sfpow34 * xr34[3];

        k_34_4(x, l3);

        x[0] = fabs(xr[0]) - sfpow * pow43[l3[0]];
        x[1] = fabs(xr[1]) - sfpow * pow43[l3[1]];
        x[2] = fabs(xr[2]) - sfpow * pow43[l3[2]];
        x[3] = fabs(xr[3]) - sfpow * pow43[l3[3]];
        xfsf += (x[0] * x[0] + x[1] * x[1]) + (x[2] * x[2] + x[3] * x[3]);

        xr += 4;
        xr34 += 4;
    }
    if (remaining) {
        x[0] = sfpow34 * xr34[0];
        x[1] = sfpow34 * xr34[1];

        k_34_2(x, l3);

        x[0] = fabs(xr[0]) - sfpow * pow43[l3[0]];
        x[1] = fabs(xr[1]) - sfpow * pow43[l3[1]];
        xfsf += x[0] * x[0] + x[1] * x[1];
    }
    return xfsf;
}



static  FLOAT
calc_sfb_noise_ISO(const FLOAT * xr, const FLOAT * xr34, unsigned int bw, int sf)
{
    DOUBLEX x[4];
    int     l3[4];
    const int SF = valid_sf(sf);
    const FLOAT sfpow = POW20(SF); /*pow(2.0,sf/4.0); */
    const FLOAT sfpow34 = IPOW20(SF); /*pow(sfpow,-3.0/4.0); */

    FLOAT   xfsf = 0;
    int     j = bw >> 1;
    int     remaining = j % 2;
    assert(bw >= 0);
    for (j >>= 1; j > 0; --j) {
        x[0] = sfpow34 * xr34[0];
        x[1] = sfpow34 * xr34[1];
        x[2] = sfpow34 * xr34[2];
        x[3] = sfpow34 * xr34[3];

        k_iso_4(x, l3);

        x[0] = fabs(xr[0]) - sfpow * pow43[l3[0]];
        x[1] = fabs(xr[1]) - sfpow * pow43[l3[1]];
        x[2] = fabs(xr[2]) - sfpow * pow43[l3[2]];
        x[3] = fabs(xr[3]) - sfpow * pow43[l3[3]];

        xfsf += (x[0] * x[0] + x[1] * x[1]) + (x[2] * x[2] + x[3] * x[3]);

        xr += 4;
        xr34 += 4;
    }
    if (remaining) {
        x[0] = sfpow34 * xr34[0];
        x[1] = sfpow34 * xr34[1];

        k_iso_2(x, l3);

        x[0] = fabs(xr[0]) - sfpow * pow43[l3[0]];
        x[1] = fabs(xr[1]) - sfpow * pow43[l3[1]];
        xfsf += x[0] * x[0] + x[1] * x[1];
    }
    return xfsf;
}



/* the find_scalefac* routines calculate
 * a quantization step size which would
 * introduce as much noise as is allowed.
 * The larger the step size the more
 * quantization noise we'll get. The
 * scalefactors are there to lower the
 * global step size, allowing limited
 * differences in quantization step sizes
 * per band (shaping the noise).
 */

static int
find_scalefac_x34(const FLOAT * xr, const FLOAT * xr34, FLOAT l3_xmin, int bw, int sf_min)
{
    int     sf = 128, sf_ok = 10000, delsf = 128, i;
    for (i = 0; i < 8; ++i) {
        delsf >>= 1;
        if (sf <= sf_min) {
            sf += delsf;
        }
        else {
            if ((sf < 255 && calc_sfb_noise_x34(xr, xr34, bw, sf + 1) > l3_xmin)
                || calc_sfb_noise_x34(xr, xr34, bw, sf) > l3_xmin
                || calc_sfb_noise_x34(xr, xr34, bw, sf - 1) > l3_xmin) {
                /* distortion.  try a smaller scalefactor */
                sf -= delsf;
            }
            else {
                sf_ok = sf;
                sf += delsf;
            }
        }
    }
    /*  returning a scalefac without distortion, if possible
     */
    if (sf_ok <= 255) {
        sf = sf_ok;
    }
    return sf;
}



static int
find_scalefac_ISO(const FLOAT * xr, const FLOAT * xr34, FLOAT l3_xmin, int bw, int sf_min)
{
    int     sf = 128, sf_ok = 10000, delsf = 128, i;
    for (i = 0; i < 8; ++i) {
        delsf >>= 1;
        if (sf <= sf_min) {
            sf += delsf;
        }
        else {
            if ((sf < 255 && calc_sfb_noise_ISO(xr, xr34, bw, sf + 1) > l3_xmin)
                || calc_sfb_noise_ISO(xr, xr34, bw, sf) > l3_xmin
                || calc_sfb_noise_ISO(xr, xr34, bw, sf - 1) > l3_xmin) {
                /* distortion.  try a smaller scalefactor */
                sf -= delsf;
            }
            else {
                sf_ok = sf;
                sf += delsf;