src/lib/dl/ext/lame/quantize_pvt.c

   1 /*
   2  *      quantize_pvt source file
   3  *
   4  *      Copyright (c) 1999-2002 Takehiro Tominaga
   5  *      Copyright (c) 2000-2011 Robert Hegemann
   6  *      Copyright (c) 2001 Naoki Shibata
   7  *      Copyright (c) 2002-2005 Gabriel Bouvigne
   8  *
   9  * This library is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Library General Public
  11  * License as published by the Free Software Foundation; either
  12  * version 2 of the License, or (at your option) any later version.
  13  *
  14  * This library is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Library General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Library General Public
  20  * License along with this library; if not, write to the
  21  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  22  * Boston, MA 02111-1307, USA.
  23  */
  24
  25 /* $Id: quantize_pvt.c,v 1.169 2011/05/24 20:45:55 robert Exp $ */
  26 #ifdef HAVE_CONFIG_H
  27 # include <config.h>
  28 #endif
  29
  30
  31 #include "lame.h"
  32 #include "machine.h"
  33 #include "encoder.h"
  34 #include "util.h"
  35 #include "quantize_pvt.h"
  36 #include "reservoir.h"
  37 #include "lame-analysis.h"
  38 #include <float.h>
  39
  40
  41 #define NSATHSCALE 100  /* Assuming dynamic range=96dB, this value should be 92 */
  42
  43 /*
  44   The following table is used to implement the scalefactor
  45   partitioning for MPEG2 as described in section
  46   2.4.3.2 of the IS. The indexing corresponds to the
  47   way the tables are presented in the IS:
  48
  49   [table_number][row_in_table][column of nr_of_sfb]
  50 */
  51 const int nr_of_sfb_block[6][3][4] = {
  52     {
  53      {6, 5, 5, 5},
  54      {9, 9, 9, 9},
  55      {6, 9, 9, 9}
  56      },
  57     {
  58      {6, 5, 7, 3},
  59      {9, 9, 12, 6},
  60      {6, 9, 12, 6}
  61      },
  62     {
  63      {11, 10, 0, 0},
  64      {18, 18, 0, 0},
  65      {15, 18, 0, 0}
  66      },
  67     {
  68      {7, 7, 7, 0},
  69      {12, 12, 12, 0},
  70      {6, 15, 12, 0}
  71      },
  72     {
  73      {6, 6, 6, 3},
  74      {12, 9, 9, 6},
  75      {6, 12, 9, 6}
  76      },
  77     {
  78      {8, 8, 5, 0},
  79      {15, 12, 9, 0},
  80      {6, 18, 9, 0}
  81      }
  82 };
  83
  84
  85 /* Table B.6: layer3 preemphasis */
  86 const int pretab[SBMAX_l] = {
  87     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  88     1, 1, 1, 1, 2, 2, 3, 3, 3, 2, 0
  89 };
  90
  91 /*
  92   Here are MPEG1 Table B.8 and MPEG2 Table B.1
  93   -- Layer III scalefactor bands.
  94   Index into this using a method such as:
  95     idx  = fr_ps->header->sampling_frequency
  96            + (fr_ps->header->version * 3)
  97 */
  98
  99
 100 const scalefac_struct sfBandIndex[9] = {
 101     {                   /* Table B.2.b: 22.05 kHz */
 102      {0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 116, 140, 168, 200, 238, 284, 336, 396, 464,
 103       522, 576},
 104      {0, 4, 8, 12, 18, 24, 32, 42, 56, 74, 100, 132, 174, 192}
 105      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 106      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 107      },
 108     {                   /* Table B.2.c: 24 kHz */ /* docs: 332. mpg123(broken): 330 */
 109      {0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 114, 136, 162, 194, 232, 278, 332, 394, 464,
 110       540, 576},
 111      {0, 4, 8, 12, 18, 26, 36, 48, 62, 80, 104, 136, 180, 192}
 112      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 113      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 114      },
 115     {                   /* Table B.2.a: 16 kHz */
 116      {0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 116, 140, 168, 200, 238, 284, 336, 396, 464,
 117       522, 576},
 118      {0, 4, 8, 12, 18, 26, 36, 48, 62, 80, 104, 134, 174, 192}
 119      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 120      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 121      },
 122     {                   /* Table B.8.b: 44.1 kHz */
 123      {0, 4, 8, 12, 16, 20, 24, 30, 36, 44, 52, 62, 74, 90, 110, 134, 162, 196, 238, 288, 342, 418,
 124       576},
 125      {0, 4, 8, 12, 16, 22, 30, 40, 52, 66, 84, 106, 136, 192}
 126      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 127      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 128      },
 129     {                   /* Table B.8.c: 48 kHz */
 130      {0, 4, 8, 12, 16, 20, 24, 30, 36, 42, 50, 60, 72, 88, 106, 128, 156, 190, 230, 276, 330, 384,
 131       576},
 132      {0, 4, 8, 12, 16, 22, 28, 38, 50, 64, 80, 100, 126, 192}
 133      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 134      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 135      },
 136     {                   /* Table B.8.a: 32 kHz */
 137      {0, 4, 8, 12, 16, 20, 24, 30, 36, 44, 54, 66, 82, 102, 126, 156, 194, 240, 296, 364, 448, 550,
 138       576},
 139      {0, 4, 8, 12, 16, 22, 30, 42, 58, 78, 104, 138, 180, 192}
 140      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 141      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 142      },
 143     {                   /* MPEG-2.5 11.025 kHz */
 144      {0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 116, 140, 168, 200, 238, 284, 336, 396, 464,
 145       522, 576},
 146      {0 / 3, 12 / 3, 24 / 3, 36 / 3, 54 / 3, 78 / 3, 108 / 3, 144 / 3, 186 / 3, 240 / 3, 312 / 3,
 147       402 / 3, 522 / 3, 576 / 3}
 148      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 149      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 150      },
 151     {                   /* MPEG-2.5 12 kHz */
 152      {0, 6, 12, 18, 24, 30, 36, 44, 54, 66, 80, 96, 116, 140, 168, 200, 238, 284, 336, 396, 464,
 153       522, 576},
 154      {0 / 3, 12 / 3, 24 / 3, 36 / 3, 54 / 3, 78 / 3, 108 / 3, 144 / 3, 186 / 3, 240 / 3, 312 / 3,
 155       402 / 3, 522 / 3, 576 / 3}
 156      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 157      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 158      },
 159     {                   /* MPEG-2.5 8 kHz */
 160      {0, 12, 24, 36, 48, 60, 72, 88, 108, 132, 160, 192, 232, 280, 336, 400, 476, 566, 568, 570,
 161       572, 574, 576},
 162      {0 / 3, 24 / 3, 48 / 3, 72 / 3, 108 / 3, 156 / 3, 216 / 3, 288 / 3, 372 / 3, 480 / 3, 486 / 3,
 163       492 / 3, 498 / 3, 576 / 3}
 164      , {0, 0, 0, 0, 0, 0, 0} /*  sfb21 pseudo sub bands */
 165      , {0, 0, 0, 0, 0, 0, 0} /*  sfb12 pseudo sub bands */
 166      }
 167 };
 168
 169
 170
 171 FLOAT   pow20[Q_MAX + Q_MAX2 + 1];
 172 FLOAT   ipow20[Q_MAX];
 173 FLOAT   pow43[PRECALC_SIZE];
 174 /* initialized in first call to iteration_init */
 175 #ifdef TAKEHIRO_IEEE754_HACK
 176 FLOAT   adj43asm[PRECALC_SIZE];
 177 #else
 178 FLOAT   adj43[PRECALC_SIZE];
 179 #endif
 180
 181 /*
 182 compute the ATH for each scalefactor band
 183 cd range:  0..96db
 184
 185 Input:  3.3kHz signal  32767 amplitude  (3.3kHz is where ATH is smallest = -5db)
 186 longblocks:  sfb=12   en0/bw=-11db    max_en0 = 1.3db
 187 shortblocks: sfb=5           -9db              0db
 188
 189 Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
 190 longblocks:  amp=1      sfb=12   en0/bw=-103 db      max_en0 = -92db
 191             amp=32767   sfb=12           -12 db                 -1.4db
 192
 193 Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
 194 shortblocks: amp=1      sfb=5   en0/bw= -99                    -86
 195             amp=32767   sfb=5           -9  db                  4db
 196
 197
 198 MAX energy of largest wave at 3.3kHz = 1db
 199 AVE energy of largest wave at 3.3kHz = -11db
 200 Let's take AVE:  -11db = maximum signal in sfb=12.
 201 Dynamic range of CD: 96db.  Therefor energy of smallest audible wave
 202 in sfb=12  = -11  - 96 = -107db = ATH at 3.3kHz.
 203
 204 ATH formula for this wave: -5db.  To adjust to LAME scaling, we need
 205 ATH = ATH_formula  - 103  (db)
 206 ATH = ATH * 2.5e-10      (ener)
 207
 208 */
 209
 210 static  FLOAT
 211 ATHmdct(SessionConfig_t const *cfg, FLOAT f)
 212 {
 213     FLOAT   ath;
 214
 215     ath = ATHformula(cfg, f);
 216
 217     if (cfg->ATHfixpoint > 0) {
 218         ath -= cfg->ATHfixpoint;
 219     }
 220     else {
 221         ath -= NSATHSCALE;
 222     }
 223     ath += cfg->ATH_offset_db;
 224
 225     /* modify the MDCT scaling for the ATH and convert to energy */
 226     ath = powf(10.0f, ath * 0.1f);
 227     return ath;
 228 }
 229
 230 static void
 231 compute_ath(lame_internal_flags const* gfc)
 232 {
 233     SessionConfig_t const *const cfg = &gfc->cfg;
 234     FLOAT  *const ATH_l = gfc->ATH->l;
 235     FLOAT  *const ATH_psfb21 = gfc->ATH->psfb21;
 236     FLOAT  *const ATH_s = gfc->ATH->s;
 237     FLOAT  *const ATH_psfb12 = gfc->ATH->psfb12;
 238     int     sfb, i, start, end;
 239     FLOAT   ATH_f;
 240     FLOAT const samp_freq = cfg->samplerate_out;
 241
 242     for (sfb = 0; sfb < SBMAX_l; sfb++) {
 243         start = gfc->scalefac_band.l[sfb];
 244         end = gfc->scalefac_band.l[sfb + 1];
 245         ATH_l[sfb] = FLOAT_MAX;
 246         for (i = start; i < end; i++) {
 247             FLOAT const freq = i * samp_freq / (2 * 576);
 248             ATH_f = ATHmdct(cfg, freq); /* freq in kHz */
 249             ATH_l[sfb] = Min(ATH_l[sfb], ATH_f);
 250         }
 251     }
 252
 253     for (sfb = 0; sfb < PSFB21; sfb++) {
 254         start = gfc->scalefac_band.psfb21[sfb];
 255         end = gfc->scalefac_band.psfb21[sfb + 1];
 256         ATH_psfb21[sfb] = FLOAT_MAX;
 257         for (i = start; i < end; i++) {
 258             FLOAT const freq = i * samp_freq / (2 * 576);
 259             ATH_f = ATHmdct(cfg, freq); /* freq in kHz */
 260             ATH_psfb21[sfb] = Min(ATH_psfb21[sfb], ATH_f);
 261         }
 262     }
 263
 264     for (sfb = 0; sfb < SBMAX_s; sfb++) {
 265         start = gfc->scalefac_band.s[sfb];
 266         end = gfc->scalefac_band.s[sfb + 1];
 267         ATH_s[sfb] = FLOAT_MAX;
 268         for (i = start; i < end; i++) {
 269             FLOAT const freq = i * samp_freq / (2 * 192);
 270             ATH_f = ATHmdct(cfg, freq); /* freq in kHz */
 271             ATH_s[sfb] = Min(ATH_s[sfb], ATH_f);
 272         }
 273         ATH_s[sfb] *= (gfc->scalefac_band.s[sfb + 1] - gfc->scalefac_band.s[sfb]);
 274     }
 275
 276     for (sfb = 0; sfb < PSFB12; sfb++) {
 277         start = gfc->scalefac_band.psfb12[sfb];
 278         end = gfc->scalefac_band.psfb12[sfb + 1];
 279         ATH_psfb12[sfb] = FLOAT_MAX;
 280         for (i = start; i < end; i++) {
 281             FLOAT const freq = i * samp_freq / (2 * 192);
 282             ATH_f = ATHmdct(cfg, freq); /* freq in kHz */
 283             ATH_psfb12[sfb] = Min(ATH_psfb12[sfb], ATH_f);
 284         }
 285         /*not sure about the following */
 286         ATH_psfb12[sfb] *= (gfc->scalefac_band.s[13] - gfc->scalefac_band.s[12]);
 287     }
 288
 289
 290     /*  no-ATH mode:
 291      *  reduce ATH to -200 dB
 292      */
 293
 294     if (cfg->noATH) {
 295         for (sfb = 0; sfb < SBMAX_l; sfb++) {
 296             ATH_l[sfb] = 1E-20;
 297         }
 298         for (sfb = 0; sfb < PSFB21; sfb++) {
 299             ATH_psfb21[sfb] = 1E-20;
 300         }
 301         for (sfb = 0; sfb < SBMAX_s; sfb++) {
 302             ATH_s[sfb] = 1E-20;
 303         }
 304         for (sfb = 0; sfb < PSFB12; sfb++) {
 305             ATH_psfb12[sfb] = 1E-20;
 306         }
 307     }
 308
 309     /*  work in progress, don't rely on it too much
 310      */
 311     gfc->ATH->floor = 10. * log10(ATHmdct(cfg, -1.));
 312
 313     /*
 314        {   FLOAT g=10000, t=1e30, x;
 315        for ( f = 100; f < 10000; f++ ) {
 316        x = ATHmdct( cfg, f );
 317        if ( t > x ) t = x, g = f;
 318        }
 319        printf("min=%g\n", g);
 320        } */
 321 }
 322
 323
 324 static float const payload_long[2][4] =
 325 { {-0.000f, -0.000f, -0.000f, +0.000f}
 326 , {-0.500f, -0.250f, -0.025f, +0.500f}
 327 };
 328 static float const payload_short[2][4] =
 329 { {-0.000f, -0.000f, -0.000f, +0.000f}
 330 , {-2.000f, -1.000f, -0.050f, +0.500f}
 331 };
 332
 333 /************************************************************************/
 334 /*  initialization for iteration_loop */
 335 /************************************************************************/
 336 void
 337 iteration_init(lame_internal_flags * gfc)
 338 {
 339     SessionConfig_t const *const cfg = &gfc->cfg;
 340     III_side_info_t *const l3_side = &gfc->l3_side;
 341     FLOAT   adjust, db;
 342     int     i, sel;
 343
 344     if (gfc->iteration_init_init == 0) {
 345         gfc->iteration_init_init = 1;
 346
 347         l3_side->main_data_begin = 0;
 348         compute_ath(gfc);
 349
 350         pow43[0] = 0.0;
 351         for (i = 1; i < PRECALC_SIZE; i++)
 352             pow43[i] = pow((FLOAT) i, 4.0 / 3.0);
 353
 354 #ifdef TAKEHIRO_IEEE754_HACK
 355         adj43asm[0] = 0.0;
 356         for (i = 1; i < PRECALC_SIZE; i++)
 357             adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]), 0.75);
 358 #else
 359         for (i = 0; i < PRECALC_SIZE - 1; i++)
 360             adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
 361         adj43[i] = 0.5;
 362 #endif
 363         for (i = 0; i < Q_MAX; i++)
 364             ipow20[i] = pow(2.0, (double) (i - 210) * -0.1875);
 365         for (i = 0; i <= Q_MAX + Q_MAX2; i++)
 366             pow20[i] = pow(2.0, (double) (i - 210 - Q_MAX2) * 0.25);
 367
 368         huffman_init(gfc);
 369         init_xrpow_core_init(gfc);
 370
 371         sel = 1;/* RH: all modes like vbr-new (cfg->vbr == vbr_mt || cfg->vbr == vbr_mtrh) ? 1 : 0;*/
 372
 373         /* long */
 374         db = cfg->adjust_bass_db + payload_long[sel][0];
 375         adjust = powf(10.f, db * 0.1f);
 376         for (i = 0; i <= 6; ++i) {
 377             gfc->sv_qnt.longfact[i] = adjust;
 378         }
 379         db = cfg->adjust_alto_db + payload_long[sel][1];
 380         adjust = powf(10.f, db * 0.1f);
 381         for (; i <= 13; ++i) {
 382             gfc->sv_qnt.longfact[i] = adjust;
 383         }
 384         db = cfg->adjust_treble_db + payload_long[sel][2];
 385         adjust = powf(10.f, db * 0.1f);
 386         for (; i <= 20; ++i) {
 387             gfc->sv_qnt.longfact[i] = adjust;
 388         }
 389         db = cfg->adjust_sfb21_db + payload_long[sel][3];
 390         adjust = powf(10.f, db * 0.1f);
 391         for (; i < SBMAX_l; ++i) {
 392             gfc->sv_qnt.longfact[i] = adjust;
 393         }
 394
 395         /* short */
 396         db = cfg->adjust_bass_db + payload_short[sel][0];
 397         adjust = powf(10.f, db * 0.1f);
 398         for (i = 0; i <= 2; ++i) {
 399             gfc->sv_qnt.shortfact[i] = adjust;
 400         }
 401         db = cfg->adjust_alto_db + payload_short[sel][1];
 402         adjust = powf(10.f, db * 0.1f);
 403         for (; i <= 6; ++i) {
 404             gfc->sv_qnt.shortfact[i] = adjust;
 405         }
 406         db = cfg->adjust_treble_db + payload_short[sel][2];
 407         adjust = powf(10.f, db * 0.1f);
 408         for (; i <= 11; ++i) {
 409             gfc->sv_qnt.shortfact[i] = adjust;
 410         }
 411         db = cfg->adjust_sfb21_db + payload_short[sel][3];
 412         adjust = powf(10.f, db * 0.1f);
 413         for (; i < SBMAX_s; ++i) {
 414             gfc->sv_qnt.shortfact[i] = adjust;
 415         }
 416     }
 417 }
 418
 419
 420
 421
 422
 423 /************************************************************************
 424  * allocate bits among 2 channels based on PE
 425  * mt 6/99
 426  * bugfixes rh 8/01: often allocated more than the allowed 4095 bits
 427  ************************************************************************/
 428 int
 429 on_pe(lame_internal_flags * gfc, const FLOAT pe[][2], int targ_bits[2], int mean_bits, int gr, int cbr)
 430 {
 431     SessionConfig_t const *const cfg = &gfc->cfg;
 432     int     extra_bits = 0, tbits, bits;
 433     int     add_bits[2] = {0, 0};
 434     int     max_bits;        /* maximum allowed bits for this granule */
 435     int     ch;
 436
 437     /* allocate targ_bits for granule */
 438     ResvMaxBits(gfc, mean_bits, &tbits, &extra_bits, cbr);
 439     max_bits = tbits + extra_bits;
 440     if (max_bits > MAX_BITS_PER_GRANULE) /* hard limit per granule */
 441         max_bits = MAX_BITS_PER_GRANULE;
 442
 443     for (bits = 0, ch = 0; ch < cfg->channels_out; ++ch) {
 444         /******************************************************************
 445          * allocate bits for each channel
 446          ******************************************************************/
 447         targ_bits[ch] = Min(MAX_BITS_PER_CHANNEL, tbits / cfg->channels_out);
 448
 449         add_bits[ch] = targ_bits[ch] * pe[gr][ch] / 700.0 - targ_bits[ch];
 450
 451         /* at most increase bits by 1.5*average */
 452         if (add_bits[ch] > mean_bits * 3 / 4)
 453             add_bits[ch] = mean_bits * 3 / 4;
 454         if (add_bits[ch] < 0)
 455             add_bits[ch] = 0;
 456
 457         if (add_bits[ch] + targ_bits[ch] > MAX_BITS_PER_CHANNEL)
 458             add_bits[ch] = Max(0, MAX_BITS_PER_CHANNEL - targ_bits[ch]);
 459
 460         bits += add_bits[ch];
 461     }
 462     if (bits > extra_bits && bits > 0) {
 463         for (ch = 0; ch < cfg->channels_out; ++ch) {
 464             add_bits[ch] = extra_bits * add_bits[ch] / bits;
 465         }
 466     }
 467
 468     for (ch = 0; ch < cfg->channels_out; ++ch) {
 469         targ_bits[ch] += add_bits[ch];
 470         extra_bits -= add_bits[ch];
 471     }
 472
 473     for (bits = 0, ch = 0; ch < cfg->channels_out; ++ch) {
 474         bits += targ_bits[ch];
 475     }
 476     if (bits > MAX_BITS_PER_GRANULE) {
 477         int     sum = 0;
 478         for (ch = 0; ch < cfg->channels_out; ++ch) {
 479             targ_bits[ch] *= MAX_BITS_PER_GRANULE;
 480             targ_bits[ch] /= bits;
 481             sum += targ_bits[ch];
 482         }
 483         assert(sum <= MAX_BITS_PER_GRANULE);
 484     }
 485
 486     return max_bits;
 487 }
 488
 489
 490
 491
 492 void
 493 reduce_side(int targ_bits[2], FLOAT ms_ener_ratio, int mean_bits, int max_bits)
 494 {
 495     int     move_bits;
 496     FLOAT   fac;
 497
 498     assert(max_bits <= MAX_BITS_PER_GRANULE);
 499     assert(targ_bits[0] + targ_bits[1] <= MAX_BITS_PER_GRANULE);
 500
 501     /*  ms_ener_ratio = 0:  allocate 66/33  mid/side  fac=.33
 502      *  ms_ener_ratio =.5:  allocate 50/50 mid/side   fac= 0 */
 503     /* 75/25 split is fac=.5 */
 504     /* float fac = .50*(.5-ms_ener_ratio[gr])/.5; */
 505     fac = .33 * (.5 - ms_ener_ratio) / .5;
 506     if (fac < 0)
 507         fac = 0;
 508     if (fac > .5)
 509         fac = .5;
 510
 511     /* number of bits to move from side channel to mid channel */
 512     /*    move_bits = fac*targ_bits[1];  */
 513     move_bits = fac * .5 * (targ_bits[0] + targ_bits[1]);
 514
 515     if (move_bits > MAX_BITS_PER_CHANNEL - targ_bits[0]) {
 516         move_bits = MAX_BITS_PER_CHANNEL - targ_bits[0];
 517     }
 518     if (move_bits < 0)
 519         move_bits = 0;
 520
 521     if (targ_bits[1] >= 125) {
 522         /* dont reduce side channel below 125 bits */
 523         if (targ_bits[1] - move_bits > 125) {
 524
 525             /* if mid channel already has 2x more than average, dont bother */
 526             /* mean_bits = bits per granule (for both channels) */
 527             if (targ_bits[0] < mean_bits)
 528                 targ_bits[0] += move_bits;
 529             targ_bits[1] -= move_bits;
 530         }
 531         else {
 532             targ_bits[0] += targ_bits[1] - 125;
 533             targ_bits[1] = 125;
 534         }
 535     }
 536
 537     move_bits = targ_bits[0] + targ_bits[1];
 538     if (move_bits > max_bits) {
 539         targ_bits[0] = (max_bits * targ_bits[0]) / move_bits;
 540         targ_bits[1] = (max_bits * targ_bits[1]) / move_bits;
 541     }
 542     assert(targ_bits[0] <= MAX_BITS_PER_CHANNEL);
 543     assert(targ_bits[1] <= MAX_BITS_PER_CHANNEL);
 544     assert(targ_bits[0] + targ_bits[1] <= MAX_BITS_PER_GRANULE);
 545 }
 546
 547
 548 /**
 549  *  Robert Hegemann 2001-04-27:
 550  *  this adjusts the ATH, keeping the original noise floor
 551  *  affects the higher frequencies more than the lower ones
 552  */
 553
 554 FLOAT
 555 athAdjust(FLOAT a, FLOAT x, FLOAT athFloor, float ATHfixpoint)
 556 {
 557     /*  work in progress
 558      */
 559     FLOAT const o = 90.30873362f;
 560     FLOAT const p = (ATHfixpoint < 1.f) ? 94.82444863f : ATHfixpoint;
 561     FLOAT   u = FAST_LOG10_X(x, 10.0f);
 562     FLOAT const v = a * a;
 563     FLOAT   w = 0.0f;
 564     u -= athFloor;      /* undo scaling */
 565     if (v > 1E-20f)
 566         w = 1.f + FAST_LOG10_X(v, 10.0f / o);
 567     if (w < 0)
 568         w = 0.f;
 569     u *= w;
 570     u += athFloor + o - p; /* redo scaling */
 571
 572     return powf(10.f, 0.1f * u);
 573 }
 574
 575
 576
 577 /*************************************************************************/
 578 /*            calc_xmin                                                  */
 579 /*************************************************************************/
 580
 581 /*
 582   Calculate the allowed distortion for each scalefactor band,
 583   as determined by the psychoacoustic model.
 584   xmin(sb) = ratio(sb) * en(sb) / bw(sb)
 585
 586   returns number of sfb's with energy > ATH
 587 */
 588
 589 int
 590 calc_xmin(lame_internal_flags const *gfc,
 591           III_psy_ratio const *const ratio, gr_info * const cod_info, FLOAT * pxmin)
 592 {
 593     SessionConfig_t const *const cfg = &gfc->cfg;
 594     int     sfb, gsfb, j = 0, ath_over = 0, k;
 595     ATH_t const *const ATH = gfc->ATH;
 596     const FLOAT *const xr = cod_info->xr;
 597     int     max_nonzero;
 598
 599     for (gsfb = 0; gsfb < cod_info->psy_lmax; gsfb++) {
 600         FLOAT   en0, xmin;
 601         FLOAT   rh1, rh2, rh3;
 602         int     width, l;
 603
 604         xmin = athAdjust(ATH->adjust_factor, ATH->l[gsfb], ATH->floor, cfg->ATHfixpoint);
 605         xmin *= gfc->sv_qnt.longfact[gsfb];
 606
 607         width = cod_info->width[gsfb];
 608         rh1 = xmin / width;
 609 #ifdef DBL_EPSILON
 610         rh2 = DBL_EPSILON;
 611 #else
 612         rh2 = 2.2204460492503131e-016;
 613 #endif
 614         en0 = 0.0;
 615         for (l = 0; l < width; ++l) {
 616             FLOAT const xa = xr[j++];
 617             FLOAT const x2 = xa * xa;
 618             en0 += x2;
 619             rh2 += (x2 < rh1) ? x2 : rh1;
 620         }
 621         if (en0 > xmin)
 622             ath_over++;
 623
 624         if (en0 < xmin) {
 625             rh3 = en0;
 626         }
 627         else if (rh2 < xmin) {
 628             rh3 = xmin;
 629         }
 630         else {
 631             rh3 = rh2;
 632         }
 633         xmin = rh3;
 634         {
 635             FLOAT const e = ratio->en.l[gsfb];
 636             if (e > 1e-12f) {
 637                 FLOAT   x;
 638                 x = en0 * ratio->thm.l[gsfb] / e;
 639                 x *= gfc->sv_qnt.longfact[gsfb];
 640                 if (xmin < x)
 641                     xmin = x;
 642             }
 643         }
 644         *pxmin++ = xmin;
 645     }                   /* end of long block loop */
 646
 647
 648
 649
 650     /*use this function to determine the highest non-zero coeff */
 651     max_nonzero = 575;
 652     if (cod_info->block_type != SHORT_TYPE) { /* NORM, START or STOP type, but not SHORT */
 653         k = 576;
 654         while (k-- && fabs(xr[k]) < 1e-12f) {
 655             max_nonzero = k;
 656         }
 657     }
 658     cod_info->max_nonzero_coeff = max_nonzero;
 659
 660
 661
 662     for (sfb = cod_info->sfb_smin; gsfb < cod_info->psymax; sfb++, gsfb += 3) {
 663         int     width, b, l;
 664         FLOAT   tmpATH;
 665
 666         tmpATH = athAdjust(ATH->adjust_factor, ATH->s[sfb], ATH->floor, cfg->ATHfixpoint);
 667         tmpATH *= gfc->sv_qnt.shortfact[sfb];
 668
 669         width = cod_info->width[gsfb];
 670         for (b = 0; b < 3; b++) {
 671             FLOAT   en0 = 0.0, xmin = tmpATH;
 672             FLOAT   rh1, rh2, rh3;
 673
 674             rh1 = tmpATH / width;
 675 #ifdef DBL_EPSILON
 676             rh2 = DBL_EPSILON;
 677 #else
 678             rh2 = 2.2204460492503131e-016;
 679 #endif
 680             for (l = 0; l < width; ++l) {
 681                 FLOAT const xa = xr[j++];
 682                 FLOAT const x2 = xa * xa;
 683                 en0 += x2;
 684                 rh2 += (x2 < rh1) ? x2 : rh1;
 685             }
 686             if (en0 > tmpATH)
 687                 ath_over++;
 688
 689             if (en0 < tmpATH) {
 690                 rh3 = en0;
 691             }
 692             else if (rh2 < tmpATH) {
 693                 rh3 = tmpATH;
 694             }
 695             else {
 696                 rh3 = rh2;
 697             }
 698             xmin = rh3;
 699             {
 700                 FLOAT const e = ratio->en.s[sfb][b];
 701                 if (e > 1e-12f) {
 702                     FLOAT   x;
 703                     x = en0 * ratio->thm.s[sfb][b] / e;
 704                     x *= gfc->sv_qnt.shortfact[sfb];
 705                     if (xmin < x)
 706                         xmin = x;
 707                 }
 708             }
 709             *pxmin++ = xmin;
 710         }               /* b */
 711         if (cfg->use_temporal_masking_effect) {
 712             if (pxmin[-3] > pxmin[-3 + 1])
 713                 pxmin[-3 + 1] += (pxmin[-3] - pxmin[-3 + 1]) * gfc->cd_psy->decay;
 714             if (pxmin[-3 + 1] > pxmin[-3 + 2])
 715                 pxmin[-3 + 2] += (pxmin[-3 + 1] - pxmin[-3 + 2]) * gfc->cd_psy->decay;
 716         }
 717     }                   /* end of short block sfb loop */
 718
 719     return ath_over;
 720 }
 721
 722
 723 static  FLOAT
 724 calc_noise_core_c(const gr_info * const cod_info, int *startline, int l, FLOAT step)
 725 {
 726     FLOAT   noise = 0;
 727     int     j = *startline;
 728     const int *const ix = cod_info->l3_enc;
 729
 730     if (j > cod_info->count1) {
 731         while (l--) {
 732             FLOAT   temp;
 733             temp = cod_info->xr[j];
 734             j++;
 735             noise += temp * temp;
 736             temp = cod_info->xr[j];
 737             j++;
 738             noise += temp * temp;
 739         }
 740     }
 741     else if (j > cod_info->big_values) {
 742         FLOAT   ix01[2];
 743         ix01[0] = 0;
 744         ix01[1] = step;
 745         while (l--) {
 746             FLOAT   temp;
 747             temp = fabs(cod_info->xr[j]) - ix01[ix[j]];
 748             j++;
 749             noise += temp * temp;
 750             temp = fabs(cod_info->xr[j]) - ix01[ix[j]];
 751             j++;
 752             noise += temp * temp;
 753         }
 754     }
 755     else {
 756         while (l--) {
 757             FLOAT   temp;
 758             temp = fabs(cod_info->xr[j]) - pow43[ix[j]] * step;
 759             j++;
 760             noise += temp * temp;
 761             temp = fabs(cod_info->xr[j]) - pow43[ix[j]] * step;
 762             j++;
 763             noise += temp * temp;
 764         }
 765     }
 766
 767     *startline = j;
 768     return noise;
 769 }
 770
 771
 772 /*************************************************************************/
 773 /*            calc_noise                                                 */
 774 /*************************************************************************/
 775
 776 /* -oo dB  =>  -1.00 */
 777 /* - 6 dB  =>  -0.97 */
 778 /* - 3 dB  =>  -0.80 */
 779 /* - 2 dB  =>  -0.64 */
 780 /* - 1 dB  =>  -0.38 */
 781 /*   0 dB  =>   0.00 */
 782 /* + 1 dB  =>  +0.49 */
 783 /* + 2 dB  =>  +1.06 */
 784 /* + 3 dB  =>  +1.68 */
 785 /* + 6 dB  =>  +3.69 */
 786 /* +10 dB  =>  +6.45 */
 787
 788 int
 789 calc_noise(gr_info const *const cod_info,
 790            FLOAT const *l3_xmin,
 791            FLOAT * distort, calc_noise_result * const res, calc_noise_data * prev_noise)
 792 {
 793     int     sfb, l, over = 0;
 794     FLOAT   over_noise_db = 0;
 795     FLOAT   tot_noise_db = 0; /*    0 dB relative to masking */
 796     FLOAT   max_noise = -20.0; /* -200 dB relative to masking */
 797     int     j = 0;
 798     const int *scalefac = cod_info->scalefac;
 799
 800     res->over_SSD = 0;
 801
 802
 803     for (sfb = 0; sfb < cod_info->psymax; sfb++) {
 804         int const s =
 805             cod_info->global_gain - (((*scalefac++) + (cod_info->preflag ? pretab[sfb] : 0))
 806                                      << (cod_info->scalefac_scale + 1))
 807             - cod_info->subblock_gain[cod_info->window[sfb]] * 8;
 808         FLOAT   noise = 0.0;
 809
 810         if (prev_noise && (prev_noise->step[sfb] == s)) {
 811
 812             /* use previously computed values */
 813             noise = prev_noise->noise[sfb];
 814             j += cod_info->width[sfb];
 815             *distort++ = noise / *l3_xmin++;
 816
 817             noise = prev_noise->noise_log[sfb];
 818
 819         }
 820         else {
 821             FLOAT const step = POW20(s);
 822             l = cod_info->width[sfb] >> 1;
 823
 824             if ((j + cod_info->width[sfb]) > cod_info->max_nonzero_coeff) {
 825                 int     usefullsize;
 826                 usefullsize = cod_info->max_nonzero_coeff - j + 1;
 827
 828                 if (usefullsize > 0)
 829                     l = usefullsize >> 1;
 830                 else
 831                     l = 0;
 832             }
 833
 834             noise = calc_noise_core_c(cod_info, &j, l, step);
 835
 836
 837             if (prev_noise) {
 838                 /* save noise values */
 839                 prev_noise->step[sfb] = s;
 840                 prev_noise->noise[sfb] = noise;
 841             }
 842
 843             noise = *distort++ = noise / *l3_xmin++;
 844
 845             /* multiplying here is adding in dB, but can overflow */
 846             noise = FAST_LOG10(Max(noise, 1E-20));
 847
 848             if (prev_noise) {
 849                 /* save noise values */
 850                 prev_noise->noise_log[sfb] = noise;
 851             }
 852         }
 853
 854         if (prev_noise) {
 855             /* save noise values */
 856             prev_noise->global_gain = cod_info->global_gain;;
 857         }
 858
 859
 860         /*tot_noise *= Max(noise, 1E-20); */
 861         tot_noise_db += noise;
 862
 863         if (noise > 0.0) {
 864             int     tmp;
 865
 866             tmp = Max((int) (noise * 10 + .5), 1);
 867             res->over_SSD += tmp * tmp;
 868
 869             over++;
 870             /* multiplying here is adding in dB -but can overflow */
 871             /*over_noise *= noise; */
 872             over_noise_db += noise;
 873         }
 874         max_noise = Max(max_noise, noise);
 875
 876     }
 877
 878     res->over_count = over;
 879     res->tot_noise = tot_noise_db;
 880     res->over_noise = over_noise_db;
 881     res->max_noise = max_noise;
 882
 883     return over;
 884 }
 885
 886
 887
 888
 889
 890
 891
 892
 893 /************************************************************************
 894  *
 895  *  set_pinfo()
 896  *
 897  *  updates plotting data
 898  *
 899  *  Mark Taylor 2000-??-??
 900  *
 901  *  Robert Hegemann: moved noise/distortion calc into it
 902  *
 903  ************************************************************************/
 904
 905 static void
 906 set_pinfo(lame_internal_flags const *gfc,
 907           gr_info * const cod_info, const III_psy_ratio * const ratio, const int gr, const int ch)
 908 {
 909     SessionConfig_t const *const cfg = &gfc->cfg;
 910     int     sfb, sfb2;
 911     int     j, i, l, start, end, bw;
 912     FLOAT   en0, en1;
 913     FLOAT const ifqstep = (cod_info->scalefac_scale == 0) ? .5 : 1.0;
 914     int const *const scalefac = cod_info->scalefac;
 915
 916     FLOAT   l3_xmin[SFBMAX], xfsf[SFBMAX];
 917     calc_noise_result noise;
 918
 919     (void) calc_xmin(gfc, ratio, cod_info, l3_xmin);
 920     (void) calc_noise(cod_info, l3_xmin, xfsf, &noise, 0);
 921
 922     j = 0;
 923     sfb2 = cod_info->sfb_lmax;
 924     if (cod_info->block_type != SHORT_TYPE && !cod_info->mixed_block_flag)
 925         sfb2 = 22;
 926     for (sfb = 0; sfb < sfb2; sfb++) {
 927         start = gfc->scalefac_band.l[sfb];
 928         end = gfc->scalefac_band.l[sfb + 1];
 929         bw = end - start;
 930         for (en0 = 0.0; j < end; j++)
 931             en0 += cod_info->xr[j] * cod_info->xr[j];
 932         en0 /= bw;
 933         /* convert to MDCT units */
 934         en1 = 1e15;     /* scaling so it shows up on FFT plot */
 935         gfc->pinfo->en[gr][ch][sfb] = en1 * en0;
 936         gfc->pinfo->xfsf[gr][ch][sfb] = en1 * l3_xmin[sfb] * xfsf[sfb] / bw;
 937
 938         if (ratio->en.l[sfb] > 0 && !cfg->ATHonly)
 939             en0 = en0 / ratio->en.l[sfb];
 940         else
 941             en0 = 0.0;
 942
 943         gfc->pinfo->thr[gr][ch][sfb] = en1 * Max(en0 * ratio->thm.l[sfb], gfc->ATH->l[sfb]);
 944
 945         /* there is no scalefactor bands >= SBPSY_l */
 946         gfc->pinfo->LAMEsfb[gr][ch][sfb] = 0;
 947         if (cod_info->preflag && sfb >= 11)
 948             gfc->pinfo->LAMEsfb[gr][ch][sfb] = -ifqstep * pretab[sfb];
 949
 950         if (sfb < SBPSY_l) {
 951             assert(scalefac[sfb] >= 0); /* scfsi should be decoded by caller side */
 952             gfc->pinfo->LAMEsfb[gr][ch][sfb] -= ifqstep * scalefac[sfb];
 953         }
 954     }                   /* for sfb */
 955
 956     if (cod_info->block_type == SHORT_TYPE) {
 957         sfb2 = sfb;
 958         for (sfb = cod_info->sfb_smin; sfb < SBMAX_s; sfb++) {
 959             start = gfc->scalefac_band.s[sfb];
 960             end = gfc->scalefac_band.s[sfb + 1];
 961             bw = end - start;
 962             for (i = 0; i < 3; i++) {
 963                 for (en0 = 0.0, l = start; l < end; l++) {
 964                     en0 += cod_info->xr[j] * cod_info->xr[j];
 965                     j++;
 966                 }
 967                 en0 = Max(en0 / bw, 1e-20);
 968                 /* convert to MDCT units */
 969                 en1 = 1e15; /* scaling so it shows up on FFT plot */
 970
 971                 gfc->pinfo->en_s[gr][ch][3 * sfb + i] = en1 * en0;
 972                 gfc->pinfo->xfsf_s[gr][ch][3 * sfb + i] = en1 * l3_xmin[sfb2] * xfsf[sfb2] / bw;
 973                 if (ratio->en.s[sfb][i] > 0)
 974                     en0 = en0 / ratio->en.s[sfb][i];
 975                 else
 976                     en0 = 0.0;
 977                 if (cfg->ATHonly || cfg->ATHshort)
 978                     en0 = 0;
 979
 980                 gfc->pinfo->thr_s[gr][ch][3 * sfb + i] =
 981                     en1 * Max(en0 * ratio->thm.s[sfb][i], gfc->ATH->s[sfb]);
 982
 983                 /* there is no scalefactor bands >= SBPSY_s */
 984                 gfc->pinfo->LAMEsfb_s[gr][ch][3 * sfb + i]
 985                     = -2.0 * cod_info->subblock_gain[i];
 986                 if (sfb < SBPSY_s) {
 987                     gfc->pinfo->LAMEsfb_s[gr][ch][3 * sfb + i] -= ifqstep * scalefac[sfb2];
 988                 }
 989                 sfb2++;
 990             }
 991         }
 992     }                   /* block type short */
 993     gfc->pinfo->LAMEqss[gr][ch] = cod_info->global_gain;
 994     gfc->pinfo->LAMEmainbits[gr][ch] = cod_info->part2_3_length + cod_info->part2_length;
 995     gfc->pinfo->LAMEsfbits[gr][ch] = cod_info->part2_length;
 996
 997     gfc->pinfo->over[gr][ch] = noise.over_count;
 998     gfc->pinfo->max_noise[gr][ch] = noise.max_noise * 10.0;
 999     gfc->pinfo->over_noise[gr][ch] = noise.over_noise * 10.0;
1000     gfc->pinfo->tot_noise[gr][ch] = noise.tot_noise * 10.0;
1001     gfc->pinfo->over_SSD[gr][ch] = noise.over_SSD;
1002 }
1003
1004
1005 /************************************************************************
1006  *
1007  *  set_frame_pinfo()
1008  *
1009  *  updates plotting data for a whole frame
1010  *
1011  *  Robert Hegemann 2000-10-21
1012  *
1013  ************************************************************************/
1014
1015 void
1016 set_frame_pinfo(lame_internal_flags * gfc, const III_psy_ratio ratio[2][2])
1017 {
1018     SessionConfig_t const *const cfg = &gfc->cfg;
1019     int     ch;
1020     int     gr;
1021
1022     /* for every granule and channel patch l3_enc and set info
1023      */
1024     for (gr = 0; gr < cfg->mode_gr; gr++) {
1025         for (ch = 0; ch < cfg->channels_out; ch++) {
1026             gr_info *const cod_info = &gfc->l3_side.tt[gr][ch];
1027             int     scalefac_sav[SFBMAX];
1028             memcpy(scalefac_sav, cod_info->scalefac, sizeof(scalefac_sav));
1029
1030             /* reconstruct the scalefactors in case SCFSI was used
1031              */
1032             if (gr == 1) {
1033                 int     sfb;
1034                 for (sfb = 0; sfb < cod_info->sfb_lmax; sfb++) {
1035                     if (cod_info->scalefac[sfb] < 0) /* scfsi */
1036                         cod_info->scalefac[sfb] = gfc->l3_side.tt[0][ch].scalefac[sfb];
1037                 }
1038             }
1039
1040             set_pinfo(gfc, cod_info, &ratio[gr][ch], gr, ch);
1041             memcpy(cod_info->scalefac, scalefac_sav, sizeof(scalefac_sav));
1042         }               /* for ch */
1043     }                   /* for gr */
1044 }