lib/lame/quantize_pvt.c

   1 /*
   2  *      quantize_pvt source file
   3  *
   4  *      Copyright (c) 1999 Takehiro TOMINAGA
   5  *
   6  * This library is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Library General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2 of the License, or (at your option) any later version.
  10  *
  11  * This library is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Library General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Library General Public
  17  * License along with this library; if not, write to the
  18  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19  * Boston, MA 02111-1307, USA.
  20  */
  21
  22 /* $Id: quantize_pvt.c,v 1.2 2006/02/09 16:55:31 kramm Exp $ */
  23 #include <stdlib.h>
  24 #include "config_static.h"
  25
  26 #include <assert.h>
  27 #include "util.h"
  28 #include "lame-analysis.h"
  29 #include "tables.h"
  30 #include "reservoir.h"
  31 #include "quantize_pvt.h"
  32
  33 #ifdef WITH_DMALLOC
  34 #include <dmalloc.h>
  35 #endif
  36
  37
  38 #define NSATHSCALE 100 // Assuming dynamic range=96dB, this value should be 92
  39
  40 const char  slen1_tab [16] = { 0, 0, 0, 0, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4 };
  41 const char  slen2_tab [16] = { 0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3, 1, 2, 3, 2, 3 };
  42
  43
  44 /*
  45   The following table is used to implement the scalefactor
  46   partitioning for MPEG2 as described in section
  47   2.4.3.2 of the IS. The indexing corresponds to the
  48   way the tables are presented in the IS:
  49
  50   [table_number][row_in_table][column of nr_of_sfb]
  51 */
  52 const int  nr_of_sfb_block [6] [3] [4] =
  53 {
  54   {
  55     {6, 5, 5, 5},
  56     {9, 9, 9, 9},
  57     {6, 9, 9, 9}
  58   },
  59   {
  60     {6, 5, 7, 3},
  61     {9, 9, 12, 6},
  62     {6, 9, 12, 6}
  63   },
  64   {
  65     {11, 10, 0, 0},
  66     {18, 18, 0, 0},
  67     {15,18,0,0}
  68   },
  69   {
  70     {7, 7, 7, 0},
  71     {12, 12, 12, 0},
  72     {6, 15, 12, 0}
  73   },
  74   {
  75     {6, 6, 6, 3},
  76     {12, 9, 9, 6},
  77     {6, 12, 9, 6}
  78   },
  79   {
  80     {8, 8, 5, 0},
  81     {15,12,9,0},
  82     {6,18,9,0}
  83   }
  84 };
  85
  86
  87 /* Table B.6: layer3 preemphasis */
  88 const char  pretab [SBMAX_l] =
  89 {
  90     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  91     1, 1, 1, 1, 2, 2, 3, 3, 3, 2, 0
  92 };
  93
  94 /*
  95   Here are MPEG1 Table B.8 and MPEG2 Table B.1
  96   -- Layer III scalefactor bands.
  97   Index into this using a method such as:
  98     idx  = fr_ps->header->sampling_frequency
  99            + (fr_ps->header->version * 3)
 100 */
 101
 102
 103 const scalefac_struct sfBandIndex[9] =
 104 {
 105   { /* Table B.2.b: 22.05 kHz */
 106     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
 107     {0,4,8,12,18,24,32,42,56,74,100,132,174,192}
 108   },
 109   { /* Table B.2.c: 24 kHz */                 /* docs: 332. mpg123(broken): 330 */
 110     {0,6,12,18,24,30,36,44,54,66,80,96,114,136,162,194,232,278, 332, 394,464,540,576},
 111     {0,4,8,12,18,26,36,48,62,80,104,136,180,192}
 112   },
 113   { /* Table B.2.a: 16 kHz */
 114     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
 115     {0,4,8,12,18,26,36,48,62,80,104,134,174,192}
 116   },
 117   { /* Table B.8.b: 44.1 kHz */
 118     {0,4,8,12,16,20,24,30,36,44,52,62,74,90,110,134,162,196,238,288,342,418,576},
 119     {0,4,8,12,16,22,30,40,52,66,84,106,136,192}
 120   },
 121   { /* Table B.8.c: 48 kHz */
 122     {0,4,8,12,16,20,24,30,36,42,50,60,72,88,106,128,156,190,230,276,330,384,576},
 123     {0,4,8,12,16,22,28,38,50,64,80,100,126,192}
 124   },
 125   { /* Table B.8.a: 32 kHz */
 126     {0,4,8,12,16,20,24,30,36,44,54,66,82,102,126,156,194,240,296,364,448,550,576},
 127     {0,4,8,12,16,22,30,42,58,78,104,138,180,192}
 128   },
 129   { /* MPEG-2.5 11.025 kHz */
 130     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
 131     {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}
 132   },
 133   { /* MPEG-2.5 12 kHz */
 134     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
 135     {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}
 136   },
 137   { /* MPEG-2.5 8 kHz */
 138     {0,12,24,36,48,60,72,88,108,132,160,192,232,280,336,400,476,566,568,570,572,574,576},
 139     {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}
 140   }
 141 };
 142
 143
 144
 145 FLOAT8 pow20[Q_MAX];
 146 FLOAT8 ipow20[Q_MAX];
 147 FLOAT8 iipow20[Q_MAX];
 148 FLOAT8 *iipow20_;
 149 FLOAT8 pow43[PRECALC_SIZE];
 150 /* initialized in first call to iteration_init */
 151 FLOAT8 adj43asm[PRECALC_SIZE];
 152 FLOAT8 adj43[PRECALC_SIZE];
 153
 154 /************************************************************************/
 155 /*  initialization for iteration_loop */
 156 /************************************************************************/
 157 void
 158 iteration_init( lame_global_flags *gfp)
 159 {
 160   lame_internal_flags *gfc=gfp->internal_flags;
 161   III_side_info_t * const l3_side = &gfc->l3_side;
 162   int i;
 163
 164   if ( gfc->iteration_init_init==0 ) {
 165     gfc->iteration_init_init=1;
 166
 167     l3_side->main_data_begin = 0;
 168     compute_ath(gfp,gfc->ATH->l,gfc->ATH->s);
 169
 170     pow43[0] = 0.0;
 171     for(i=1;i<PRECALC_SIZE;i++)
 172         pow43[i] = pow((FLOAT8)i, 4.0/3.0);
 173
 174     adj43asm[0] = 0.0;
 175     for (i = 1; i < PRECALC_SIZE; i++)
 176       adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]),0.75);
 177     for (i = 0; i < PRECALC_SIZE-1; i++)
 178         adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
 179     adj43[i] = 0.5;
 180     iipow20_ = &iipow20[210];
 181     for (i = 0; i < Q_MAX; i++) {
 182         iipow20[i] = pow(2.0, (double)(i - 210) * 0.1875);
 183         ipow20[i] = pow(2.0, (double)(i - 210) * -0.1875);
 184         pow20[i] = pow(2.0, (double)(i - 210) * 0.25);
 185     }
 186     huffman_init(gfc);
 187   }
 188 }
 189
 190
 191
 192
 193
 194 /*
 195 compute the ATH for each scalefactor band
 196 cd range:  0..96db
 197
 198 Input:  3.3kHz signal  32767 amplitude  (3.3kHz is where ATH is smallest = -5db)
 199 longblocks:  sfb=12   en0/bw=-11db    max_en0 = 1.3db
 200 shortblocks: sfb=5           -9db              0db
 201
 202 Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
 203 longblocks:  amp=1      sfb=12   en0/bw=-103 db      max_en0 = -92db
 204             amp=32767   sfb=12           -12 db                 -1.4db
 205
 206 Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
 207 shortblocks: amp=1      sfb=5   en0/bw= -99                    -86
 208             amp=32767   sfb=5           -9  db                  4db
 209
 210
 211 MAX energy of largest wave at 3.3kHz = 1db
 212 AVE energy of largest wave at 3.3kHz = -11db
 213 Let's take AVE:  -11db = maximum signal in sfb=12.
 214 Dynamic range of CD: 96db.  Therefor energy of smallest audible wave
 215 in sfb=12  = -11  - 96 = -107db = ATH at 3.3kHz.
 216
 217 ATH formula for this wave: -5db.  To adjust to LAME scaling, we need
 218 ATH = ATH_formula  - 103  (db)
 219 ATH = ATH * 2.5e-10      (ener)
 220
 221 */
 222
 223 FLOAT8 ATHmdct( lame_global_flags *gfp, FLOAT8 f )
 224 {
 225     lame_internal_flags *gfc = gfp->internal_flags;
 226     FLOAT8 ath;
 227
 228     ath = ATHformula( f , gfp );
 229
 230     if (gfc->nsPsy.use) {
 231         ath -= NSATHSCALE;
 232     } else {
 233         ath -= 114;
 234     }
 235
 236     /*  modify the MDCT scaling for the ATH  */
 237     ath -= gfp->ATHlower;
 238
 239     /* convert to energy */
 240     ath = pow( 10.0, ath/10.0 );
 241     return ath;
 242 }
 243
 244
 245 void compute_ath( lame_global_flags *gfp, FLOAT8 ATH_l[], FLOAT8 ATH_s[] )
 246 {
 247     lame_internal_flags *gfc = gfp->internal_flags;
 248     int sfb, i, start, end;
 249     FLOAT8 ATH_f;
 250     FLOAT8 samp_freq = gfp->out_samplerate;
 251
 252     for (sfb = 0; sfb < SBMAX_l; sfb++) {
 253         start = gfc->scalefac_band.l[ sfb ];
 254         end   = gfc->scalefac_band.l[ sfb+1 ];
 255         ATH_l[sfb]=FLOAT8_MAX;
 256         for (i = start ; i < end; i++) {
 257             FLOAT8 freq = i*samp_freq/(2*576);
 258             ATH_f = ATHmdct( gfp, freq );  /* freq in kHz */
 259             ATH_l[sfb] = Min( ATH_l[sfb], ATH_f );
 260         }
 261     }
 262
 263     for (sfb = 0; sfb < SBMAX_s; sfb++){
 264         start = gfc->scalefac_band.s[ sfb ];
 265         end   = gfc->scalefac_band.s[ sfb+1 ];
 266         ATH_s[sfb] = FLOAT8_MAX;
 267         for (i = start ; i < end; i++) {
 268             FLOAT8 freq = i*samp_freq/(2*192);
 269             ATH_f = ATHmdct( gfp, freq );    /* freq in kHz */
 270             ATH_s[sfb] = Min( ATH_s[sfb], ATH_f );
 271         }
 272     }
 273
 274     /*  no-ATH mode:
 275      *  reduce ATH to -200 dB
 276      */
 277
 278     if (gfp->noATH) {
 279         for (sfb = 0; sfb < SBMAX_l; sfb++) {
 280             ATH_l[sfb] = 1E-37;
 281         }
 282         for (sfb = 0; sfb < SBMAX_s; sfb++) {
 283             ATH_s[sfb] = 1E-37;
 284         }
 285     }
 286
 287     /*  work in progress, don't rely on it too much
 288      */
 289     gfc->ATH-> floor = 10. * log10( ATHmdct( gfp, -1. ) );
 290
 291     /*
 292     {   FLOAT8 g=10000, t=1e30, x;
 293         for ( f = 100; f < 10000; f++ ) {
 294             x = ATHmdct( gfp, f );
 295             if ( t > x ) t = x, g = f;
 296         }
 297         printf("min=%g\n", g);
 298     }*/
 299 }
 300
 301
 302
 303
 304
 305 /* convert from L/R <-> Mid/Side, src == dst allowed */
 306 void ms_convert(FLOAT8 dst[2][576], FLOAT8 src[2][576])
 307 {
 308     FLOAT8 l;
 309     FLOAT8 r;
 310     int i;
 311     for (i = 0; i < 576; ++i) {
 312         l = src[0][i];
 313         r = src[1][i];
 314         dst[0][i] = (l+r) * (FLOAT8)(SQRT2*0.5);
 315         dst[1][i] = (l-r) * (FLOAT8)(SQRT2*0.5);
 316     }
 317 }
 318
 319
 320
 321 /************************************************************************
 322  * allocate bits among 2 channels based on PE
 323  * mt 6/99
 324  * bugfixes rh 8/01: often allocated more than the allowed 4095 bits
 325  ************************************************************************/
 326 int on_pe( lame_global_flags *gfp, FLOAT8 pe[][2], III_side_info_t *l3_side,
 327            int targ_bits[2], int mean_bits, int gr )
 328 {
 329     lame_internal_flags * gfc = gfp->internal_flags;
 330     gr_info *   cod_info;
 331     int     extra_bits, tbits, bits;
 332     int     add_bits[2];
 333     int     max_bits;  /* maximum allowed bits for this granule */
 334     int     ch;
 335
 336     /* allocate targ_bits for granule */
 337     ResvMaxBits( gfp, mean_bits, &tbits, &extra_bits );
 338     max_bits = tbits + extra_bits;
 339     if (max_bits > MAX_BITS) /* hard limit per granule */
 340         max_bits = MAX_BITS;
 341
 342     for ( bits = 0, ch = 0; ch < gfc->channels_out; ++ch ) {
 343         /******************************************************************
 344          * allocate bits for each channel
 345          ******************************************************************/
 346         cod_info = &l3_side->gr[gr].ch[ch].tt;
 347
 348         targ_bits[ch] = Min( MAX_BITS, tbits/gfc->channels_out );
 349
 350         if (gfc->nsPsy.use) {
 351             add_bits[ch] = targ_bits[ch] * pe[gr][ch] / 700.0 - targ_bits[ch];
 352         }
 353         else {
 354             add_bits[ch] = (pe[gr][ch]-750) / 1.4;
 355             /* short blocks us a little extra, no matter what the pe */
 356             if ( cod_info->block_type == SHORT_TYPE ) {
 357                 if (add_bits[ch] < mean_bits/4)
 358                     add_bits[ch] = mean_bits/4;
 359             }
 360         }
 361
 362         /* at most increase bits by 1.5*average */
 363         if (add_bits[ch] > mean_bits*3/4)
 364             add_bits[ch] = mean_bits*3/4;
 365         if (add_bits[ch] < 0)
 366             add_bits[ch] = 0;
 367
 368         if (add_bits[ch]+targ_bits[ch] > MAX_BITS)
 369             add_bits[ch] = Max( 0, MAX_BITS-targ_bits[ch] );
 370
 371         bits += add_bits[ch];
 372     }
 373     if (bits > extra_bits) {
 374         for ( ch = 0; ch < gfc->channels_out; ++ch ) {
 375             add_bits[ch] = extra_bits * add_bits[ch] / bits;
 376         }
 377     }
 378
 379     for ( ch = 0; ch < gfc->channels_out; ++ch ) {
 380         targ_bits[ch] += add_bits[ch];
 381         extra_bits    -= add_bits[ch];
 382         assert( targ_bits[ch] <= MAX_BITS );
 383     }
 384     assert( max_bits <= MAX_BITS );
 385     return max_bits;
 386 }
 387
 388
 389
 390
 391 void reduce_side(int targ_bits[2],FLOAT8 ms_ener_ratio,int mean_bits,int max_bits)
 392 {
 393   int move_bits;
 394   FLOAT fac;
 395
 396
 397   /*  ms_ener_ratio = 0:  allocate 66/33  mid/side  fac=.33
 398    *  ms_ener_ratio =.5:  allocate 50/50 mid/side   fac= 0 */
 399   /* 75/25 split is fac=.5 */
 400   /* float fac = .50*(.5-ms_ener_ratio[gr])/.5;*/
 401   fac = .33*(.5-ms_ener_ratio)/.5;
 402   if (fac<0) fac=0;
 403   if (fac>.5) fac=.5;
 404
 405     /* number of bits to move from side channel to mid channel */
 406     /*    move_bits = fac*targ_bits[1];  */
 407     move_bits = fac*.5*(targ_bits[0]+targ_bits[1]);
 408
 409     if (move_bits > MAX_BITS - targ_bits[0]) {
 410         move_bits = MAX_BITS - targ_bits[0];
 411     }
 412     if (move_bits<0) move_bits=0;
 413
 414     if (targ_bits[1] >= 125) {
 415       /* dont reduce side channel below 125 bits */
 416       if (targ_bits[1]-move_bits > 125) {
 417
 418         /* if mid channel already has 2x more than average, dont bother */
 419         /* mean_bits = bits per granule (for both channels) */
 420         if (targ_bits[0] < mean_bits)
 421           targ_bits[0] += move_bits;
 422         targ_bits[1] -= move_bits;
 423       } else {
 424         targ_bits[0] += targ_bits[1] - 125;
 425         targ_bits[1] = 125;
 426       }
 427     }
 428
 429     move_bits=targ_bits[0]+targ_bits[1];
 430     if (move_bits > max_bits) {
 431       targ_bits[0]=(max_bits*targ_bits[0])/move_bits;
 432       targ_bits[1]=(max_bits*targ_bits[1])/move_bits;
 433     }
 434     assert (targ_bits[0] <= MAX_BITS);
 435     assert (targ_bits[1] <= MAX_BITS);
 436 }
 437
 438
 439 /**
 440  *  Robert Hegemann 2001-04-27:
 441  *  this adjusts the ATH, keeping the original noise floor
 442  *  affects the higher frequencies more than the lower ones
 443  */
 444
 445 FLOAT8 athAdjust( FLOAT8 a, FLOAT8 x, FLOAT8 athFloor )
 446 {
 447     /*  work in progress
 448      */
 449     FLOAT8 const o = 90.30873362;
 450     FLOAT8 const p = 94.82444863;
 451     FLOAT8 u = 10. * log10(x);
 452     FLOAT8 v = a*a;
 453     FLOAT8 w = 0.0;
 454     u -= athFloor;                                  // undo scaling
 455     if ( v > 1E-20 ) w = 1. + 10. * log10(v) / o;
 456     if ( w < 0  )    w = 0.;
 457     u *= w;
 458     u += athFloor + o-p;                            // redo scaling
 459
 460     return pow( 10., 0.1*u );
 461 }
 462
 463
 464 /*************************************************************************/
 465 /*            calc_xmin                                                  */
 466 /*************************************************************************/
 467
 468 /*
 469   Calculate the allowed distortion for each scalefactor band,
 470   as determined by the psychoacoustic model.
 471   xmin(sb) = ratio(sb) * en(sb) / bw(sb)
 472
 473   returns number of sfb's with energy > ATH
 474 */
 475 int calc_xmin(
 476         lame_global_flags *gfp,
 477         const FLOAT8                xr [576],
 478         const III_psy_ratio * const ratio,
 479         const gr_info       * const cod_info,
 480               III_psy_xmin  * const l3_xmin )
 481 {
 482     lame_internal_flags *gfc = gfp->internal_flags;
 483     int sfb,j,start, end, bw,l, b, ath_over=0;
 484     FLOAT8 en0, xmin, ener, tmpATH;
 485     ATH_t * ATH = gfc->ATH;
 486
 487     if (cod_info->block_type == SHORT_TYPE) {
 488
 489         for ( j = 0, sfb = 0; sfb < SBMAX_s; sfb++ ) {
 490             if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh )
 491                 tmpATH = athAdjust( ATH->adjust, ATH->s[sfb], ATH->floor );
 492             else
 493                 tmpATH = ATH->adjust * ATH->s[sfb];
 494             start = gfc->scalefac_band.s[ sfb ];
 495             end   = gfc->scalefac_band.s[ sfb + 1 ];
 496             bw = end - start;
 497
 498             for ( b = 0; b < 3; b++ ) {
 499                 for (en0 = 0.0, l = start; l < end; l++) {
 500                     ener = xr[j++];
 501                     ener = ener * ener;
 502                     en0 += ener;
 503                 }
 504                 en0 /= bw;
 505
 506                 if (gfp->ATHonly || gfp->ATHshort) {
 507                     xmin = tmpATH;
 508                 }
 509                 else {
 510                     xmin = ratio->en.s[sfb][b];
 511                     if (xmin > 0.0)
 512                         xmin = en0 * ratio->thm.s[sfb][b] * gfc->masking_lower / xmin;
 513                     if (xmin < tmpATH)
 514                         xmin = tmpATH;
 515                 }
 516                 xmin *= bw;
 517
 518                 if (gfc->nsPsy.use) {
 519                     if      (sfb <=  5) xmin *= gfc->nsPsy.bass;
 520                     else if (sfb <= 10) xmin *= gfc->nsPsy.alto;
 521                     else                xmin *= gfc->nsPsy.treble;
 522                     if ((gfp->VBR == vbr_off || gfp->VBR == vbr_abr) && gfp->quality <= 1)
 523                         xmin *= 0.001;
 524                 }
 525                 l3_xmin->s[sfb][b] = xmin;
 526                 if (en0 > tmpATH) ath_over++;
 527             }   /* b */
 528         }   /* sfb */
 529
 530         if (gfp->useTemporal) {
 531             for (sfb = 0; sfb < SBMAX_s; sfb++ ) {
 532                 for ( b = 1; b < 3; b++ ) {
 533                     xmin = l3_xmin->s[sfb][b] * (1.0 - gfc->decay)
 534                          + l3_xmin->s[sfb][b-1] * gfc->decay;
 535                     if (l3_xmin->s[sfb][b] < xmin)
 536                         l3_xmin->s[sfb][b] = xmin;
 537                 }
 538             }
 539         }
 540
 541     }   /* end of short block case */
 542     else {
 543
 544         for ( sfb = 0; sfb < SBMAX_l; sfb++ ){
 545             if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh )
 546                 tmpATH = athAdjust( ATH->adjust, ATH->l[sfb], ATH->floor );
 547             else
 548                 tmpATH = ATH->adjust * ATH->l[sfb];
 549             start = gfc->scalefac_band.l[ sfb ];
 550             end   = gfc->scalefac_band.l[ sfb+1 ];
 551             bw = end - start;
 552
 553             for (en0 = 0.0, l = start; l < end; l++ ) {
 554                 ener = xr[l] * xr[l];
 555                 en0 += ener;
 556             }
 557             /* why is it different from short blocks <?> */
 558             if ( !gfc->nsPsy.use ) en0 /= bw;
 559
 560             if (gfp->ATHonly) {
 561                 xmin = tmpATH;
 562             }
 563             else {
 564                 xmin = ratio->en.l[sfb];
 565                 if (xmin > 0.0)
 566                     xmin = en0 * ratio->thm.l[sfb] * gfc->masking_lower / xmin;
 567                 if (xmin < tmpATH)
 568                     xmin = tmpATH;
 569             }
 570             /* why is it different from short blocks <?> */
 571             if ( !gfc->nsPsy.use ) {
 572                 xmin *= bw;
 573             }
 574             else {
 575                 if      (sfb <=  6) xmin *= gfc->nsPsy.bass;
 576                 else if (sfb <= 13) xmin *= gfc->nsPsy.alto;
 577                 else if (sfb <= 20) xmin *= gfc->nsPsy.treble;
 578                 else                xmin *= gfc->nsPsy.sfb21;
 579                 if ((gfp->VBR == vbr_off || gfp->VBR == vbr_abr) && gfp->quality <= 1)
 580                     xmin *= 0.001;
 581             }
 582             l3_xmin->l[sfb] = xmin;
 583             if (en0 > tmpATH) ath_over++;
 584         }   /* sfb */
 585
 586     }   /* end of long block case */
 587
 588     return ath_over;
 589 }
 590
 591
 592
 593
 594
 595
 596
 597 /*************************************************************************/
 598 /*            calc_noise                                                 */
 599 /*************************************************************************/
 600
 601 // -oo dB  =>  -1.00
 602 // - 6 dB  =>  -0.97
 603 // - 3 dB  =>  -0.80
 604 // - 2 dB  =>  -0.64
 605 // - 1 dB  =>  -0.38
 606 //   0 dB  =>   0.00
 607 // + 1 dB  =>  +0.49
 608 // + 2 dB  =>  +1.06
 609 // + 3 dB  =>  +1.68
 610 // + 6 dB  =>  +3.69
 611 // +10 dB  =>  +6.45
 612
 613 double penalties ( double noise )
 614 {
 615     return log ( 0.368 + 0.632 * noise * noise * noise );
 616 }
 617
 618 /*  mt 5/99:  Function: Improved calc_noise for a single channel   */
 619
 620 int  calc_noise(
 621         const lame_internal_flags           * const gfc,
 622         const FLOAT8                    xr [576],
 623         const int                       ix [576],
 624         const gr_info           * const cod_info,
 625         const III_psy_xmin      * const l3_xmin,
 626         const III_scalefac_t    * const scalefac,
 627               III_psy_xmin      * xfsf,
 628               calc_noise_result * const res )
 629 {
 630   int sfb,start, end, j,l, i, over=0;
 631   FLOAT8 sum;
 632
 633   int count=0;
 634   FLOAT8 noise,noise_db;
 635   FLOAT8 over_noise_db = 0;
 636   FLOAT8 tot_noise_db  = 0;     /*    0 dB relative to masking */
 637   FLOAT8 max_noise  = 1E-20; /* -200 dB relative to masking */
 638   double klemm_noise = 1E-37;
 639
 640   if (cod_info->block_type == SHORT_TYPE) {
 641     int max_index = gfc->sfb21_extra ? SBMAX_s : SBPSY_s;
 642
 643     for ( j=0, sfb = 0; sfb < max_index; sfb++ ) {
 644          start = gfc->scalefac_band.s[ sfb ];
 645          end   = gfc->scalefac_band.s[ sfb+1 ];
 646          for ( i = 0; i < 3; i++ ) {
 647             FLOAT8 step;
 648             int s;
 649
 650             s = (scalefac->s[sfb][i] << (cod_info->scalefac_scale + 1))
 651                 + cod_info->subblock_gain[i] * 8;
 652             s = cod_info->global_gain - s;
 653
 654             assert(s<Q_MAX);
 655             assert(s>=0);
 656             step = POW20(s);
 657             sum = 0.0;
 658             l = start;
 659             do {
 660               FLOAT8 temp;
 661               temp = pow43[ix[j]];
 662               temp *= step;
 663               temp -= fabs(xr[j]);
 664               ++j;
 665               sum += temp * temp;
 666               l++;
 667             } while (l < end);
 668             noise = xfsf->s[sfb][i]  = sum / l3_xmin->s[sfb][i];
 669
 670             max_noise    = Max(max_noise,noise);
 671             klemm_noise += penalties (noise);
 672
 673             noise_db=10*log10(Max(noise,1E-20));
 674             tot_noise_db += noise_db;
 675
 676             if (noise > 1) {
 677                 over++;
 678                 over_noise_db += noise_db;
 679             }
 680             count++;
 681         }
 682     }
 683   }else{ /* cod_info->block_type == SHORT_TYPE */
 684     int max_index = gfc->sfb21_extra ? SBMAX_l : SBPSY_l;
 685
 686     for ( sfb = 0; sfb < max_index; sfb++ ) {
 687         FLOAT8 step;
 688         int s = scalefac->l[sfb];
 689
 690         if (cod_info->preflag)
 691             s += pretab[sfb];
 692
 693         s = cod_info->global_gain - (s << (cod_info->scalefac_scale + 1));
 694         assert(s<Q_MAX);
 695         assert(s>=0);
 696         step = POW20(s);
 697
 698         start = gfc->scalefac_band.l[ sfb ];
 699         end   = gfc->scalefac_band.l[ sfb+1 ];
 700
 701         for ( sum = 0.0, l = start; l < end; l++ ) {
 702           FLOAT8 temp;
 703           temp = fabs(xr[l]) - pow43[ix[l]] * step;
 704           sum += temp * temp;
 705         }
 706         noise = xfsf->l[sfb] = sum / l3_xmin->l[sfb];
 707         max_noise=Max(max_noise,noise);
 708         klemm_noise += penalties (noise);
 709
 710         noise_db=10*log10(Max(noise,1E-20));
 711         /* multiplying here is adding in dB, but can overflow */
 712         //tot_noise *= Max(noise, 1E-20);
 713         tot_noise_db += noise_db;
 714
 715         if (noise > 1) {
 716           over++;
 717           /* multiplying here is adding in dB -but can overflow */
 718           //over_noise *= noise;
 719           over_noise_db += noise_db;
 720         }
 721
 722         count++;
 723     }
 724   } /* cod_info->block_type == SHORT_TYPE */
 725
 726   /* normalization at this point by "count" is not necessary, since
 727    * the values are only used to compare with previous values */
 728   res->over_count = over;
 729
 730   /* convert to db. DO NOT CHANGE THESE */
 731   /* tot_noise = is really the average over each sfb of:
 732      [noise(db) - allowed_noise(db)]
 733
 734      and over_noise is the same average, only over only the
 735      bands with noise > allowed_noise.
 736
 737   */
 738
 739   //res->tot_noise   = 10.*log10(Max(1e-20,tot_noise ));
 740   //res->over_noise  = 10.*log10(Max(1e+00,over_noise));
 741   res->tot_noise   = tot_noise_db;
 742   res->over_noise  = over_noise_db;
 743   res->max_noise   = 10.*log10(Max(1e-20,max_noise ));
 744   res->klemm_noise = 10.*log10(Max(1e-20,klemm_noise));
 745
 746   return over;
 747 }
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762 /************************************************************************
 763  *
 764  *  set_pinfo()
 765  *
 766  *  updates plotting data
 767  *
 768  *  Mark Taylor 2000-??-??
 769  *
 770  *  Robert Hegemann: moved noise/distortion calc into it
 771  *
 772  ************************************************************************/
 773
 774 static
 775 void set_pinfo (
 776         lame_global_flags *gfp,
 777         const gr_info        * const cod_info,
 778         const III_psy_ratio  * const ratio,
 779         const III_scalefac_t * const scalefac,
 780         const FLOAT8                 xr[576],
 781         const int                    l3_enc[576],
 782         const int                    gr,
 783         const int                    ch )
 784 {
 785     lame_internal_flags *gfc=gfp->internal_flags;
 786     int sfb;
 787     int j,i,l,start,end,bw;
 788     FLOAT8 en0,en1;
 789     FLOAT ifqstep = ( cod_info->scalefac_scale == 0 ) ? .5 : 1.0;
 790
 791
 792     III_psy_xmin l3_xmin;
 793     calc_noise_result noise;
 794     III_psy_xmin xfsf;
 795
 796     calc_xmin (gfp,xr, ratio, cod_info, &l3_xmin);
 797
 798     calc_noise (gfc, xr, l3_enc, cod_info, &l3_xmin, scalefac, &xfsf, &noise);
 799
 800     if (cod_info->block_type == SHORT_TYPE) {
 801         for (j=0, sfb = 0; sfb < SBMAX_s; sfb++ )  {
 802             start = gfc->scalefac_band.s[ sfb ];
 803             end   = gfc->scalefac_band.s[ sfb + 1 ];
 804             bw = end - start;
 805             for ( i = 0; i < 3; i++ ) {
 806                 for ( en0 = 0.0, l = start; l < end; l++ ) {
 807                     en0 += xr[j] * xr[j];
 808                     ++j;
 809                 }
 810                 en0=Max(en0/bw,1e-20);
 811
 812
 813 #if 0
 814 {
 815     double tot1,tot2;
 816     if (sfb<SBMAX_s-1) {
 817         if (sfb==0) {
 818             tot1=0;
 819             tot2=0;
 820         }
 821         tot1 += en0;
 822         tot2 += ratio->en.s[sfb][i];
 823
 824         DEBUGF("%i %i sfb=%i mdct=%f fft=%f  fft-mdct=%f db \n",
 825                 gr,ch,sfb,
 826                 10*log10(Max(1e-25,ratio->en.s[sfb][i])),
 827                 10*log10(Max(1e-25,en0)),
 828                 10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.s[sfb][i]))));
 829
 830         if (sfb==SBMAX_s-2) {
 831             DEBUGF("%i %i toti %f %f ratio=%f db \n",gr,ch,
 832                     10*log10(Max(1e-25,tot2)),
 833                     10*log(Max(1e-25,tot1)),
 834                     10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
 835
 836         }
 837     }
 838     /*
 839         masking: multiplied by en0/fft_energy
 840         average seems to be about -135db.
 841      */
 842 }
 843 #endif
 844
 845
 846                 /* convert to MDCT units */
 847                 en1=1e15;  /* scaling so it shows up on FFT plot */
 848                 gfc->pinfo->xfsf_s[gr][ch][3*sfb+i]
 849                     = en1*xfsf.s[sfb][i]*l3_xmin.s[sfb][i]/bw;
 850                 gfc->pinfo->en_s[gr][ch][3*sfb+i] = en1*en0;
 851
 852                 if (ratio->en.s[sfb][i]>0)
 853                     en0 = en0/ratio->en.s[sfb][i];
 854                 else
 855                     en0=0;
 856                 if (gfp->ATHonly || gfp->ATHshort)
 857                     en0=0;
 858
 859                 gfc->pinfo->thr_s[gr][ch][3*sfb+i] =
 860                         en1*Max(en0*ratio->thm.s[sfb][i],gfc->ATH->s[sfb]);
 861
 862
 863                 /* there is no scalefactor bands >= SBPSY_s */
 864                 if (sfb < SBPSY_s) {
 865                     gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=
 866                                             -ifqstep*scalefac->s[sfb][i];
 867                 } else {
 868                     gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=0;
 869                 }
 870                 gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i] -=
 871                                              2*cod_info->subblock_gain[i];
 872             }
 873         }
 874     } else {
 875         for ( sfb = 0; sfb < SBMAX_l; sfb++ )   {
 876             start = gfc->scalefac_band.l[ sfb ];
 877             end   = gfc->scalefac_band.l[ sfb+1 ];
 878             bw = end - start;
 879             for ( en0 = 0.0, l = start; l < end; l++ )
 880                 en0 += xr[l] * xr[l];
 881             if (!gfc->nsPsy.use) en0/=bw;
 882       /*
 883     DEBUGF("diff  = %f \n",10*log10(Max(ratio[gr][ch].en.l[sfb],1e-20))
 884                             -(10*log10(en0)+150));
 885        */
 886
 887 #if 0
 888  {
 889     double tot1,tot2;
 890     if (sfb==0) {
 891         tot1=0;
 892         tot2=0;
 893     }
 894     tot1 += en0;
 895     tot2 += ratio->en.l[sfb];
 896
 897
 898     DEBUGF("%i %i sfb=%i mdct=%f fft=%f  fft-mdct=%f db \n",
 899             gr,ch,sfb,
 900             10*log10(Max(1e-25,ratio->en.l[sfb])),
 901             10*log10(Max(1e-25,en0)),
 902             10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.l[sfb]))));
 903
 904     if (sfb==SBMAX_l-1) {
 905         DEBUGF("%i %i toti %f %f ratio=%f db \n",
 906             gr,ch,
 907             10*log10(Max(1e-25,tot2)),
 908             10*log(Max(1e-25,tot1)),
 909             10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
 910     }
 911     /*
 912         masking: multiplied by en0/fft_energy
 913         average seems to be about -147db.
 914      */
 915 }
 916 #endif
 917
 918
 919             /* convert to MDCT units */
 920             en1=1e15;  /* scaling so it shows up on FFT plot */
 921             gfc->pinfo->xfsf[gr][ch][sfb] =  en1*xfsf.l[sfb]*l3_xmin.l[sfb]/bw;
 922             gfc->pinfo->en[gr][ch][sfb] = en1*en0;
 923             if (ratio->en.l[sfb]>0)
 924                 en0 = en0/ratio->en.l[sfb];
 925             else
 926                 en0=0;
 927             if (gfp->ATHonly)
 928                 en0=0;
 929             gfc->pinfo->thr[gr][ch][sfb] =
 930                              en1*Max(en0*ratio->thm.l[sfb],gfc->ATH->l[sfb]);
 931
 932             /* there is no scalefactor bands >= SBPSY_l */
 933             if (sfb<SBPSY_l) {
 934                 if (scalefac->l[sfb]<0)  /* scfsi! */
 935                     gfc->pinfo->LAMEsfb[gr][ch][sfb] =
 936                                             gfc->pinfo->LAMEsfb[0][ch][sfb];
 937                 else
 938                     gfc->pinfo->LAMEsfb[gr][ch][sfb] = -ifqstep*scalefac->l[sfb];
 939             }else{
 940                 gfc->pinfo->LAMEsfb[gr][ch][sfb] = 0;
 941             }
 942
 943             if (cod_info->preflag && sfb>=11)
 944                 gfc->pinfo->LAMEsfb[gr][ch][sfb] -= ifqstep*pretab[sfb];
 945         } /* for sfb */
 946     } /* block type long */
 947     gfc->pinfo->LAMEqss     [gr][ch] = cod_info->global_gain;
 948     gfc->pinfo->LAMEmainbits[gr][ch] = cod_info->part2_3_length;
 949     gfc->pinfo->LAMEsfbits  [gr][ch] = cod_info->part2_length;
 950
 951     gfc->pinfo->over      [gr][ch] = noise.over_count;
 952     gfc->pinfo->max_noise [gr][ch] = noise.max_noise;
 953     gfc->pinfo->over_noise[gr][ch] = noise.over_noise;
 954     gfc->pinfo->tot_noise [gr][ch] = noise.tot_noise;
 955 }
 956
 957
 958 /************************************************************************
 959  *
 960  *  set_frame_pinfo()
 961  *
 962  *  updates plotting data for a whole frame
 963  *
 964  *  Robert Hegemann 2000-10-21
 965  *
 966  ************************************************************************/
 967
 968 void set_frame_pinfo(
 969         lame_global_flags *gfp,
 970         FLOAT8          xr       [2][2][576],
 971         III_psy_ratio   ratio    [2][2],
 972         int             l3_enc   [2][2][576],
 973         III_scalefac_t  scalefac [2][2] )
 974 {
 975     lame_internal_flags *gfc=gfp->internal_flags;
 976     unsigned int          sfb;
 977         int                   ch;
 978         int                   gr;
 979     int                   act_l3enc[576];
 980     III_scalefac_t        act_scalefac [2];
 981     int scsfi[2] = {0,0};
 982
 983
 984     gfc->masking_lower = 1.0;
 985
 986     /* reconstruct the scalefactors in case SCSFI was used
 987      */
 988     for (ch = 0; ch < gfc->channels_out; ch ++) {
 989         for (sfb = 0; sfb < SBMAX_l; sfb ++) {
 990             if (scalefac[1][ch].l[sfb] == -1) {/* scfsi */
 991                 act_scalefac[ch].l[sfb] = scalefac[0][ch].l[sfb];
 992                 scsfi[ch] = 1;
 993             } else {
 994                 act_scalefac[ch].l[sfb] = scalefac[1][ch].l[sfb];
 995             }
 996         }
 997     }
 998
 999     /* for every granule and channel patch l3_enc and set info
1000      */
1001     for (gr = 0; gr < gfc->mode_gr; gr ++) {
1002         for (ch = 0; ch < gfc->channels_out; ch ++) {
1003             int i;
1004             gr_info *cod_info = &gfc->l3_side.gr[gr].ch[ch].tt;
1005
1006             /* revert back the sign of l3enc */
1007             for ( i = 0; i < 576; i++) {
1008                 if (xr[gr][ch][i] < 0)
1009                     act_l3enc[i] = -l3_enc[gr][ch][i];
1010                 else
1011                     act_l3enc[i] = +l3_enc[gr][ch][i];
1012             }
1013             if (gr == 1 && scsfi[ch])
1014                 set_pinfo (gfp, cod_info, &ratio[gr][ch], &act_scalefac[ch],
1015                         xr[gr][ch], act_l3enc, gr, ch);
1016             else
1017                 set_pinfo (gfp, cod_info, &ratio[gr][ch], &scalefac[gr][ch],
1018                         xr[gr][ch], act_l3enc, gr, ch);
1019         } /* for ch */
1020     }    /* for gr */
1021 }
1022
1023
1024
1025
1026 /*********************************************************************
1027  * nonlinear quantization of xr
1028  * More accurate formula than the ISO formula.  Takes into account
1029  * the fact that we are quantizing xr -> ix, but we want ix^4/3 to be
1030  * as close as possible to x^4/3.  (taking the nearest int would mean
1031  * ix is as close as possible to xr, which is different.)
1032  * From Segher Boessenkool <segher@eastsite.nl>  11/1999
1033  * ASM optimization from
1034  *    Mathew Hendry <scampi@dial.pipex.com> 11/1999
1035  *    Acy Stapp <AStapp@austin.rr.com> 11/1999
1036  *    Takehiro Tominaga <tominaga@isoternet.org> 11/1999
1037  * 9/00: ASM code removed in favor of IEEE754 hack.  If you need
1038  * the ASM code, check CVS circa Aug 2000.
1039  *********************************************************************/
1040
1041
1042 #ifdef TAKEHIRO_IEEE754_HACK
1043
1044 typedef union {
1045     float f;
1046     int i;
1047 } fi_union;
1048
1049 #define MAGIC_FLOAT (65536*(128))
1050 #define MAGIC_INT 0x4b000000
1051
1052 void quantize_xrpow(const FLOAT8 *xp, int *pi, FLOAT8 istep)
1053 {
1054     /* quantize on xr^(3/4) instead of xr */
1055     int j;
1056     fi_union *fi;
1057
1058     fi = (fi_union *)pi;
1059     for (j = 576 / 4 - 1; j >= 0; --j) {
1060         double x0 = istep * xp[0];
1061         double x1 = istep * xp[1];
1062         double x2 = istep * xp[2];
1063         double x3 = istep * xp[3];
1064
1065         x0 += MAGIC_FLOAT; fi[0].f = x0;
1066         x1 += MAGIC_FLOAT; fi[1].f = x1;
1067         x2 += MAGIC_FLOAT; fi[2].f = x2;
1068         x3 += MAGIC_FLOAT; fi[3].f = x3;
1069
1070         fi[0].f = x0 + (adj43asm - MAGIC_INT)[fi[0].i];
1071         fi[1].f = x1 + (adj43asm - MAGIC_INT)[fi[1].i];
1072         fi[2].f = x2 + (adj43asm - MAGIC_INT)[fi[2].i];
1073         fi[3].f = x3 + (adj43asm - MAGIC_INT)[fi[3].i];
1074
1075         fi[0].i -= MAGIC_INT;
1076         fi[1].i -= MAGIC_INT;
1077         fi[2].i -= MAGIC_INT;
1078         fi[3].i -= MAGIC_INT;
1079         fi += 4;
1080         xp += 4;
1081     }
1082 }
1083
1084 #  define ROUNDFAC -0.0946
1085 void quantize_xrpow_ISO(const FLOAT8 *xp, int *pi, FLOAT8 istep)
1086 {
1087     /* quantize on xr^(3/4) instead of xr */
1088     int j;
1089     fi_union *fi;
1090
1091     fi = (fi_union *)pi;
1092     for (j=576/4 - 1;j>=0;j--) {
1093         fi[0].f = istep * xp[0] + (ROUNDFAC + MAGIC_FLOAT);
1094         fi[1].f = istep * xp[1] + (ROUNDFAC + MAGIC_FLOAT);
1095         fi[2].f = istep * xp[2] + (ROUNDFAC + MAGIC_FLOAT);
1096         fi[3].f = istep * xp[3] + (ROUNDFAC + MAGIC_FLOAT);
1097
1098         fi[0].i -= MAGIC_INT;
1099         fi[1].i -= MAGIC_INT;
1100         fi[2].i -= MAGIC_INT;
1101         fi[3].i -= MAGIC_INT;
1102         fi+=4;
1103         xp+=4;
1104     }
1105 }
1106
1107 #else
1108
1109 /*********************************************************************
1110  * XRPOW_FTOI is a macro to convert floats to ints.
1111  * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x]
1112  *                                         ROUNDFAC= -0.0946
1113  *
1114  * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]
1115  *                                   ROUNDFAC=0.4054
1116  *
1117  * Note: using floor() or (int) is extermely slow. On machines where
1118  * the TAKEHIRO_IEEE754_HACK code above does not work, it is worthwile
1119  * to write some ASM for XRPOW_FTOI().
1120  *********************************************************************/
1121 #define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
1122 #define QUANTFAC(rx)  adj43[rx]
1123 #define ROUNDFAC 0.4054
1124
1125
1126 void quantize_xrpow(const FLOAT8 *xr, int *ix, FLOAT8 istep) {
1127     /* quantize on xr^(3/4) instead of xr */
1128     /* from Wilfried.Behne@t-online.de.  Reported to be 2x faster than
1129        the above code (when not using ASM) on PowerPC */
1130     int j;
1131
1132     for ( j = 576/8; j > 0; --j) {
1133         FLOAT8  x1, x2, x3, x4, x5, x6, x7, x8;
1134         int     rx1, rx2, rx3, rx4, rx5, rx6, rx7, rx8;
1135         x1 = *xr++ * istep;
1136         x2 = *xr++ * istep;
1137         XRPOW_FTOI(x1, rx1);
1138         x3 = *xr++ * istep;
1139         XRPOW_FTOI(x2, rx2);
1140         x4 = *xr++ * istep;
1141         XRPOW_FTOI(x3, rx3);
1142         x5 = *xr++ * istep;
1143         XRPOW_FTOI(x4, rx4);
1144         x6 = *xr++ * istep;
1145         XRPOW_FTOI(x5, rx5);
1146         x7 = *xr++ * istep;
1147         XRPOW_FTOI(x6, rx6);
1148         x8 = *xr++ * istep;
1149         XRPOW_FTOI(x7, rx7);
1150         x1 += QUANTFAC(rx1);
1151         XRPOW_FTOI(x8, rx8);
1152         x2 += QUANTFAC(rx2);
1153         XRPOW_FTOI(x1,*ix++);
1154         x3 += QUANTFAC(rx3);
1155         XRPOW_FTOI(x2,*ix++);
1156         x4 += QUANTFAC(rx4);
1157         XRPOW_FTOI(x3,*ix++);
1158         x5 += QUANTFAC(rx5);
1159         XRPOW_FTOI(x4,*ix++);
1160         x6 += QUANTFAC(rx6);
1161         XRPOW_FTOI(x5,*ix++);
1162         x7 += QUANTFAC(rx7);
1163         XRPOW_FTOI(x6,*ix++);
1164         x8 += QUANTFAC(rx8);
1165         XRPOW_FTOI(x7,*ix++);
1166         XRPOW_FTOI(x8,*ix++);
1167     }
1168 }
1169
1170
1171
1172
1173
1174
1175 void quantize_xrpow_ISO( const FLOAT8 *xr, int *ix, FLOAT8 istep )
1176 {
1177     /* quantize on xr^(3/4) instead of xr */
1178     const FLOAT8 compareval0 = (1.0 - 0.4054)/istep;
1179     int j;
1180     /* depending on architecture, it may be worth calculating a few more
1181        compareval's.
1182
1183        eg.  compareval1 = (2.0 - 0.4054/istep);
1184        .. and then after the first compare do this ...
1185        if compareval1>*xr then ix = 1;
1186
1187        On a pentium166, it's only worth doing the one compare (as done here),
1188        as the second compare becomes more expensive than just calculating
1189        the value. Architectures with slow FP operations may want to add some
1190        more comparevals. try it and send your diffs statistically speaking
1191
1192        73% of all xr*istep values give ix=0
1193        16% will give 1
1194        4%  will give 2
1195     */
1196     for (j=576;j>0;j--) {
1197         if (compareval0 > *xr) {
1198             *(ix++) = 0;
1199             xr++;
1200         } else {
1201             /*    *(ix++) = (int)( istep*(*(xr++))  + 0.4054); */
1202             XRPOW_FTOI(  istep*(*(xr++))  + ROUNDFAC , *(ix++) );
1203         }
1204     }
1205 }
1206
1207 #endif