#include <assert.h>#include "util.h"#include "gtkanal.h"#include "tables.h"#include "reservoir.h"#include "quantize-pvt.h"/* some problems found with -O2 and above, gcc 2.95 */#if (defined(__GNUC__) && defined(__i386__))#undef TAKEHIRO_IEEE754_HACK#endifconst int slen1_tab[16] = { 0, 0, 0, 0, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4 };const int slen2_tab[16] = { 0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3, 1, 2, 3, 2, 3 };/*  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]*/unsigned 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 */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}  }};FLOAT8 pow20[Q_MAX];FLOAT8 ipow20[Q_MAX];FLOAT8 pow43[PRECALC_SIZE];/* initialized in first call to iteration_init */#ifndef TAKEHIRO_IEEE754_HACKFLOAT8 adj43[PRECALC_SIZE];#endifFLOAT8 adj43asm[PRECALC_SIZE];/************************************************************************//*  initialization for iteration_loop *//************************************************************************/voiditeration_init( lame_global_flags *gfp,III_side_info_t *l3_side, int l3_enc[2][2][576]){  lame_internal_flags *gfc=gfp->internal_flags;  gr_info *cod_info;  int ch, gr, i;  if ( gfc->iteration_init_init==0 ) {    gfc->iteration_init_init=1;    l3_side->main_data_begin = 0;    compute_ath(gfp,gfc->ATH_l,gfc->ATH_s);    for(i=0;i<PRECALC_SIZE;i++)        pow43[i] = pow((FLOAT8)i, 4.0/3.0);#ifndef TAKEHIRO_IEEE754_HACK    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    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);    for (i = 0; i < Q_MAX; i++) {	ipow20[i] = pow(2.0, (double)(i - 210) * -0.1875);	pow20[i] = pow(2.0, (double)(i - 210) * 0.25);    }  }    /* some intializations. */  for ( gr = 0; gr < gfc->mode_gr; gr++ ){    for ( ch = 0; ch < gfc->stereo; ch++ ){      cod_info = (gr_info *) &(l3_side->gr[gr].ch[ch]);      if (cod_info->block_type == SHORT_TYPE)        {	  cod_info->sfb_lmax = 0; /* No sb*/	  cod_info->sfb_smax = 0;        }      else	{	  /* MPEG 1 doesnt use last scalefactor band */	  cod_info->sfb_lmax = SBPSY_l;	  cod_info->sfb_smax = SBPSY_s;    /* No sb */	}    }  }  huffman_init();}/* compute the ATH for each scalefactor band cd range:  0..96dbInput:  3.3kHz signal  32767 amplitude  (3.3kHz is where ATH is smallest = -5db)longblocks:  sfb=12   en0/bw=-11db    max_en0 = 1.3dbshortblocks: sfb=5           -9db              0dbInput:  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 = 1dbAVE energy of largest wave at 3.3kHz = -11dbLet'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 needATH = ATH_formula  - 103  (db)ATH = ATH * 2.5e-10      (ener)*/FLOAT8 ATHmdct(lame_global_flags *gfp,FLOAT8 f){  lame_internal_flags *gfc=gfp->internal_flags;  FLOAT8 ath;    ath = ATHformula(f);	    /* convert to energy */  ath -= 114;    /* MDCT scaling.  From tests by macik and MUS420 code */  /* purpose of RH_QUALITY_CONTROL:   * at higher quality lower ATH masking abilities   => needs more bits   * at lower quality increase ATH masking abilities => needs less bits   * works together with adjusted masking lowering of GPSYCHO thresholds   * (Robert.Hegemann@gmx.de 2000-01-30)   */  if (gfp->VBR!=vbr_off)     {      ath -= gfc->ATH_lower;      ath = Min(gfp->VBR_q-62,ath);    }      ath = pow( 10.0, ath/10.0 );  return ath;} void compute_ath(lame_global_flags *gfp,FLOAT8 ATH_l[],FLOAT8 ATH_s[]){  lame_internal_flags *gfc=gfp->internal_flags;  int sfb,i,start,end;  FLOAT8 ATH_f;  FLOAT8 samp_freq = gfp->out_samplerate/1000.0;  for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {    start = gfc->scalefac_band.l[ sfb ];    end   = gfc->scalefac_band.l[ sfb+1 ];    ATH_l[sfb]=1e99;    for (i=start ; i < end; i++) {      FLOAT8 freq = samp_freq*i/(2*576);      assert( freq < 25 );      ATH_f = ATHmdct(gfp,freq);  /* freq in kHz */      ATH_l[sfb]=Min(ATH_l[sfb],ATH_f);    }    /*    DEBUGF("sfb=%2i freq(khz): %5.2f ..%5.2f  ATH=%6.2f %6.2f  %6.2f   \n",sfb,samp_freq*start/(2*576),samp_freq*end/(2*576),10*log10(ATH_l[sfb]),10*log10( ATHmdct(gfp,samp_freq*start/(2*576)))  ,10*log10(ATHmdct(gfp,samp_freq*end/(2*576))));    */  }  for ( sfb = 0; sfb < SBMAX_s; sfb++ ){    start = gfc->scalefac_band.s[ sfb ];    end   = gfc->scalefac_band.s[ sfb+1 ];    ATH_s[sfb]=1e99;    for (i=start ; i < end; i++) {      FLOAT8 freq = samp_freq*i/(2*192);      assert( freq < 25 );      ATH_f = ATHmdct(gfp,freq);    /* freq in kHz */      ATH_s[sfb]=Min(ATH_s[sfb],ATH_f);    }  }  /* in no ATH mode leave ATH for the last scalefactor band in    * because VBR mode needs it   */  if (gfp->noATH) {    for ( sfb = 0; sfb < SBMAX_l-1; sfb++ ) {      ATH_l[sfb]=1E-20;    }    for ( sfb = 0; sfb < SBMAX_s-1; sfb++ ) {      ATH_s[sfb]=1E-20;    }  }}/* convert from L/R <-> Mid/Side */void ms_convert(FLOAT8 xr[2][576],FLOAT8 xr_org[2][576]){  int i;  for ( i = 0; i < 576; i++ ) {    FLOAT8 l = xr_org[0][i];    FLOAT8 r = xr_org[1][i];    xr[0][i] = (l+r)*(SQRT2*0.5);    xr[1][i] = (l-r)*(SQRT2*0.5);  }}/************************************************************************ * allocate bits among 2 channels based on PE * mt 6/99 ************************************************************************/int on_pe(lame_global_flags *gfp,FLOAT8 pe[2][2],III_side_info_t *l3_side,int targ_bits[2],int mean_bits, int gr){  lame_internal_flags *gfc=gfp->internal_flags;  gr_info *cod_info;  int extra_bits,tbits,bits;  int add_bits[2];   int ch;  int max_bits;  /* maximum allowed bits for this granule */  /* allocate targ_bits for granule */  ResvMaxBits(gfp, mean_bits, &tbits, &extra_bits);  max_bits=tbits+extra_bits;  bits=0;  for (ch=0 ; ch < gfc->stereo ; ch ++) {    /******************************************************************     * allocate bits for each channel      ******************************************************************/    cod_info = &l3_side->gr[gr].ch[ch].tt;        targ_bits[ch]=Min(4095,tbits/gfc->stereo);        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] > .75*mean_bits) add_bits[ch]=mean_bits*.75;    if (add_bits[ch] < 0) add_bits[ch]=0;    if ((targ_bits[ch]+add_bits[ch]) > 4095)       add_bits[ch]=Max(0,4095-targ_bits[ch]);    bits += add_bits[ch];  }  if (bits > extra_bits)    for (ch=0 ; ch < gfc->stereo ; ch ++) {      add_bits[ch] = (extra_bits*add_bits[ch])/bits;    }  for (ch=0 ; ch < gfc->stereo ; ch ++) {    targ_bits[ch] = targ_bits[ch] + add_bits[ch];    extra_bits -= add_bits[ch];  }  return max_bits;}void reduce_side(int targ_bits[2],FLOAT8 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 + targ_bits[0]) > 4095) {      move_bits = 4095 - 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) {	/* if mid channel already has 2x more than average, dont bother */	/* mean_bits = bits per granule (for both channels) */	if (targ_bits[0] < mean_bits)	  targ_bits[0] += move_bits;	targ_bits[1] -= move_bits;      } else {	targ_bits[0] += targ_bits[1] - 125;	targ_bits[1] = 125;      }    }        move_bits=targ_bits[0]+targ_bits[1];    if (move_bits > max_bits) {      targ_bits[0]=(max_bits*targ_bits[0])/move_bits;      targ_bits[1]=(max_bits*targ_bits[1])/move_bits;    }}/***************************************************************************  *         inner_loop                                                      *  ***************************************************************************  * The code selects the best global gain for a particular set of scalefacs */ intinner_loop( lame_global_flags *gfp,FLOAT8 xrpow[576],	    int l3_enc[576], int max_bits,	    gr_info *cod_info){    int bits;    assert( max_bits >= 0 );    cod_info->global_gain--;    do    {      cod_info->global_gain++;      bits = count_bits(gfp,l3_enc, xrpow, cod_info);    }    while ( bits > max_bits );    return bits;}/*************************************************************************//*            scale_bitcount                                             *//*************************************************************************//* Also calculates the number of bits necessary to code the scalefactors. */int scale_bitcount( III_scalefac_t *scalefac, gr_info *cod_info){    int i, k, sfb, max_slen1 = 0, max_slen2 = 0, /*a, b, */ ep = 2;    /* maximum values */    static const int slen1[16] = { 1, 1, 1, 1, 8, 2, 2, 2, 4, 4, 4, 8, 8, 8,16,16 };    static const int slen2[16] = { 1, 2, 4, 8, 1, 2, 4, 8, 2, 4, 8, 2, 4, 8, 4, 8 };    /* number of bits used to encode scalefacs */    static const int slen1_value[16] = {0,	18, 36, 54, 54, 36, 54, 72, 54, 72, 90, 72, 90,108,108,126    };    static const int slen2_value[16] = {0,	10, 20, 30, 33, 21, 31, 41, 32, 42, 52, 43, 53, 63, 64, 74    };    const int *tab;    if ( cod_info->block_type == SHORT_TYPE )    {            tab = slen1_value;            /* a = 18; b = 18;  */            for ( i = 0; i < 3; i++ )            {                for ( sfb = 0; sfb < 6; sfb++ )                    if (scalefac->s[sfb][i] > max_slen1 )                        max_slen1 = scalefac->s[sfb][i];                for (sfb = 6; sfb < SBPSY_s; sfb++ )                    if ( scalefac->s[sfb][i] > max_slen2 )                        max_slen2 = scalefac->s[sfb][i];            }    }    else    { /* block_type == 1,2,or 3 */        tab = slen2_value;        /* a = 11; b = 10;   */        for ( sfb = 0; sfb < 11; sfb++ )            if ( scalefac->l[sfb] > max_slen1 )                max_slen1 = scalefac->l[sfb];	if (!cod_info->preflag) {	    for ( sfb = 11; sfb < SBPSY_l; sfb++ )		if (scalefac->l[sfb] < pretab[sfb])		    break;	    if (sfb == SBPSY_l) {		cod_info->preflag = 1;		for ( sfb = 11; sfb < SBPSY_l; sfb++ )		    scalefac->l[sfb] -= pretab[sfb];	    }	}        for ( sfb = 11; sfb < SBPSY_l; sfb++ )            if ( scalefac->l[sfb] > max_slen2 )                max_slen2 = scalefac->l[sfb];    }    /* from Takehiro TOMINAGA <tominaga@isoternet.org> 10/99     * loop over *all* posible values of scalefac_compress to find the     * one which uses the smallest number of bits.  ISO would stop     * at first valid index */    cod_info->part2_length = LARGE_BITS;    for ( k = 0; k < 16; k++ )    {        if ( (max_slen1 < slen1[k]) && (max_slen2 < slen2[k]) &&             ((int)cod_info->part2_length > tab[k])) {	  cod_info->part2_length=tab[k];	  cod_info->scalefac_compress=k;	  ep=0;  /* we found a suitable scalefac_compress */	}    }    return ep;}/*  table of largest scalefactor values for MPEG2*/static const unsigned int max_range_sfac_tab[6][4] ={ { 15, 15, 7,  7}, { 15, 15, 7,  0}, { 7,  3,  0,  0}, { 15, 31, 31, 0}, { 7,  7,  7,  0}, { 3,  3,  0,  0}};/*************************************************************************//*            scale_bitcount_lsf                                         *//*************************************************************************//* Also counts the number of bits to encode the scalefacs but for MPEG 2 */ /* Lower sampling frequencies  (24, 22.05 and 16 kHz.)                   */ /*  This is reverse-engineered from section 2.4.3.2 of the MPEG2 IS,     *//* "Audio Decoding Layer III"                                            */int scale_bitcount_lsf(III_scalefac_t *scalefac, gr_info *cod_info){    int table_number, row_in_table, partition, nr_sfb, window, over;    int i, sfb, max_sfac[ 4 ];    unsigned int *partition_table;    /*      Set partition table. Note that should try to use table one,      but do not yet...    */    if ( cod_info->preflag )