Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

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