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
iteration_init(lame_global_flags * gfp)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
ATHmdct(lame_global_flags * gfp,FLOAT8 f)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
compute_ath(lame_global_flags * gfp,FLOAT8 ATH_l[],FLOAT8 ATH_s[])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 */
ms_convert(FLOAT8 dst[2][576],FLOAT8 src[2][576])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 ************************************************************************/
on_pe(lame_global_flags * gfp,FLOAT8 pe[2][2],III_side_info_t * l3_side,int targ_bits[2],int mean_bits,int gr)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
reduce_side(int targ_bits[2],FLOAT8 ms_ener_ratio,int mean_bits,int max_bits)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 */
calc_xmin(lame_global_flags * gfp,const FLOAT8 xr[576],const III_psy_ratio * const ratio,const gr_info * const cod_info,III_psy_xmin * const l3_xmin)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
penalties(double noise)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
calc_noise(const lame_internal_flags * const gfc,const FLOAT8 xr[576],const int ix[576],const gr_info * const cod_info,const III_psy_xmin * const l3_xmin,const III_scalefac_t * const scalefac,III_psy_xmin * xfsf,calc_noise_result * const res)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
set_pinfo(lame_global_flags * gfp,const gr_info * const cod_info,const III_psy_ratio * const ratio,const III_scalefac_t * const scalefac,const FLOAT8 xr[576],const int l3_enc[576],const int gr,const int ch)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
set_frame_pinfo(lame_global_flags * gfp,FLOAT8 xr[2][2][576],III_psy_ratio ratio[2][2],int l3_enc[2][2][576],III_scalefac_t scalefac[2][2])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
quantize_xrpow(const FLOAT8 xp[576],int pi[576],FLOAT8 istep)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
quantize_xrpow_ISO(const FLOAT8 xp[576],int pi[576],FLOAT8 istep)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
quantize_xrpow(const FLOAT8 xr[576],int ix[576],FLOAT8 istep)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
quantize_xrpow_ISO(const FLOAT8 xr[576],int ix[576],FLOAT8 istep)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
1172