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