quantize_pvt.c

音频编码

/*
 *	quantize_pvt source file
 *
 *	Copyright (c) 1999-2002 Takehiro Tominaga
 *	Copyright (c) 2000-2002 Robert Hegemann
 *	Copyright (c) 2001 Naoki Shibata
 *	Copyright (c) 2002-2005 Gabriel Bouvigne
 *
 * 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: quantize_pvt.c,v 1.131.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 "lame-analysis.h"
#include "tables.h"
#include "reservoir.h"
#include "quantize_pvt.h"

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


#define NSATHSCALE 100 /* Assuming dynamic range=96dB, this value should be 92 */

/*
  The following table is used to implement the scalefactor
  partitioning for MPEG2 as described in section
  2.4.3.2 of the IS. The indexing corresponds to the
  way the tables are presented in the IS:

  [table_number][row_in_table][column of nr_of_sfb]
*/
const int  nr_of_sfb_block [6] [3] [4] =
{
  {
    {6, 5, 5, 5},
    {9, 9, 9, 9},
    {6, 9, 9, 9}
  },
  {
    {6, 5, 7, 3},
    {9, 9, 12, 6},
    {6, 9, 12, 6}
  },
  {
    {11, 10, 0, 0},
    {18, 18, 0, 0},
    {15,18,0,0}
  },
  {
    {7, 7, 7, 0},
    {12, 12, 12, 0},
    {6, 15, 12, 0}
  },
  {
    {6, 6, 6, 3},
    {12, 9, 9, 6},
    {6, 12, 9, 6}
  },
  {
    {8, 8, 5, 0},
    {15,12,9,0},
    {6,18,9,0}
  }
};


/* Table B.6: layer3 preemphasis */
const int  pretab [SBMAX_l] =
{
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    1, 1, 1, 1, 2, 2, 3, 3, 3, 2, 0
};

/*
  Here are MPEG1 Table B.8 and MPEG2 Table B.1
  -- Layer III scalefactor bands. 
  Index into this using a method such as:
    idx  = fr_ps->header->sampling_frequency
           + (fr_ps->header->version * 3)
*/


const scalefac_struct sfBandIndex[9] =
{
  { /* Table B.2.b: 22.05 kHz */
    {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
    {0,4,8,12,18,24,32,42,56,74,100,132,174,192}
  },
  { /* Table B.2.c: 24 kHz */                 /* docs: 332. mpg123(broken): 330 */
    {0,6,12,18,24,30,36,44,54,66,80,96,114,136,162,194,232,278, 332, 394,464,540,576},
    {0,4,8,12,18,26,36,48,62,80,104,136,180,192}
  },
  { /* Table B.2.a: 16 kHz */
    {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
    {0,4,8,12,18,26,36,48,62,80,104,134,174,192}
  },
  { /* Table B.8.b: 44.1 kHz */
    {0,4,8,12,16,20,24,30,36,44,52,62,74,90,110,134,162,196,238,288,342,418,576},
    {0,4,8,12,16,22,30,40,52,66,84,106,136,192}
  },
  { /* Table B.8.c: 48 kHz */
    {0,4,8,12,16,20,24,30,36,42,50,60,72,88,106,128,156,190,230,276,330,384,576},
    {0,4,8,12,16,22,28,38,50,64,80,100,126,192}
  },
  { /* Table B.8.a: 32 kHz */
    {0,4,8,12,16,20,24,30,36,44,54,66,82,102,126,156,194,240,296,364,448,550,576},
    {0,4,8,12,16,22,30,42,58,78,104,138,180,192}
  },
  { /* MPEG-2.5 11.025 kHz */
    {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
    {0/3,12/3,24/3,36/3,54/3,78/3,108/3,144/3,186/3,240/3,312/3,402/3,522/3,576/3}
  },
  { /* MPEG-2.5 12 kHz */
    {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
    {0/3,12/3,24/3,36/3,54/3,78/3,108/3,144/3,186/3,240/3,312/3,402/3,522/3,576/3}
  },
  { /* MPEG-2.5 8 kHz */
    {0,12,24,36,48,60,72,88,108,132,160,192,232,280,336,400,476,566,568,570,572,574,576},
    {0/3,24/3,48/3,72/3,108/3,156/3,216/3,288/3,372/3,480/3,486/3,492/3,498/3,576/3}
  }
};



FLOAT pow20[Q_MAX+Q_MAX2+1];
FLOAT ipow20[Q_MAX];
FLOAT pow43[PRECALC_SIZE];
/* initialized in first call to iteration_init */
#ifdef TAKEHIRO_IEEE754_HACK
FLOAT adj43asm[PRECALC_SIZE];
#else
FLOAT adj43[PRECALC_SIZE];
#endif

/* 
compute the ATH for each scalefactor band 
cd range:  0..96db

Input:  3.3kHz signal  32767 amplitude  (3.3kHz is where ATH is smallest = -5db)
longblocks:  sfb=12   en0/bw=-11db    max_en0 = 1.3db
shortblocks: sfb=5           -9db              0db

Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
longblocks:  amp=1      sfb=12   en0/bw=-103 db      max_en0 = -92db
            amp=32767   sfb=12           -12 db                 -1.4db 

Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
shortblocks: amp=1      sfb=5   en0/bw= -99                    -86 
            amp=32767   sfb=5           -9  db                  4db 


MAX energy of largest wave at 3.3kHz = 1db
AVE energy of largest wave at 3.3kHz = -11db
Let's take AVE:  -11db = maximum signal in sfb=12.  
Dynamic range of CD: 96db.  Therefor energy of smallest audible wave 
in sfb=12  = -11  - 96 = -107db = ATH at 3.3kHz.  

ATH formula for this wave: -5db.  To adjust to LAME scaling, we need
ATH = ATH_formula  - 103  (db)
ATH = ATH * 2.5e-10      (ener)

*/

static FLOAT ATHmdct( lame_global_flags *gfp, FLOAT f )
{
    FLOAT ath;
  
    ath = ATHformula( f , gfp );
	  
    if (gfp->psymodel == PSY_NSPSYTUNE) {
        ath -= NSATHSCALE;
    } else {
        ath -= 114;
    }
    
    /* modify the MDCT scaling for the ATH and convert to energy */
    ath = pow( 10.0, ath/10.0 + gfp->ATHlower);
    return ath;
}

static void compute_ath( lame_global_flags *gfp )
{
    FLOAT *ATH_l = gfp->internal_flags->ATH->l;
    FLOAT *ATH_psfb21 = gfp->internal_flags->ATH->psfb21;
    FLOAT *ATH_s = gfp->internal_flags->ATH->s;
    FLOAT *ATH_psfb12 = gfp->internal_flags->ATH->psfb12;
    lame_internal_flags *gfc = gfp->internal_flags;
    int sfb, i, start, end;
    FLOAT ATH_f;
    FLOAT samp_freq = gfp->out_samplerate;

    for (sfb = 0; sfb < SBMAX_l; sfb++) {
        start = gfc->scalefac_band.l[ sfb ];
        end   = gfc->scalefac_band.l[ sfb+1 ];
        ATH_l[sfb]=FLOAT_MAX;
        for (i = start ; i < end; i++) {
            FLOAT freq = i*samp_freq/(2*576);
            ATH_f = ATHmdct( gfp, freq );  /* freq in kHz */
            ATH_l[sfb] = Min( ATH_l[sfb], ATH_f );
        }
	    if (gfp->psymodel == PSY_GPSYCHO)
	        ATH_l[sfb] *=
		        (gfc->scalefac_band.l[sfb+1] - gfc->scalefac_band.l[sfb]);
    }

    for (sfb = 0; sfb < PSFB21; sfb++) {
        start = gfc->scalefac_band.psfb21[ sfb ];
        end = gfc->scalefac_band.psfb21[ sfb+1 ];
        ATH_psfb21[sfb]=FLOAT_MAX;
        for (i = start ; i < end; i++) {
            FLOAT freq = i*samp_freq/(2*576);
            ATH_f = ATHmdct( gfp, freq );  /* freq in kHz */
            ATH_psfb21[sfb] = Min( ATH_psfb21[sfb], ATH_f );
        }
    }

    for (sfb = 0; sfb < SBMAX_s; sfb++){
        start = gfc->scalefac_band.s[ sfb ];
        end   = gfc->scalefac_band.s[ sfb+1 ];
        ATH_s[sfb] = FLOAT_MAX;
        for (i = start ; i < end; i++) {
            FLOAT freq = i*samp_freq/(2*192);
            ATH_f = ATHmdct( gfp, freq );    /* freq in kHz */
            ATH_s[sfb] = Min( ATH_s[sfb], ATH_f );
        }
    	ATH_s[sfb] *=
	        (gfc->scalefac_band.s[sfb+1] - gfc->scalefac_band.s[sfb]);
    }

    for (sfb = 0; sfb < PSFB12; sfb++) {
        start = gfc->scalefac_band.psfb12[ sfb ];
        end = gfc->scalefac_band.psfb12[ sfb+1 ];
        ATH_psfb12[sfb]=FLOAT_MAX;
        for (i = start ; i < end; i++) {
            FLOAT freq = i*samp_freq/(2*192);
            ATH_f = ATHmdct( gfp, freq );  /* freq in kHz */
            ATH_psfb12[sfb] = Min( ATH_psfb12[sfb], ATH_f );
        }
        /*not sure about the following*/
        ATH_psfb12[sfb] *=
	        (gfc->scalefac_band.s[13] - gfc->scalefac_band.s[12]);
    }


    /*  no-ATH mode:
     *  reduce ATH to -200 dB
     */
    
    if (gfp->noATH) {
        for (sfb = 0; sfb < SBMAX_l; sfb++) {
            ATH_l[sfb] = 1E-37;
        }
        for (sfb = 0; sfb < PSFB21; sfb++) {
            ATH_psfb21[sfb] = 1E-37;
        }
        for (sfb = 0; sfb < SBMAX_s; sfb++) {
            ATH_s[sfb] = 1E-37;
        }
        for (sfb = 0; sfb < PSFB12; sfb++) {
            ATH_psfb12[sfb] = 1E-37;
        }
    }
    
    /*  work in progress, don't rely on it too much
     */
    gfc->ATH-> floor = 10. * log10( ATHmdct( gfp, -1. ) );
    
    /*
    {   FLOAT g=10000, t=1e30, x;
        for ( f = 100; f < 10000; f++ ) {
            x = ATHmdct( gfp, f );
            if ( t > x ) t = x, g = f;
        }
        printf("min=%g\n", g);
    }*/
}




/************************************************************************/
/*  initialization for iteration_loop */
/************************************************************************/
void
iteration_init( lame_global_flags *gfp)
{
  lame_internal_flags *gfc=gfp->internal_flags;
  III_side_info_t * const l3_side = &gfc->l3_side;
  int i;

  if ( gfc->iteration_init_init==0 ) {
    gfc->iteration_init_init=1;

    l3_side->main_data_begin = 0;
    compute_ath(gfp);

    pow43[0] = 0.0;
    for(i=1;i<PRECALC_SIZE;i++)
        pow43[i] = pow((FLOAT)i, 4.0/3.0);

#ifdef TAKEHIRO_IEEE754_HACK
    adj43asm[0] = 0.0;
    for (i = 1; i < PRECALC_SIZE; i++)
      adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]),0.75);
#else
    for (i = 0; i < PRECALC_SIZE-1; i++)
	adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
    adj43[i] = 0.5;
#endif
    for (i = 0; i < Q_MAX; i++)
	ipow20[i] = pow(2.0, (double)(i - 210) * -0.1875);
    for (i = 0; i <= Q_MAX+Q_MAX2; i++)
	pow20[i] = pow(2.0, (double)(i - 210 - Q_MAX2) * 0.25);

    huffman_init(gfc);
    quantize_init(gfc);
    init_xrpow_core_init(gfc);

    if (gfp->psymodel == PSY_NSPSYTUNE) {
	    FLOAT bass, alto, treble, sfb21;

        i = (gfp->exp_nspsytune >> 2) & 63;
        if (i >= 32)
            i -= 64;
        bass = pow(10, i / 4.0 / 10.0);

        i = (gfp->exp_nspsytune >> 8) & 63;
        if (i >= 32)
            i -= 64;
        alto = pow(10, i / 4.0 / 10.0);

        i = (gfp->exp_nspsytune >> 14) & 63;
        if (i >= 32)
            i -= 64;
        treble = pow(10, i / 4.0 / 10.0);

        /*  to be compatible with Naoki's original code, the next 6 bits
         *  define only the amount of changing treble for sfb21 */
        i = (gfp->exp_nspsytune >> 20) & 63;
        if (i >= 32)
            i -= 64;
        sfb21 = treble * pow(10, i / 4.0 / 10.0);

	for (i = 0; i < SBMAX_l; i++) {
	    FLOAT f;
	    if      (i <=  6) f = bass;
	    else if (i <= 13) f = alto;
	    else if (i <= 20) f = treble;
	    else              f = sfb21;

	    gfc->nsPsy.longfact[i] = f;
	}
	for (i = 0; i < SBMAX_s; i++) {
	    FLOAT f;
	    if      (i <=  5) f = bass;
	    else if (i <= 10) f = alto;
	    else if (i <= 11) f = treble;
	    else              f = sfb21;

	    gfc->nsPsy.shortfact[i] = f;
	}
    } else {
	for (i = 0; i < SBMAX_l; i++)
	    gfc->nsPsy.longfact[i] = 1.0;
	for (i = 0; i < SBMAX_s; i++)
	    gfc->nsPsy.shortfact[i] = 1.0;
    }
  }
}





/************************************************************************
 * allocate bits among 2 channels based on PE
 * mt 6/99
 * bugfixes rh 8/01: often allocated more than the allowed 4095 bits
 ************************************************************************/
int on_pe( lame_global_flags *gfp, FLOAT pe[][2], III_side_info_t *l3_side,
           int targ_bits[2], int mean_bits, int gr, int cbr )
{
    lame_internal_flags * gfc = gfp->internal_flags;
    gr_info *   cod_info;
    int     extra_bits, tbits, bits;
    int     add_bits[2]; 
    int     max_bits;  /* maximum allowed bits for this granule */
    int     ch;

    /* allocate targ_bits for granule */
    ResvMaxBits( gfp, mean_bits, &tbits, &extra_bits, cbr);
    max_bits = tbits + extra_bits;
    if (max_bits > MAX_BITS) /* hard limit per granule */
        max_bits = MAX_BITS;
    
    for ( bits = 0, ch = 0; ch < gfc->channels_out; ++ch ) {
        /******************************************************************
         * allocate bits for each channel 
         ******************************************************************/
        cod_info = &l3_side->tt[gr][ch];
    
        targ_bits[ch] = Min( MAX_BITS, tbits/gfc->channels_out );
    
        if (gfp->psymodel == PSY_NSPSYTUNE) {
            add_bits[ch] = targ_bits[ch] * pe[gr][ch] / 700.0 - targ_bits[ch];
        } 
        else {
            add_bits[ch] = (pe[gr][ch]-750) / 1.4;
            /* short blocks us a little extra, no matter what the pe */
            if ( cod_info->block_type == SHORT_TYPE ) {
	        if (add_bits[ch] < mean_bits/4) 
                    add_bits[ch] = mean_bits/4;
            }
        }

        /* at most increase bits by 1.5*average */
        if (add_bits[ch] > mean_bits*3/4) 
            add_bits[ch] = mean_bits*3/4;
        if (add_bits[ch] < 0) 
            add_bits[ch] = 0;

        if (add_bits[ch]+targ_bits[ch] > MAX_BITS) 
	    add_bits[ch] = Max( 0, MAX_BITS-targ_bits[ch] );

        bits += add_bits[ch];
    }
    if (bits > extra_bits) {
        for ( ch = 0; ch < gfc->channels_out; ++ch ) {
            add_bits[ch] = extra_bits * add_bits[ch] / bits;
        }
    }

    for ( ch = 0; ch < gfc->channels_out; ++ch ) {
        targ_bits[ch] += add_bits[ch];
        extra_bits    -= add_bits[ch];
        assert( targ_bits[ch] <= MAX_BITS );
    }
    assert( max_bits <= MAX_BITS );
    return max_bits;
}




void reduce_side(int targ_bits[2],FLOAT ms_ener_ratio,int mean_bits,int max_bits)
{
  int move_bits;
  FLOAT fac;


  /*  ms_ener_ratio = 0:  allocate 66/33  mid/side  fac=.33  
   *  ms_ener_ratio =.5:  allocate 50/50 mid/side   fac= 0 */
  /* 75/25 split is fac=.5 */
  /* float fac = .50*(.5-ms_ener_ratio[gr])/.5;*/
  fac = .33*(.5-ms_ener_ratio)/.5;
  if (fac<0) fac=0;
  if (fac>.5) fac=.5;
  
    /* number of bits to move from side channel to mid channel */
    /*    move_bits = fac*targ_bits[1];  */
    move_bits = fac*.5*(targ_bits[0]+targ_bits[1]);  

    if (move_bits > MAX_BITS - targ_bits[0]) {
        move_bits = MAX_BITS - targ_bits[0];
    }
    if (move_bits<0) move_bits=0;
    
    if (targ_bits[1] >= 125) {
      /* dont reduce side channel below 125 bits */
      if (targ_bits[1]-move_bits > 125) {