2 * quantize_pvt source file
4 * Copyright (c) 1999 Takehiro TOMINAGA
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.
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.
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.
22 /* $Id: quantize_pvt.c,v 1.1 2002/04/28 17:30:26 kramm Exp $ */
23 #include "config_static.h"
27 #include "lame-analysis.h"
29 #include "reservoir.h"
30 #include "quantize_pvt.h"
37 #define NSATHSCALE 100 // Assuming dynamic range=96dB, this value should be 92
39 const char slen1_tab [16] = { 0, 0, 0, 0, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4 };
40 const char slen2_tab [16] = { 0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3, 1, 2, 3, 2, 3 };
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:
49 [table_number][row_in_table][column of nr_of_sfb]
51 const int nr_of_sfb_block [6] [3] [4] =
86 /* Table B.6: layer3 preemphasis */
87 const char pretab [SBMAX_l] =
89 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
90 1, 1, 1, 1, 2, 2, 3, 3, 3, 2, 0
94 Here are MPEG1 Table B.8 and MPEG2 Table B.1
95 -- Layer III scalefactor bands.
96 Index into this using a method such as:
97 idx = fr_ps->header->sampling_frequency
98 + (fr_ps->header->version * 3)
102 const scalefac_struct sfBandIndex[9] =
104 { /* Table B.2.b: 22.05 kHz */
105 {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
106 {0,4,8,12,18,24,32,42,56,74,100,132,174,192}
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,540,576},
110 {0,4,8,12,18,26,36,48,62,80,104,136,180,192}
112 { /* Table B.2.a: 16 kHz */
113 {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
114 {0,4,8,12,18,26,36,48,62,80,104,134,174,192}
116 { /* Table B.8.b: 44.1 kHz */
117 {0,4,8,12,16,20,24,30,36,44,52,62,74,90,110,134,162,196,238,288,342,418,576},
118 {0,4,8,12,16,22,30,40,52,66,84,106,136,192}
120 { /* Table B.8.c: 48 kHz */
121 {0,4,8,12,16,20,24,30,36,42,50,60,72,88,106,128,156,190,230,276,330,384,576},
122 {0,4,8,12,16,22,28,38,50,64,80,100,126,192}
124 { /* Table B.8.a: 32 kHz */
125 {0,4,8,12,16,20,24,30,36,44,54,66,82,102,126,156,194,240,296,364,448,550,576},
126 {0,4,8,12,16,22,30,42,58,78,104,138,180,192}
128 { /* MPEG-2.5 11.025 kHz */
129 {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
130 {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 { /* MPEG-2.5 12 kHz */
133 {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
134 {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 { /* MPEG-2.5 8 kHz */
137 {0,12,24,36,48,60,72,88,108,132,160,192,232,280,336,400,476,566,568,570,572,574,576},
138 {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}
145 FLOAT8 ipow20[Q_MAX];
146 FLOAT8 iipow20[Q_MAX];
148 FLOAT8 pow43[PRECALC_SIZE];
149 /* initialized in first call to iteration_init */
150 FLOAT8 adj43asm[PRECALC_SIZE];
151 FLOAT8 adj43[PRECALC_SIZE];
153 /************************************************************************/
154 /* initialization for iteration_loop */
155 /************************************************************************/
157 iteration_init( lame_global_flags *gfp)
159 lame_internal_flags *gfc=gfp->internal_flags;
160 III_side_info_t * const l3_side = &gfc->l3_side;
163 if ( gfc->iteration_init_init==0 ) {
164 gfc->iteration_init_init=1;
166 l3_side->main_data_begin = 0;
167 compute_ath(gfp,gfc->ATH->l,gfc->ATH->s);
170 for(i=1;i<PRECALC_SIZE;i++)
171 pow43[i] = pow((FLOAT8)i, 4.0/3.0);
174 for (i = 1; i < PRECALC_SIZE; i++)
175 adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]),0.75);
176 for (i = 0; i < PRECALC_SIZE-1; i++)
177 adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
179 iipow20_ = &iipow20[210];
180 for (i = 0; i < Q_MAX; i++) {
181 iipow20[i] = pow(2.0, (double)(i - 210) * 0.1875);
182 ipow20[i] = pow(2.0, (double)(i - 210) * -0.1875);
183 pow20[i] = pow(2.0, (double)(i - 210) * 0.25);
194 compute the ATH for each scalefactor band
197 Input: 3.3kHz signal 32767 amplitude (3.3kHz is where ATH is smallest = -5db)
198 longblocks: sfb=12 en0/bw=-11db max_en0 = 1.3db
199 shortblocks: sfb=5 -9db 0db
201 Input: 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
202 longblocks: amp=1 sfb=12 en0/bw=-103 db max_en0 = -92db
203 amp=32767 sfb=12 -12 db -1.4db
205 Input: 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
206 shortblocks: amp=1 sfb=5 en0/bw= -99 -86
207 amp=32767 sfb=5 -9 db 4db
210 MAX energy of largest wave at 3.3kHz = 1db
211 AVE energy of largest wave at 3.3kHz = -11db
212 Let's take AVE: -11db = maximum signal in sfb=12.
213 Dynamic range of CD: 96db. Therefor energy of smallest audible wave
214 in sfb=12 = -11 - 96 = -107db = ATH at 3.3kHz.
216 ATH formula for this wave: -5db. To adjust to LAME scaling, we need
217 ATH = ATH_formula - 103 (db)
218 ATH = ATH * 2.5e-10 (ener)
222 FLOAT8 ATHmdct( lame_global_flags *gfp, FLOAT8 f )
224 lame_internal_flags *gfc = gfp->internal_flags;
227 ath = ATHformula( f , gfp );
229 if (gfc->nsPsy.use) {
235 /* modify the MDCT scaling for the ATH */
236 ath -= gfp->ATHlower;
238 /* convert to energy */
239 ath = pow( 10.0, ath/10.0 );
244 void compute_ath( lame_global_flags *gfp, FLOAT8 ATH_l[], FLOAT8 ATH_s[] )
246 lame_internal_flags *gfc = gfp->internal_flags;
247 int sfb, i, start, end;
249 FLOAT8 samp_freq = gfp->out_samplerate;
251 for (sfb = 0; sfb < SBMAX_l; sfb++) {
252 start = gfc->scalefac_band.l[ sfb ];
253 end = gfc->scalefac_band.l[ sfb+1 ];
254 ATH_l[sfb]=FLOAT8_MAX;
255 for (i = start ; i < end; i++) {
256 FLOAT8 freq = i*samp_freq/(2*576);
257 ATH_f = ATHmdct( gfp, freq ); /* freq in kHz */
258 ATH_l[sfb] = Min( ATH_l[sfb], ATH_f );
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] = FLOAT8_MAX;
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 );
274 * reduce ATH to -200 dB
278 for (sfb = 0; sfb < SBMAX_l; sfb++) {
281 for (sfb = 0; sfb < SBMAX_s; sfb++) {
286 /* work in progress, don't rely on it too much
288 gfc->ATH-> floor = 10. * log10( ATHmdct( gfp, -1. ) );
291 { FLOAT8 g=10000, t=1e30, x;
292 for ( f = 100; f < 10000; f++ ) {
293 x = ATHmdct( gfp, f );
294 if ( t > x ) t = x, g = f;
296 printf("min=%g\n", g);
304 /* convert from L/R <-> Mid/Side, src == dst allowed */
305 void ms_convert(FLOAT8 dst[2][576], FLOAT8 src[2][576])
310 for (i = 0; i < 576; ++i) {
313 dst[0][i] = (l+r) * (FLOAT8)(SQRT2*0.5);
314 dst[1][i] = (l-r) * (FLOAT8)(SQRT2*0.5);
320 /************************************************************************
321 * allocate bits among 2 channels based on PE
323 * bugfixes rh 8/01: often allocated more than the allowed 4095 bits
324 ************************************************************************/
325 int on_pe( lame_global_flags *gfp, FLOAT8 pe[][2], III_side_info_t *l3_side,
326 int targ_bits[2], int mean_bits, int gr )
328 lame_internal_flags * gfc = gfp->internal_flags;
330 int extra_bits, tbits, bits;
332 int max_bits; /* maximum allowed bits for this granule */
335 /* allocate targ_bits for granule */
336 ResvMaxBits( gfp, mean_bits, &tbits, &extra_bits );
337 max_bits = tbits + extra_bits;
338 if (max_bits > MAX_BITS) /* hard limit per granule */
341 for ( bits = 0, ch = 0; ch < gfc->channels_out; ++ch ) {
342 /******************************************************************
343 * allocate bits for each channel
344 ******************************************************************/
345 cod_info = &l3_side->gr[gr].ch[ch].tt;
347 targ_bits[ch] = Min( MAX_BITS, tbits/gfc->channels_out );
349 if (gfc->nsPsy.use) {
350 add_bits[ch] = targ_bits[ch] * pe[gr][ch] / 700.0 - targ_bits[ch];
353 add_bits[ch] = (pe[gr][ch]-750) / 1.4;
354 /* short blocks us a little extra, no matter what the pe */
355 if ( cod_info->block_type == SHORT_TYPE ) {
356 if (add_bits[ch] < mean_bits/4)
357 add_bits[ch] = mean_bits/4;
361 /* at most increase bits by 1.5*average */
362 if (add_bits[ch] > mean_bits*3/4)
363 add_bits[ch] = mean_bits*3/4;
364 if (add_bits[ch] < 0)
367 if (add_bits[ch]+targ_bits[ch] > MAX_BITS)
368 add_bits[ch] = Max( 0, MAX_BITS-targ_bits[ch] );
370 bits += add_bits[ch];
372 if (bits > extra_bits) {
373 for ( ch = 0; ch < gfc->channels_out; ++ch ) {
374 add_bits[ch] = extra_bits * add_bits[ch] / bits;
378 for ( ch = 0; ch < gfc->channels_out; ++ch ) {
379 targ_bits[ch] += add_bits[ch];
380 extra_bits -= add_bits[ch];
381 assert( targ_bits[ch] <= MAX_BITS );
383 assert( max_bits <= MAX_BITS );
390 void reduce_side(int targ_bits[2],FLOAT8 ms_ener_ratio,int mean_bits,int max_bits)
396 /* ms_ener_ratio = 0: allocate 66/33 mid/side fac=.33
397 * ms_ener_ratio =.5: allocate 50/50 mid/side fac= 0 */
398 /* 75/25 split is fac=.5 */
399 /* float fac = .50*(.5-ms_ener_ratio[gr])/.5;*/
400 fac = .33*(.5-ms_ener_ratio)/.5;
404 /* number of bits to move from side channel to mid channel */
405 /* move_bits = fac*targ_bits[1]; */
406 move_bits = fac*.5*(targ_bits[0]+targ_bits[1]);
408 if (move_bits > MAX_BITS - targ_bits[0]) {
409 move_bits = MAX_BITS - targ_bits[0];
411 if (move_bits<0) move_bits=0;
413 if (targ_bits[1] >= 125) {
414 /* dont reduce side channel below 125 bits */
415 if (targ_bits[1]-move_bits > 125) {
417 /* if mid channel already has 2x more than average, dont bother */
418 /* mean_bits = bits per granule (for both channels) */
419 if (targ_bits[0] < mean_bits)
420 targ_bits[0] += move_bits;
421 targ_bits[1] -= move_bits;
423 targ_bits[0] += targ_bits[1] - 125;
428 move_bits=targ_bits[0]+targ_bits[1];
429 if (move_bits > max_bits) {
430 targ_bits[0]=(max_bits*targ_bits[0])/move_bits;
431 targ_bits[1]=(max_bits*targ_bits[1])/move_bits;
433 assert (targ_bits[0] <= MAX_BITS);
434 assert (targ_bits[1] <= MAX_BITS);
439 * Robert Hegemann 2001-04-27:
440 * this adjusts the ATH, keeping the original noise floor
441 * affects the higher frequencies more than the lower ones
444 FLOAT8 athAdjust( FLOAT8 a, FLOAT8 x, FLOAT8 athFloor )
448 FLOAT8 const o = 90.30873362;
449 FLOAT8 const p = 94.82444863;
450 FLOAT8 u = 10. * log10(x);
453 u -= athFloor; // undo scaling
454 if ( v > 1E-20 ) w = 1. + 10. * log10(v) / o;
457 u += athFloor + o-p; // redo scaling
459 return pow( 10., 0.1*u );
463 /*************************************************************************/
465 /*************************************************************************/
468 Calculate the allowed distortion for each scalefactor band,
469 as determined by the psychoacoustic model.
470 xmin(sb) = ratio(sb) * en(sb) / bw(sb)
472 returns number of sfb's with energy > ATH
475 lame_global_flags *gfp,
476 const FLOAT8 xr [576],
477 const III_psy_ratio * const ratio,
478 const gr_info * const cod_info,
479 III_psy_xmin * const l3_xmin )
481 lame_internal_flags *gfc = gfp->internal_flags;
482 int sfb,j,start, end, bw,l, b, ath_over=0;
483 FLOAT8 en0, xmin, ener, tmpATH;
484 ATH_t * ATH = gfc->ATH;
486 if (cod_info->block_type == SHORT_TYPE) {
488 for ( j = 0, sfb = 0; sfb < SBMAX_s; sfb++ ) {
489 if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh )
490 tmpATH = athAdjust( ATH->adjust, ATH->s[sfb], ATH->floor );
492 tmpATH = ATH->adjust * ATH->s[sfb];
493 start = gfc->scalefac_band.s[ sfb ];
494 end = gfc->scalefac_band.s[ sfb + 1 ];
497 for ( b = 0; b < 3; b++ ) {
498 for (en0 = 0.0, l = start; l < end; l++) {
505 if (gfp->ATHonly || gfp->ATHshort) {
509 xmin = ratio->en.s[sfb][b];
511 xmin = en0 * ratio->thm.s[sfb][b] * gfc->masking_lower / xmin;
517 if (gfc->nsPsy.use) {
518 if (sfb <= 5) xmin *= gfc->nsPsy.bass;
519 else if (sfb <= 10) xmin *= gfc->nsPsy.alto;
520 else xmin *= gfc->nsPsy.treble;
521 if ((gfp->VBR == vbr_off || gfp->VBR == vbr_abr) && gfp->quality <= 1)
524 l3_xmin->s[sfb][b] = xmin;
525 if (en0 > tmpATH) ath_over++;
529 if (gfp->useTemporal) {
530 for (sfb = 0; sfb < SBMAX_s; sfb++ ) {
531 for ( b = 1; b < 3; b++ ) {
532 xmin = l3_xmin->s[sfb][b] * (1.0 - gfc->decay)
533 + l3_xmin->s[sfb][b-1] * gfc->decay;
534 if (l3_xmin->s[sfb][b] < xmin)
535 l3_xmin->s[sfb][b] = xmin;
540 } /* end of short block case */
543 for ( sfb = 0; sfb < SBMAX_l; sfb++ ){
544 if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh )
545 tmpATH = athAdjust( ATH->adjust, ATH->l[sfb], ATH->floor );
547 tmpATH = ATH->adjust * ATH->l[sfb];
548 start = gfc->scalefac_band.l[ sfb ];
549 end = gfc->scalefac_band.l[ sfb+1 ];
552 for (en0 = 0.0, l = start; l < end; l++ ) {
553 ener = xr[l] * xr[l];
556 /* why is it different from short blocks <?> */
557 if ( !gfc->nsPsy.use ) en0 /= bw;
563 xmin = ratio->en.l[sfb];
565 xmin = en0 * ratio->thm.l[sfb] * gfc->masking_lower / xmin;
569 /* why is it different from short blocks <?> */
570 if ( !gfc->nsPsy.use ) {
574 if (sfb <= 6) xmin *= gfc->nsPsy.bass;
575 else if (sfb <= 13) xmin *= gfc->nsPsy.alto;
576 else if (sfb <= 20) xmin *= gfc->nsPsy.treble;
577 else xmin *= gfc->nsPsy.sfb21;
578 if ((gfp->VBR == vbr_off || gfp->VBR == vbr_abr) && gfp->quality <= 1)
581 l3_xmin->l[sfb] = xmin;
582 if (en0 > tmpATH) ath_over++;
585 } /* end of long block case */
596 /*************************************************************************/
598 /*************************************************************************/
612 double penalties ( double noise )
614 return log ( 0.368 + 0.632 * noise * noise * noise );
617 /* mt 5/99: Function: Improved calc_noise for a single channel */
620 const lame_internal_flags * const gfc,
621 const FLOAT8 xr [576],
623 const gr_info * const cod_info,
624 const III_psy_xmin * const l3_xmin,
625 const III_scalefac_t * const scalefac,
627 calc_noise_result * const res )
629 int sfb,start, end, j,l, i, over=0;
633 FLOAT8 noise,noise_db;
634 FLOAT8 over_noise_db = 0;
635 FLOAT8 tot_noise_db = 0; /* 0 dB relative to masking */
636 FLOAT8 max_noise = 1E-20; /* -200 dB relative to masking */
637 double klemm_noise = 1E-37;
639 if (cod_info->block_type == SHORT_TYPE) {
640 int max_index = gfc->sfb21_extra ? SBMAX_s : SBPSY_s;
642 for ( j=0, sfb = 0; sfb < max_index; sfb++ ) {
643 start = gfc->scalefac_band.s[ sfb ];
644 end = gfc->scalefac_band.s[ sfb+1 ];
645 for ( i = 0; i < 3; i++ ) {
649 s = (scalefac->s[sfb][i] << (cod_info->scalefac_scale + 1))
650 + cod_info->subblock_gain[i] * 8;
651 s = cod_info->global_gain - s;
667 noise = xfsf->s[sfb][i] = sum / l3_xmin->s[sfb][i];
669 max_noise = Max(max_noise,noise);
670 klemm_noise += penalties (noise);
672 noise_db=10*log10(Max(noise,1E-20));
673 tot_noise_db += noise_db;
677 over_noise_db += noise_db;
682 }else{ /* cod_info->block_type == SHORT_TYPE */
683 int max_index = gfc->sfb21_extra ? SBMAX_l : SBPSY_l;
685 for ( sfb = 0; sfb < max_index; sfb++ ) {
687 int s = scalefac->l[sfb];
689 if (cod_info->preflag)
692 s = cod_info->global_gain - (s << (cod_info->scalefac_scale + 1));
697 start = gfc->scalefac_band.l[ sfb ];
698 end = gfc->scalefac_band.l[ sfb+1 ];
700 for ( sum = 0.0, l = start; l < end; l++ ) {
702 temp = fabs(xr[l]) - pow43[ix[l]] * step;
705 noise = xfsf->l[sfb] = sum / l3_xmin->l[sfb];
706 max_noise=Max(max_noise,noise);
707 klemm_noise += penalties (noise);
709 noise_db=10*log10(Max(noise,1E-20));
710 /* multiplying here is adding in dB, but can overflow */
711 //tot_noise *= Max(noise, 1E-20);
712 tot_noise_db += noise_db;
716 /* multiplying here is adding in dB -but can overflow */
717 //over_noise *= noise;
718 over_noise_db += noise_db;
723 } /* cod_info->block_type == SHORT_TYPE */
725 /* normalization at this point by "count" is not necessary, since
726 * the values are only used to compare with previous values */
727 res->over_count = over;
729 /* convert to db. DO NOT CHANGE THESE */
730 /* tot_noise = is really the average over each sfb of:
731 [noise(db) - allowed_noise(db)]
733 and over_noise is the same average, only over only the
734 bands with noise > allowed_noise.
738 //res->tot_noise = 10.*log10(Max(1e-20,tot_noise ));
739 //res->over_noise = 10.*log10(Max(1e+00,over_noise));
740 res->tot_noise = tot_noise_db;
741 res->over_noise = over_noise_db;
742 res->max_noise = 10.*log10(Max(1e-20,max_noise ));
743 res->klemm_noise = 10.*log10(Max(1e-20,klemm_noise));
761 /************************************************************************
765 * updates plotting data
767 * Mark Taylor 2000-??-??
769 * Robert Hegemann: moved noise/distortion calc into it
771 ************************************************************************/
775 lame_global_flags *gfp,
776 const gr_info * const cod_info,
777 const III_psy_ratio * const ratio,
778 const III_scalefac_t * const scalefac,
779 const FLOAT8 xr[576],
780 const int l3_enc[576],
784 lame_internal_flags *gfc=gfp->internal_flags;
786 int j,i,l,start,end,bw;
788 FLOAT ifqstep = ( cod_info->scalefac_scale == 0 ) ? .5 : 1.0;
791 III_psy_xmin l3_xmin;
792 calc_noise_result noise;
795 calc_xmin (gfp,xr, ratio, cod_info, &l3_xmin);
797 calc_noise (gfc, xr, l3_enc, cod_info, &l3_xmin, scalefac, &xfsf, &noise);
799 if (cod_info->block_type == SHORT_TYPE) {
800 for (j=0, sfb = 0; sfb < SBMAX_s; sfb++ ) {
801 start = gfc->scalefac_band.s[ sfb ];
802 end = gfc->scalefac_band.s[ sfb + 1 ];
804 for ( i = 0; i < 3; i++ ) {
805 for ( en0 = 0.0, l = start; l < end; l++ ) {
806 en0 += xr[j] * xr[j];
809 en0=Max(en0/bw,1e-20);
821 tot2 += ratio->en.s[sfb][i];
823 DEBUGF("%i %i sfb=%i mdct=%f fft=%f fft-mdct=%f db \n",
825 10*log10(Max(1e-25,ratio->en.s[sfb][i])),
826 10*log10(Max(1e-25,en0)),
827 10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.s[sfb][i]))));
829 if (sfb==SBMAX_s-2) {
830 DEBUGF("%i %i toti %f %f ratio=%f db \n",gr,ch,
831 10*log10(Max(1e-25,tot2)),
832 10*log(Max(1e-25,tot1)),
833 10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
838 masking: multiplied by en0/fft_energy
839 average seems to be about -135db.
845 /* convert to MDCT units */
846 en1=1e15; /* scaling so it shows up on FFT plot */
847 gfc->pinfo->xfsf_s[gr][ch][3*sfb+i]
848 = en1*xfsf.s[sfb][i]*l3_xmin.s[sfb][i]/bw;
849 gfc->pinfo->en_s[gr][ch][3*sfb+i] = en1*en0;
851 if (ratio->en.s[sfb][i]>0)
852 en0 = en0/ratio->en.s[sfb][i];
855 if (gfp->ATHonly || gfp->ATHshort)
858 gfc->pinfo->thr_s[gr][ch][3*sfb+i] =
859 en1*Max(en0*ratio->thm.s[sfb][i],gfc->ATH->s[sfb]);
862 /* there is no scalefactor bands >= SBPSY_s */
864 gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=
865 -ifqstep*scalefac->s[sfb][i];
867 gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=0;
869 gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i] -=
870 2*cod_info->subblock_gain[i];
874 for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
875 start = gfc->scalefac_band.l[ sfb ];
876 end = gfc->scalefac_band.l[ sfb+1 ];
878 for ( en0 = 0.0, l = start; l < end; l++ )
879 en0 += xr[l] * xr[l];
880 if (!gfc->nsPsy.use) en0/=bw;
882 DEBUGF("diff = %f \n",10*log10(Max(ratio[gr][ch].en.l[sfb],1e-20))
883 -(10*log10(en0)+150));
894 tot2 += ratio->en.l[sfb];
897 DEBUGF("%i %i sfb=%i mdct=%f fft=%f fft-mdct=%f db \n",
899 10*log10(Max(1e-25,ratio->en.l[sfb])),
900 10*log10(Max(1e-25,en0)),
901 10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.l[sfb]))));
903 if (sfb==SBMAX_l-1) {
904 DEBUGF("%i %i toti %f %f ratio=%f db \n",
906 10*log10(Max(1e-25,tot2)),
907 10*log(Max(1e-25,tot1)),
908 10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
911 masking: multiplied by en0/fft_energy
912 average seems to be about -147db.
918 /* convert to MDCT units */
919 en1=1e15; /* scaling so it shows up on FFT plot */
920 gfc->pinfo->xfsf[gr][ch][sfb] = en1*xfsf.l[sfb]*l3_xmin.l[sfb]/bw;
921 gfc->pinfo->en[gr][ch][sfb] = en1*en0;
922 if (ratio->en.l[sfb]>0)
923 en0 = en0/ratio->en.l[sfb];
928 gfc->pinfo->thr[gr][ch][sfb] =
929 en1*Max(en0*ratio->thm.l[sfb],gfc->ATH->l[sfb]);
931 /* there is no scalefactor bands >= SBPSY_l */
933 if (scalefac->l[sfb]<0) /* scfsi! */
934 gfc->pinfo->LAMEsfb[gr][ch][sfb] =
935 gfc->pinfo->LAMEsfb[0][ch][sfb];
937 gfc->pinfo->LAMEsfb[gr][ch][sfb] = -ifqstep*scalefac->l[sfb];
939 gfc->pinfo->LAMEsfb[gr][ch][sfb] = 0;
942 if (cod_info->preflag && sfb>=11)
943 gfc->pinfo->LAMEsfb[gr][ch][sfb] -= ifqstep*pretab[sfb];
945 } /* block type long */
946 gfc->pinfo->LAMEqss [gr][ch] = cod_info->global_gain;
947 gfc->pinfo->LAMEmainbits[gr][ch] = cod_info->part2_3_length;
948 gfc->pinfo->LAMEsfbits [gr][ch] = cod_info->part2_length;
950 gfc->pinfo->over [gr][ch] = noise.over_count;
951 gfc->pinfo->max_noise [gr][ch] = noise.max_noise;
952 gfc->pinfo->over_noise[gr][ch] = noise.over_noise;
953 gfc->pinfo->tot_noise [gr][ch] = noise.tot_noise;
957 /************************************************************************
961 * updates plotting data for a whole frame
963 * Robert Hegemann 2000-10-21
965 ************************************************************************/
967 void set_frame_pinfo(
968 lame_global_flags *gfp,
969 FLOAT8 xr [2][2][576],
970 III_psy_ratio ratio [2][2],
971 int l3_enc [2][2][576],
972 III_scalefac_t scalefac [2][2] )
974 lame_internal_flags *gfc=gfp->internal_flags;
979 III_scalefac_t act_scalefac [2];
980 int scsfi[2] = {0,0};
983 gfc->masking_lower = 1.0;
985 /* reconstruct the scalefactors in case SCSFI was used
987 for (ch = 0; ch < gfc->channels_out; ch ++) {
988 for (sfb = 0; sfb < SBMAX_l; sfb ++) {
989 if (scalefac[1][ch].l[sfb] == -1) {/* scfsi */
990 act_scalefac[ch].l[sfb] = scalefac[0][ch].l[sfb];
993 act_scalefac[ch].l[sfb] = scalefac[1][ch].l[sfb];
998 /* for every granule and channel patch l3_enc and set info
1000 for (gr = 0; gr < gfc->mode_gr; gr ++) {
1001 for (ch = 0; ch < gfc->channels_out; ch ++) {
1003 gr_info *cod_info = &gfc->l3_side.gr[gr].ch[ch].tt;
1005 /* revert back the sign of l3enc */
1006 for ( i = 0; i < 576; i++) {
1007 if (xr[gr][ch][i] < 0)
1008 act_l3enc[i] = -l3_enc[gr][ch][i];
1010 act_l3enc[i] = +l3_enc[gr][ch][i];
1012 if (gr == 1 && scsfi[ch])
1013 set_pinfo (gfp, cod_info, &ratio[gr][ch], &act_scalefac[ch],
1014 xr[gr][ch], act_l3enc, gr, ch);
1016 set_pinfo (gfp, cod_info, &ratio[gr][ch], &scalefac[gr][ch],
1017 xr[gr][ch], act_l3enc, gr, ch);
1025 /*********************************************************************
1026 * nonlinear quantization of xr
1027 * More accurate formula than the ISO formula. Takes into account
1028 * the fact that we are quantizing xr -> ix, but we want ix^4/3 to be
1029 * as close as possible to x^4/3. (taking the nearest int would mean
1030 * ix is as close as possible to xr, which is different.)
1031 * From Segher Boessenkool <segher@eastsite.nl> 11/1999
1032 * ASM optimization from
1033 * Mathew Hendry <scampi@dial.pipex.com> 11/1999
1034 * Acy Stapp <AStapp@austin.rr.com> 11/1999
1035 * Takehiro Tominaga <tominaga@isoternet.org> 11/1999
1036 * 9/00: ASM code removed in favor of IEEE754 hack. If you need
1037 * the ASM code, check CVS circa Aug 2000.
1038 *********************************************************************/
1041 #ifdef TAKEHIRO_IEEE754_HACK
1048 #define MAGIC_FLOAT (65536*(128))
1049 #define MAGIC_INT 0x4b000000
1051 void quantize_xrpow(const FLOAT8 *xp, int *pi, FLOAT8 istep)
1053 /* quantize on xr^(3/4) instead of xr */
1057 fi = (fi_union *)pi;
1058 for (j = 576 / 4 - 1; j >= 0; --j) {
1059 double x0 = istep * xp[0];
1060 double x1 = istep * xp[1];
1061 double x2 = istep * xp[2];
1062 double x3 = istep * xp[3];
1064 x0 += MAGIC_FLOAT; fi[0].f = x0;
1065 x1 += MAGIC_FLOAT; fi[1].f = x1;
1066 x2 += MAGIC_FLOAT; fi[2].f = x2;
1067 x3 += MAGIC_FLOAT; fi[3].f = x3;
1069 fi[0].f = x0 + (adj43asm - MAGIC_INT)[fi[0].i];
1070 fi[1].f = x1 + (adj43asm - MAGIC_INT)[fi[1].i];
1071 fi[2].f = x2 + (adj43asm - MAGIC_INT)[fi[2].i];
1072 fi[3].f = x3 + (adj43asm - MAGIC_INT)[fi[3].i];
1074 fi[0].i -= MAGIC_INT;
1075 fi[1].i -= MAGIC_INT;
1076 fi[2].i -= MAGIC_INT;
1077 fi[3].i -= MAGIC_INT;
1083 # define ROUNDFAC -0.0946
1084 void quantize_xrpow_ISO(const FLOAT8 *xp, int *pi, FLOAT8 istep)
1086 /* quantize on xr^(3/4) instead of xr */
1090 fi = (fi_union *)pi;
1091 for (j=576/4 - 1;j>=0;j--) {
1092 fi[0].f = istep * xp[0] + (ROUNDFAC + MAGIC_FLOAT);
1093 fi[1].f = istep * xp[1] + (ROUNDFAC + MAGIC_FLOAT);
1094 fi[2].f = istep * xp[2] + (ROUNDFAC + MAGIC_FLOAT);
1095 fi[3].f = istep * xp[3] + (ROUNDFAC + MAGIC_FLOAT);
1097 fi[0].i -= MAGIC_INT;
1098 fi[1].i -= MAGIC_INT;
1099 fi[2].i -= MAGIC_INT;
1100 fi[3].i -= MAGIC_INT;
1108 /*********************************************************************
1109 * XRPOW_FTOI is a macro to convert floats to ints.
1110 * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x]
1113 * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]
1116 * Note: using floor() or (int) is extermely slow. On machines where
1117 * the TAKEHIRO_IEEE754_HACK code above does not work, it is worthwile
1118 * to write some ASM for XRPOW_FTOI().
1119 *********************************************************************/
1120 #define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
1121 #define QUANTFAC(rx) adj43[rx]
1122 #define ROUNDFAC 0.4054
1125 void quantize_xrpow(const FLOAT8 *xr, int *ix, FLOAT8 istep) {
1126 /* quantize on xr^(3/4) instead of xr */
1127 /* from Wilfried.Behne@t-online.de. Reported to be 2x faster than
1128 the above code (when not using ASM) on PowerPC */
1131 for ( j = 576/8; j > 0; --j) {
1132 FLOAT8 x1, x2, x3, x4, x5, x6, x7, x8;
1133 int rx1, rx2, rx3, rx4, rx5, rx6, rx7, rx8;
1136 XRPOW_FTOI(x1, rx1);
1138 XRPOW_FTOI(x2, rx2);
1140 XRPOW_FTOI(x3, rx3);
1142 XRPOW_FTOI(x4, rx4);
1144 XRPOW_FTOI(x5, rx5);
1146 XRPOW_FTOI(x6, rx6);
1148 XRPOW_FTOI(x7, rx7);
1149 x1 += QUANTFAC(rx1);
1150 XRPOW_FTOI(x8, rx8);
1151 x2 += QUANTFAC(rx2);
1152 XRPOW_FTOI(x1,*ix++);
1153 x3 += QUANTFAC(rx3);
1154 XRPOW_FTOI(x2,*ix++);
1155 x4 += QUANTFAC(rx4);
1156 XRPOW_FTOI(x3,*ix++);
1157 x5 += QUANTFAC(rx5);
1158 XRPOW_FTOI(x4,*ix++);
1159 x6 += QUANTFAC(rx6);
1160 XRPOW_FTOI(x5,*ix++);
1161 x7 += QUANTFAC(rx7);
1162 XRPOW_FTOI(x6,*ix++);
1163 x8 += QUANTFAC(rx8);
1164 XRPOW_FTOI(x7,*ix++);
1165 XRPOW_FTOI(x8,*ix++);
1174 void quantize_xrpow_ISO( const FLOAT8 *xr, int *ix, FLOAT8 istep )
1176 /* quantize on xr^(3/4) instead of xr */
1177 const FLOAT8 compareval0 = (1.0 - 0.4054)/istep;
1179 /* depending on architecture, it may be worth calculating a few more
1182 eg. compareval1 = (2.0 - 0.4054/istep);
1183 .. and then after the first compare do this ...
1184 if compareval1>*xr then ix = 1;
1186 On a pentium166, it's only worth doing the one compare (as done here),
1187 as the second compare becomes more expensive than just calculating
1188 the value. Architectures with slow FP operations may want to add some
1189 more comparevals. try it and send your diffs statistically speaking
1191 73% of all xr*istep values give ix=0
1195 for (j=576;j>0;j--) {
1196 if (compareval0 > *xr) {
1200 /* *(ix++) = (int)( istep*(*(xr++)) + 0.4054); */
1201 XRPOW_FTOI( istep*(*(xr++)) + ROUNDFAC , *(ix++) );