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