xref: /plan9/sys/src/games/mp3enc/vbrquantize.c (revision 8f5875f3e9b20916b4c52ad4336922bc8653eb7b)
1 /*
2  *	MP3 quantization
3  *
4  *	Copyright (c) 1999 Mark Taylor
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: vbrquantize.c,v 1.41 2001/03/12 07:26:08 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 "l3side.h"
31 #include "quantize.h"
32 #include "reservoir.h"
33 #include "quantize_pvt.h"
34 #include "lame-analysis.h"
35 
36 #ifdef WITH_DMALLOC
37 #include <dmalloc.h>
38 #endif
39 
40 #undef MAXQUANTERROR
41 
42 
43 typedef union {
44     float f;
45     int   i;
46 } fi_union;
47 
48 #define MAGIC_FLOAT (65536*(128))
49 #define MAGIC_INT    0x4b000000
50 
51 #ifdef TAKEHIRO_IEEE754_HACK
52 
53 #ifdef MAXQUANTERROR
54 #define DUFFBLOCK() do { \
55         xp = xr34[0] * sfpow34_p1; \
56         xe = xr34[0] * sfpow34_eq; \
57         xm = xr34[0] * sfpow34_m1; \
58         if (xm > IXMAX_VAL)  \
59             return -1; \
60         xp += MAGIC_FLOAT; \
61         xe += MAGIC_FLOAT; \
62         xm += MAGIC_FLOAT; \
63         fi[0].f = xp; \
64         fi[1].f = xe; \
65         fi[2].f = xm; \
66         fi[0].f = xp + (adj43asm - MAGIC_INT)[fi[0].i]; \
67         fi[1].f = xe + (adj43asm - MAGIC_INT)[fi[1].i]; \
68         fi[2].f = xm + (adj43asm - MAGIC_INT)[fi[2].i]; \
69         fi[0].i -= MAGIC_INT; \
70         fi[1].i -= MAGIC_INT; \
71         fi[2].i -= MAGIC_INT; \
72         x0 = fabs(xr[0]); \
73         xp = x0 - pow43[fi[0].i] * sfpow_p1; \
74         xe = x0 - pow43[fi[1].i] * sfpow_eq; \
75         xm = x0 - pow43[fi[2].i] * sfpow_m1; \
76         xp *= xp; \
77         xe *= xe; \
78         xm *= xm; \
79         xfsf_eq = Max(xfsf_eq, xe); \
80         xfsf_p1 = Max(xfsf_p1, xp); \
81         xfsf_m1 = Max(xfsf_m1, xm); \
82         ++xr; \
83         ++xr34; \
84     } while(0)
85 #else
86 #define DUFFBLOCK() do { \
87         xp = xr34[0] * sfpow34_p1; \
88         xe = xr34[0] * sfpow34_eq; \
89         xm = xr34[0] * sfpow34_m1; \
90         if (xm > IXMAX_VAL)  \
91             return -1; \
92         xp += MAGIC_FLOAT; \
93         xe += MAGIC_FLOAT; \
94         xm += MAGIC_FLOAT; \
95         fi[0].f = xp; \
96         fi[1].f = xe; \
97         fi[2].f = xm; \
98         fi[0].f = xp + (adj43asm - MAGIC_INT)[fi[0].i]; \
99         fi[1].f = xe + (adj43asm - MAGIC_INT)[fi[1].i]; \
100         fi[2].f = xm + (adj43asm - MAGIC_INT)[fi[2].i]; \
101         fi[0].i -= MAGIC_INT; \
102         fi[1].i -= MAGIC_INT; \
103         fi[2].i -= MAGIC_INT; \
104         x0 = fabs(xr[0]); \
105         xp = x0 - pow43[fi[0].i] * sfpow_p1; \
106         xe = x0 - pow43[fi[1].i] * sfpow_eq; \
107         xm = x0 - pow43[fi[2].i] * sfpow_m1; \
108         xfsf_p1 += xp * xp; \
109         xfsf_eq += xe * xe; \
110         xfsf_m1 += xm * xm; \
111         ++xr; \
112         ++xr34; \
113     } while(0)
114 #endif
115 
116 #else
117 
118 /*********************************************************************
119  * XRPOW_FTOI is a macro to convert floats to ints.
120  * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x]
121  *                                         ROUNDFAC= -0.0946
122  *
123  * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]
124  *                                   ROUNDFAC=0.4054
125  *********************************************************************/
126 #  define QUANTFAC(rx)  adj43[rx]
127 #  define ROUNDFAC 0.4054
128 #  define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
129 
130 
131 #endif
132 
133 
134 static FLOAT8
calc_sfb_noise(const FLOAT8 * xr,const FLOAT8 * xr34,const int bw,const int sf)135 calc_sfb_noise(const FLOAT8 *xr, const FLOAT8 *xr34, const int bw, const int sf)
136 {
137   int j;
138   fi_union fi;
139   FLOAT8 temp;
140   FLOAT8 xfsf=0;
141   FLOAT8 sfpow,sfpow34;
142 
143   sfpow = POW20(sf+210); /*pow(2.0,sf/4.0); */
144   sfpow34  = IPOW20(sf+210); /*pow(sfpow,-3.0/4.0);*/
145 
146   for ( j=0; j < bw ; ++j) {
147 #if 0
148     int ix;
149     if (xr34[j]*sfpow34 > IXMAX_VAL) return -1;
150     ix=floor( xr34[j]*sfpow34);
151     temp = fabs(xr[j])- pow43[ix]*sfpow;
152     temp *= temp;
153 
154     if (ix < IXMAX_VAL) {
155       temp2 = fabs(xr[j])- pow43[ix+1]*sfpow;
156       temp2 *=temp2;
157       if (temp2<temp) {
158 	temp=temp2;
159 	++ix;
160       }
161     }
162 #else
163     if (xr34[j]*sfpow34 > IXMAX_VAL) return -1;
164 
165 #ifdef TAKEHIRO_IEEE754_HACK
166     temp   = sfpow34*xr34[j];
167     temp  += MAGIC_FLOAT;
168     fi.f  = temp;
169     fi.f  = temp + (adj43asm - MAGIC_INT)[fi.i];
170     fi.i -= MAGIC_INT;
171 #else
172     temp = xr34[j]*sfpow34;
173     XRPOW_FTOI(temp, fi.i);
174     XRPOW_FTOI(temp + QUANTFAC(fi.i), fi.i);
175 #endif
176 
177     temp = fabs(xr[j])- pow43[fi.i]*sfpow;
178     temp *= temp;
179 
180 #endif
181 
182 #ifdef MAXQUANTERROR
183     xfsf = Max(xfsf,temp);
184 #else
185     xfsf += temp;
186 #endif
187   }
188 #ifdef MAXQUANTERROR
189   return xfsf;
190 #else
191   return xfsf;//bw;
192 #endif
193 }
194 
195 
196 
197 
198 static FLOAT8
calc_sfb_noise_ave(const FLOAT8 * xr,const FLOAT8 * xr34,const int bw,const int sf)199 calc_sfb_noise_ave(const FLOAT8 *xr, const FLOAT8 *xr34, const int bw, const int sf)
200 {
201     double xp;
202     double xe;
203     double xm;
204 #ifdef TAKEHIRO_IEEE754_HACK
205     double x0;
206 #endif
207     int xx[3], j;
208     fi_union *fi = (fi_union *)xx;
209     FLOAT8 sfpow34_eq, sfpow34_p1, sfpow34_m1;
210     FLOAT8 sfpow_eq, sfpow_p1, sfpow_m1;
211     FLOAT8 xfsf_eq = 0, xfsf_p1 = 0, xfsf_m1 = 0;
212 
213     sfpow_eq = POW20(sf + 210); /*pow(2.0,sf/4.0); */
214     sfpow_m1 = sfpow_eq * .8408964153;  /* pow(2,(sf-1)/4.0) */
215     sfpow_p1 = sfpow_eq * 1.189207115;
216 
217     sfpow34_eq = IPOW20(sf + 210); /*pow(sfpow,-3.0/4.0);*/
218     sfpow34_m1 = sfpow34_eq * 1.13878863476;       /* .84089 ^ -3/4 */
219     sfpow34_p1 = sfpow34_eq * 0.878126080187;
220 
221 #ifdef TAKEHIRO_IEEE754_HACK
222     /*
223      *  loop unrolled into "Duff's Device".   Robert Hegemann
224      */
225     j = (bw+3) / 4;
226     switch (bw % 4) {
227         default:
228         case 0: do{ DUFFBLOCK();
229         case 3:     DUFFBLOCK();
230         case 2:     DUFFBLOCK();
231         case 1:     DUFFBLOCK(); } while (--j);
232     }
233 #else
234     for (j = 0; j < bw; ++j) {
235 
236         if (xr34[j]*sfpow34_m1 > IXMAX_VAL) return -1;
237 
238         xe = xr34[j]*sfpow34_eq;
239         XRPOW_FTOI(xe, fi[0].i);
240         XRPOW_FTOI(xe + QUANTFAC(fi[0].i), fi[0].i);
241         xe = fabs(xr[j])- pow43[fi[0].i]*sfpow_eq;
242         xe *= xe;
243 
244         xp = xr34[j]*sfpow34_p1;
245         XRPOW_FTOI(xp, fi[0].i);
246         XRPOW_FTOI(xp + QUANTFAC(fi[0].i), fi[0].i);
247         xp = fabs(xr[j])- pow43[fi[0].i]*sfpow_p1;
248         xp *= xp;
249 
250         xm = xr34[j]*sfpow34_m1;
251         XRPOW_FTOI(xm, fi[0].i);
252         XRPOW_FTOI(xm + QUANTFAC(fi[0].i), fi[0].i);
253         xm = fabs(xr[j])- pow43[fi[0].i]*sfpow_m1;
254         xm *= xm;
255 
256 #ifdef MAXQUANTERROR
257         xfsf_eq = Max(xfsf,xm);
258         xfsf_p1 = Max(xfsf_p1,xp);
259         xfsf_m1 = Max(xfsf_m1,xm);
260 #else
261         xfsf_eq += xe;
262         xfsf_p1 += xp;
263         xfsf_m1 += xm;
264 #endif
265     }
266 #endif
267 
268     if (xfsf_eq < xfsf_p1)
269         xfsf_eq = xfsf_p1;
270     if (xfsf_eq < xfsf_m1)
271         xfsf_eq = xfsf_m1;
272 #ifdef MAXQUANTERROR
273     return xfsf_eq;
274 #else
275     return xfsf_eq;//bw;
276 #endif
277 }
278 
279 
280 
281 static int
find_scalefac(const FLOAT8 * xr,const FLOAT8 * xr34,const int sfb,const FLOAT8 l3_xmin,const int bw)282 find_scalefac(const FLOAT8 *xr, const FLOAT8 *xr34, const int sfb,
283 		     const FLOAT8 l3_xmin, const int bw)
284 {
285   FLOAT8 xfsf;
286   int i,sf,sf_ok,delsf;
287 
288   /* search will range from sf:  -209 -> 45  */
289   sf = -82;
290   delsf = 128;
291 
292   sf_ok=10000;
293   for (i=0; i<7; i++) {
294     delsf /= 2;
295     xfsf = calc_sfb_noise(xr,xr34,bw,sf);
296 
297     if (xfsf < 0) {
298       /* scalefactors too small */
299       sf += delsf;
300     }else{
301       if (sf_ok==10000) sf_ok=sf;
302       if (xfsf > l3_xmin)  {
303 	/* distortion.  try a smaller scalefactor */
304 	sf -= delsf;
305       }else{
306 	sf_ok = sf;
307 	sf += delsf;
308       }
309     }
310   }
311   assert(sf_ok!=10000);
312 #if 0
313   assert(delsf==1);  /* when for loop goes up to 7 */
314 #endif
315 
316   return sf;
317 }
318 
319 static int
find_scalefac_ave(const FLOAT8 * xr,const FLOAT8 * xr34,const int sfb,const FLOAT8 l3_xmin,const int bw)320 find_scalefac_ave(const FLOAT8 *xr, const FLOAT8 *xr34, const int sfb,
321 		     const FLOAT8 l3_xmin, const int bw)
322 {
323   FLOAT8 xfsf;
324   int i,sf,sf_ok,delsf;
325 
326   /* search will range from sf:  -209 -> 45  */
327   sf = -82;
328   delsf = 128;
329 
330   sf_ok=10000;
331   for (i=0; i<7; i++) {
332     delsf /= 2;
333     xfsf = calc_sfb_noise_ave(xr,xr34,bw,sf);
334 
335     if (xfsf < 0) {
336       /* scalefactors too small */
337       sf += delsf;
338     }else{
339       if (sf_ok==10000) sf_ok=sf;
340       if (xfsf > l3_xmin)  {
341 	/* distortion.  try a smaller scalefactor */
342 	sf -= delsf;
343       }else{
344 	sf_ok = sf;
345 	sf += delsf;
346       }
347     }
348   }
349   assert(sf_ok!=10000);
350 #if 0
351   assert(delsf==1);  /* when for loop goes up to 7 */
352 #endif
353 
354   return sf;
355 }
356 
357 
358 
359 /*  ???
360         How do the following tables look like for MPEG-2-LSF
361                                                             ??? */
362 
363 
364 
365 static const int max_range_short[SBPSY_s] =
366 {15, 15, 15, 15, 15, 15, 7, 7, 7, 7, 7, 7 };
367 
368 static const int max_range_long[SBPSY_l] =
369 {15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7};
370 
371 static const int max_range_short_lsf[SBPSY_s] =
372 {15, 15, 15, 15, 15, 15, 7, 7, 7, 7, 7, 7 };
373 
374 /*static const int max_range_short_lsf_pretab[SBPSY_s] =
375 {}*/
376 
377 static const int max_range_long_lsf[SBPSY_l] =
378 {15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7};
379 
380 static const int max_range_long_lsf_pretab[SBPSY_l] =
381 { 7,7,7,7,7,7, 3,3,3,3,3, 0,0,0,0, 0,0,0, 0,0,0 };
382 
383 
384 static int
compute_scalefacs_short_lsf(int sf[SBPSY_s][3],gr_info * cod_info,int scalefac[SBPSY_s][3],int sbg[3])385 compute_scalefacs_short_lsf (
386     int sf[SBPSY_s][3],gr_info *cod_info, int scalefac[SBPSY_s][3],int sbg[3])
387 {
388     int maxrange, maxrange1, maxrange2, maxover;
389     int sfb, i;
390     int ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
391 
392     maxover   = 0;
393     maxrange1 = max_range_short_lsf[0];
394     maxrange2 = max_range_short_lsf[6];
395 
396 
397     for (i=0; i<3; ++i) {
398         int maxsf1 = 0, maxsf2 = 0, minsf = 1000;
399         /* see if we should use subblock gain */
400         for (sfb = 0; sfb < SBPSY_s; sfb++) {
401             if (sfb < 6) {
402                 if (maxsf1 < -sf[sfb][i])
403                     maxsf1 = -sf[sfb][i];
404             } else {
405                 if (maxsf2 < -sf[sfb][i])
406                     maxsf2 = -sf[sfb][i];
407             }
408             if (minsf > -sf[sfb][i])
409                 minsf = -sf[sfb][i];
410         }
411 
412         /* boost subblock gain as little as possible so we can
413          * reach maxsf1 with scalefactors
414          * 8*sbg >= maxsf1
415          */
416         maxsf1 = Max (maxsf1-maxrange1*ifqstep, maxsf2-maxrange2*ifqstep);
417         sbg[i] = 0;
418         if (minsf > 0)
419             sbg[i] = floor (.125*minsf + .001);
420         if (maxsf1 > 0)
421             sbg[i] = Max (sbg[i], (maxsf1/8 + (maxsf1 % 8 != 0)));
422 
423         if (sbg[i] > 7)
424             sbg[i] = 7;
425 
426 
427         for (sfb = 0; sfb < SBPSY_s; sfb++) {
428             sf[sfb][i] += 8*sbg[i];
429 
430             if (sf[sfb][i] < 0) {
431                 maxrange = sfb < 6 ? maxrange1 : maxrange2;
432 
433                 scalefac[sfb][i]
434                      = -sf[sfb][i]/ifqstep + (-sf[sfb][i]%ifqstep != 0);
435 
436                 if (scalefac[sfb][i] > maxrange)
437                     scalefac[sfb][i] = maxrange;
438 
439                 if (maxover < -(sf[sfb][i] + scalefac[sfb][i]*ifqstep))
440                     maxover = -(sf[sfb][i] + scalefac[sfb][i]*ifqstep);
441             }
442         }
443     }
444 
445     return maxover;
446 }
447 
448 static int
compute_scalefacs_long_lsf(int sf[SBPSY_l],const gr_info * const cod_info,int scalefac[SBPSY_l])449 compute_scalefacs_long_lsf (
450               int             sf       [SBPSY_l],
451         const gr_info * const cod_info,
452               int             scalefac [SBPSY_l] )
453 {
454     const int * max_range = max_range_long_lsf;
455     int ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
456     int sfb;
457     int maxover;
458 
459 
460     if (cod_info->preflag) {
461         max_range = max_range_long_lsf_pretab;
462         for (sfb = 11; sfb < SBPSY_l; sfb++)
463             sf[sfb] += pretab[sfb] * ifqstep;
464     }
465 
466     maxover = 0;
467     for (sfb = 0; sfb < SBPSY_l; sfb++) {
468 
469         if (sf[sfb] < 0) {
470             /* ifqstep*scalefac >= -sf[sfb], so round UP */
471             scalefac[sfb] = -sf[sfb]/ifqstep  + (-sf[sfb] % ifqstep != 0);
472             if (scalefac[sfb] > max_range[sfb])
473                 scalefac[sfb] = max_range[sfb];
474 
475             /* sf[sfb] should now be positive: */
476             if (-(sf[sfb] + scalefac[sfb]*ifqstep) > maxover) {
477                 maxover = -(sf[sfb] + scalefac[sfb]*ifqstep);
478             }
479         }
480     }
481 
482     return maxover;
483 }
484 
485 
486 /*
487     sfb=0..5  scalefac < 16
488     sfb>5     scalefac < 8
489 
490     ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
491     ol_sf =  (cod_info->global_gain-210.0);
492     ol_sf -= 8*cod_info->subblock_gain[i];
493     ol_sf -= ifqstep*scalefac[gr][ch].s[sfb][i];
494 
495 */
496 static int
compute_scalefacs_short(int sf[SBPSY_s][3],gr_info * cod_info,int scalefac[SBPSY_s][3],int sbg[3])497 compute_scalefacs_short(int sf[SBPSY_s][3],gr_info *cod_info,
498 int scalefac[SBPSY_s][3],int sbg[3])
499 {
500   int maxrange,maxrange1,maxrange2,maxover;
501   int sfb,i;
502   int ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
503 
504   maxover=0;
505   maxrange1 = 15;
506   maxrange2 = 7;
507 
508 
509   for (i=0; i<3; ++i) {
510     int maxsf1=0,maxsf2=0,minsf=1000;
511     /* see if we should use subblock gain */
512     for ( sfb = 0; sfb < SBPSY_s; sfb++ ) {
513       if (sfb < 6) {
514 	if (-sf[sfb][i]>maxsf1) maxsf1 = -sf[sfb][i];
515       } else {
516 	if (-sf[sfb][i]>maxsf2) maxsf2 = -sf[sfb][i];
517       }
518       if (-sf[sfb][i]<minsf) minsf = -sf[sfb][i];
519     }
520 
521     /* boost subblock gain as little as possible so we can
522      * reach maxsf1 with scalefactors
523      * 8*sbg >= maxsf1
524      */
525     maxsf1 = Max(maxsf1-maxrange1*ifqstep,maxsf2-maxrange2*ifqstep);
526     sbg[i]=0;
527     if (minsf >0 ) sbg[i] = floor(.125*minsf + .001);
528     if (maxsf1 > 0)  sbg[i] = Max(sbg[i],(maxsf1/8 + (maxsf1 % 8 != 0)));
529     if (sbg[i] > 7) sbg[i]=7;
530 
531 
532     for ( sfb = 0; sfb < SBPSY_s; sfb++ ) {
533       sf[sfb][i] += 8*sbg[i];
534 
535       if (sf[sfb][i] < 0) {
536 	maxrange = sfb < 6 ? maxrange1 : maxrange2;
537         scalefac[sfb][i]= -sf[sfb][i]/ifqstep + (-sf[sfb][i]%ifqstep != 0);
538 	if (scalefac[sfb][i]>maxrange) scalefac[sfb][i]=maxrange;
539 
540 	if (-(sf[sfb][i] + scalefac[sfb][i]*ifqstep) >maxover)  {
541 	  maxover=-(sf[sfb][i] + scalefac[sfb][i]*ifqstep);
542 	}
543       }
544     }
545   }
546 
547   return maxover;
548 }
549 
550 
551 
552 /*
553 	  ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
554 	  ol_sf =  (cod_info->global_gain-210.0);
555 	  ol_sf -= ifqstep*scalefac[gr][ch].l[sfb];
556 	  if (cod_info->preflag && sfb>=11)
557 	  ol_sf -= ifqstep*pretab[sfb];
558 */
559 static int
compute_scalefacs_long(int sf[SBPSY_l],gr_info * cod_info,int scalefac[SBPSY_l])560 compute_scalefacs_long(int sf[SBPSY_l],gr_info *cod_info,int scalefac[SBPSY_l])
561 {
562   int sfb;
563   int maxover;
564   int ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
565 
566 
567   if (cod_info->preflag)
568     for ( sfb = 11; sfb < SBPSY_l; sfb++ )
569       sf[sfb] += pretab[sfb]*ifqstep;
570 
571 
572   maxover=0;
573   for ( sfb = 0; sfb < SBPSY_l; sfb++ ) {
574 
575     if (sf[sfb]<0) {
576       /* ifqstep*scalefac >= -sf[sfb], so round UP */
577       scalefac[sfb]= -sf[sfb]/ifqstep + (-sf[sfb] % ifqstep != 0);
578       if (scalefac[sfb] > max_range_long[sfb]) scalefac[sfb]=max_range_long[sfb];
579 
580       /* sf[sfb] should now be positive: */
581       if (  -(sf[sfb] + scalefac[sfb]*ifqstep)  > maxover) {
582 	maxover = -(sf[sfb] + scalefac[sfb]*ifqstep);
583       }
584     }
585   }
586 
587   return maxover;
588 }
589 
590 
591 
592 
593 
594 
595 
596 
597 /************************************************************************
598  *
599  * quantize and encode with the given scalefacs and global gain
600  *
601  * compute scalefactors, l3_enc, and return number of bits needed to encode
602  *
603  *
604  ************************************************************************/
605 
606 static int
VBR_quantize_granule(lame_global_flags * gfp,FLOAT8 xr34[576],int l3_enc[576],const III_psy_ratio * const ratio,III_scalefac_t * const scalefac,const int gr,const int ch)607 VBR_quantize_granule(
608           lame_global_flags *gfp,
609           FLOAT8                 xr34[576],
610           int                    l3_enc[576],
611     const III_psy_ratio  * const ratio,
612           III_scalefac_t * const scalefac,
613     const int gr,
614     const int ch)
615 {
616   lame_internal_flags *gfc=gfp->internal_flags;
617   int status;
618   gr_info *cod_info;
619   III_side_info_t * l3_side;
620   l3_side = &gfc->l3_side;
621   cod_info = &l3_side->gr[gr].ch[ch].tt;
622 
623 
624   /* encode scalefacs */
625   if ( gfc->is_mpeg1 )
626     status=scale_bitcount(scalefac, cod_info);
627   else
628     status=scale_bitcount_lsf(gfc,scalefac, cod_info);
629 
630   if (status!=0) {
631     return -1;
632   }
633 
634   /* quantize xr34 */
635   cod_info->part2_3_length = count_bits(gfc,l3_enc,xr34,cod_info);
636   if (cod_info->part2_3_length >= LARGE_BITS) return -2;
637   cod_info->part2_3_length += cod_info->part2_length;
638 
639 
640   if (gfc->use_best_huffman==1) {
641     best_huffman_divide(gfc, gr, ch, cod_info, l3_enc);
642   }
643   return 0;
644 }
645 
646 
647 
648 /***********************************************************************
649  *
650  *      calc_short_block_vbr_sf()
651  *      calc_long_block_vbr_sf()
652  *
653  *  Mark Taylor 2000-??-??
654  *  Robert Hegemann 2000-10-25 made functions of it
655  *
656  ***********************************************************************/
657 static const int MAX_SF_DELTA = 4;
658 
659 static int
short_block_vbr_sf(const lame_internal_flags * const gfc,const III_psy_xmin * const l3_xmin,const FLOAT8 xr34_orig[576],const FLOAT8 xr34[576],III_scalefac_t * const vbrsf)660 short_block_vbr_sf (
661     const lame_internal_flags        * const gfc,
662     const III_psy_xmin   * const l3_xmin,
663     const FLOAT8                 xr34_orig[576],
664     const FLOAT8                 xr34     [576],
665           III_scalefac_t * const vbrsf )
666 {
667      int j, sfb, b;
668     int vbrmax = -10000; /* initialize for maximum search */
669 
670     for (j = 0, sfb = 0; sfb < SBMAX_s; sfb++) {
671         for (b = 0; b < 3; b++) {
672 	    const  int start = gfc->scalefac_band.s[ sfb ];
673 	    const  int end   = gfc->scalefac_band.s[ sfb+1 ];
674 	    const  int width = end - start;
675 
676             vbrsf->s[sfb][b] = find_scalefac_ave (&xr34[j], &xr34_orig[j],
677                                               sfb, l3_xmin->s[sfb][b], width);
678             j += width;
679         }
680     }
681 
682     for (sfb = 0; sfb < SBMAX_s; sfb++) {
683         for (b = 0; b < 3; b++) {
684 	    if (sfb > 0)
685 	        if (vbrsf->s[sfb][b] > vbrsf->s[sfb-1][b]+MAX_SF_DELTA)
686                     vbrsf->s[sfb][b] = vbrsf->s[sfb-1][b]+MAX_SF_DELTA;
687 	    if (sfb < SBMAX_s-1)
688 	        if (vbrsf->s[sfb][b] > vbrsf->s[sfb+1][b]+MAX_SF_DELTA)
689                     vbrsf->s[sfb][b] = vbrsf->s[sfb+1][b]+MAX_SF_DELTA;
690             if (vbrmax < vbrsf->s[sfb][b])
691                 vbrmax = vbrsf->s[sfb][b];
692         }
693     }
694 
695     return vbrmax;
696 }
697 
698 
699 
700 static int
long_block_vbr_sf(const lame_internal_flags * const gfc,const III_psy_xmin * const l3_xmin,const FLOAT8 xr34_orig[576],const FLOAT8 xr34[576],III_scalefac_t * const vbrsf)701 long_block_vbr_sf (
702     const lame_internal_flags        * const gfc,
703     const III_psy_xmin   * const l3_xmin,
704     const FLOAT8                 xr34_orig[576],
705     const FLOAT8                 xr34     [576],
706           III_scalefac_t * const vbrsf )
707 {
708      int sfb;
709     int vbrmax = -10000; /* initialize for maximum search */
710 
711     for (sfb = 0; sfb < SBMAX_l; sfb++) {
712         const  int start = gfc->scalefac_band.l[ sfb ];
713         const  int end   = gfc->scalefac_band.l[ sfb+1 ];
714         const  int width = end - start;
715 
716         vbrsf->l[sfb] = find_scalefac_ave (&xr34[start], &xr34_orig[start],
717                                                sfb, l3_xmin->l[sfb], width);
718     }
719 
720     for (sfb = 0; sfb < SBMAX_l; sfb++) {
721         if (sfb > 0)
722 	    if (vbrsf->l[sfb] > vbrsf->l[sfb-1]+MAX_SF_DELTA)
723                 vbrsf->l[sfb] = vbrsf->l[sfb-1]+MAX_SF_DELTA;
724         if (sfb < SBMAX_l-1)
725 	    if (vbrsf->l[sfb] > vbrsf->l[sfb+1]+MAX_SF_DELTA)
726                 vbrsf->l[sfb] = vbrsf->l[sfb+1]+MAX_SF_DELTA;
727         if (vbrmax < vbrsf->l[sfb])
728             vbrmax = vbrsf->l[sfb];
729     }
730 
731     return vbrmax;
732 }
733 
734 
735     /* a variation for vbr-mtrh */
736 static int
short_block_sf(const lame_internal_flags * const gfc,const III_psy_xmin * const l3_xmin,const FLOAT8 xr34_orig[576],const FLOAT8 xr34[576],III_scalefac_t * const vbrsf)737 short_block_sf (
738     const lame_internal_flags        * const gfc,
739     const III_psy_xmin   * const l3_xmin,
740     const FLOAT8                 xr34_orig[576],
741     const FLOAT8                 xr34     [576],
742           III_scalefac_t * const vbrsf )
743 {
744      int j, sfb, b;
745     int vbrmean, vbrmin, vbrmax;
746     int sf_cache[SBMAX_s];
747 
748     for (j = 0, sfb = 0; sfb < SBMAX_s; sfb++) {
749         for (b = 0; b < 3; b++) {
750 	    const  int start = gfc->scalefac_band.s[ sfb ];
751 	    const  int end   = gfc->scalefac_band.s[ sfb+1 ];
752 	    const  int width = end - start;
753 
754             if (0 == gfc->noise_shaping_amp) {
755                 /*  the faster and sloppier mode to use at lower quality
756                  */
757                 vbrsf->s[sfb][b] = find_scalefac (&xr34[j], &xr34_orig[j], sfb,
758                                               l3_xmin->s[sfb][b], width);
759 	    }
760             else {
761                 /*  the slower and better mode to use at higher quality
762                  */
763                 vbrsf->s[sfb][b] = find_scalefac_ave (&xr34[j], &xr34_orig[j],
764                                               sfb, l3_xmin->s[sfb][b], width);
765             }
766             j += width;
767         }
768     }
769 
770     vbrmax = -10000;
771 
772     for (b = 0; b < 3; b++) {
773 
774         /*  make working copy, select_kth_int will reorder!
775          */
776         for (sfb = 0; sfb < SBMAX_s; sfb++)
777             sf_cache[sfb] = vbrsf->s[sfb][b];
778 
779         /*  find median value, take it as mean
780          */
781         vbrmean = select_kth_int (sf_cache, SBMAX_s, (SBMAX_s+1)/2);
782 
783         /*  get min value
784          */
785         vbrmin = 10000;
786         for (sfb = 0; sfb < SBMAX_s; sfb++) {
787             if (vbrmin > vbrsf->s[sfb][b])
788                 vbrmin = vbrsf->s[sfb][b];
789         }
790 
791         /*  patch sfb12
792          */
793         vbrsf->s[SBPSY_s][b] = Min (vbrsf->s[SBPSY_s][b], vbrmean);
794         vbrsf->s[SBPSY_s][b] = Max (vbrsf->s[SBPSY_s][b], vbrmin-(vbrmean-vbrmin));
795 
796         /*  cut peaks
797          */
798         for (sfb = 0; sfb < SBMAX_s; sfb++) {
799             if (vbrsf->s[sfb][b] > vbrmean+(vbrmean-vbrmin))
800                 vbrsf->s[sfb][b] = vbrmean+(vbrmean-vbrmin);
801         }
802 
803         /*  get max value
804          */
805         for (sfb = 0; sfb < SBMAX_s; sfb++) {
806             if (vbrmax < vbrsf->s[sfb][b])
807                 vbrmax = vbrsf->s[sfb][b];
808         }
809     }
810 
811     return vbrmax;
812 }
813 
814 
815     /* a variation for vbr-mtrh */
816 static int
long_block_sf(const lame_internal_flags * const gfc,const III_psy_xmin * const l3_xmin,const FLOAT8 xr34_orig[576],const FLOAT8 xr34[576],III_scalefac_t * const vbrsf)817 long_block_sf (
818     const lame_internal_flags        * const gfc,
819     const III_psy_xmin   * const l3_xmin,
820     const FLOAT8                 xr34_orig[576],
821     const FLOAT8                 xr34     [576],
822           III_scalefac_t * const vbrsf )
823 {
824      int sfb;
825     int vbrmean, vbrmin, vbrmax;
826     int sf_cache[SBMAX_l];
827 
828     for (sfb = 0; sfb < SBMAX_l; sfb++) {
829         const  int start = gfc->scalefac_band.l[ sfb ];
830         const  int end   = gfc->scalefac_band.l[ sfb+1 ];
831         const  int width = end - start;
832 
833         if (0 == gfc->noise_shaping_amp) {
834             /*  the faster and sloppier mode to use at lower quality
835              */
836             vbrsf->l[sfb] = find_scalefac (&xr34[start], &xr34_orig[start],
837                                            sfb, l3_xmin->l[sfb], width);
838         }
839         else {
840             /*  the slower and better mode to use at higher quality
841              */
842             vbrsf->l[sfb] = find_scalefac_ave (&xr34[start], &xr34_orig[start],
843                                                sfb, l3_xmin->l[sfb], width);
844         }
845     }
846 
847     /*  make working copy, select_kth_int will reorder!
848      */
849     for (sfb = 0; sfb < SBMAX_l; sfb++)
850         sf_cache[sfb] = vbrsf->l[sfb];
851 
852     /*  find median value, take it as mean
853      */
854     vbrmean = select_kth_int (sf_cache, SBMAX_l, (SBMAX_l+1)/2);
855 
856     /*  get min value
857      */
858     vbrmin = +10000;
859     for (sfb = 0; sfb < SBMAX_l; sfb++) {
860         if (vbrmin > vbrsf->l[sfb])
861             vbrmin = vbrsf->l[sfb];
862     }
863 
864     /*  patch sfb21
865      */
866     vbrsf->l[SBPSY_l] = Min (vbrsf->l[SBPSY_l], vbrmean);
867     vbrsf->l[SBPSY_l] = Max (vbrsf->l[SBPSY_l], vbrmin-(vbrmean-vbrmin));
868 
869     /*  cut peaks
870      */
871     for (sfb = 0; sfb < SBMAX_l; sfb++) {
872         if (vbrsf->l[sfb] > vbrmean+(vbrmean-vbrmin))
873             vbrsf->l[sfb] = vbrmean+(vbrmean-vbrmin);
874     }
875 
876     /*  get max value
877      */
878     vbrmax = -10000;
879     for (sfb = 0; sfb < SBMAX_l; sfb++) {
880         if (vbrmax < vbrsf->l[sfb])
881             vbrmax = vbrsf->l[sfb];
882     }
883 
884     return vbrmax;
885 }
886 
887 
888 
889 /******************************************************************
890  *
891  *  short block scalefacs
892  *
893  ******************************************************************/
894 
895 static void
short_block_scalefacs(lame_global_flags * gfp,gr_info * const cod_info,III_scalefac_t * const scalefac,III_scalefac_t * const vbrsf,int * const VBRmax)896 short_block_scalefacs (
897        lame_global_flags *gfp,
898           gr_info        * const cod_info,
899           III_scalefac_t * const scalefac,
900           III_scalefac_t * const vbrsf,
901           int            * const VBRmax )
902 {
903     lame_internal_flags *gfc=gfp->internal_flags;
904     const int * max_range;
905      int sfb, b;
906     int maxover, maxover0, maxover1, mover;
907     int v0, v1;
908     int minsfb;
909     int vbrmax = *VBRmax;
910 
911     max_range = gfc->is_mpeg1 ? max_range_short : max_range_short_lsf;
912 
913     maxover0 = 0;
914     maxover1 = 0;
915     for (sfb = 0; sfb < SBPSY_s; sfb++) {
916         for (b = 0; b < 3; b++) {
917             v0 = (vbrmax - vbrsf->s[sfb][b]) - (4*14 + 2*max_range[sfb]);
918             v1 = (vbrmax - vbrsf->s[sfb][b]) - (4*14 + 4*max_range[sfb]);
919             if (maxover0 < v0)
920                 maxover0 = v0;
921             if (maxover1 < v1)
922                 maxover1 = v1;
923         }
924     }
925 
926     if (gfc->noise_shaping == 2)
927         /* allow scalefac_scale=1 */
928         mover = Min (maxover0, maxover1);
929     else
930         mover = maxover0;
931 
932     vbrmax   -= mover;
933     maxover0 -= mover;
934     maxover1 -= mover;
935 
936     if (maxover0 == 0)
937         cod_info->scalefac_scale = 0;
938     else if (maxover1 == 0)
939         cod_info->scalefac_scale = 1;
940 
941     /* sf =  (cod_info->global_gain-210.0) */
942     cod_info->global_gain = vbrmax + 210;
943     assert(cod_info->global_gain < 256);
944 
945     if (vbr_mtrh == gfp->VBR && cod_info->global_gain > 1) {
946         /*  just to be safe, reduce global_gain by one
947          */
948         cod_info->global_gain -= 1;
949     }
950 
951     if (cod_info->global_gain > 255)
952         cod_info->global_gain = 255;
953 
954     for (sfb = 0; sfb < SBPSY_s; sfb++) {
955         for (b = 0; b < 3; b++) {
956             vbrsf->s[sfb][b] -= vbrmax;
957         }
958     }
959     if ( gfc->is_mpeg1 )
960         maxover = compute_scalefacs_short (vbrsf->s, cod_info, scalefac->s,
961                                            cod_info->subblock_gain);
962     else
963         maxover = compute_scalefacs_short_lsf (vbrsf->s, cod_info, scalefac->s,
964                                                cod_info->subblock_gain);
965 
966     assert (maxover <= 0);
967 
968     /* adjust global_gain so at least 1 subblock gain = 0 */
969     minsfb = 999; /* prepare for minimum search */
970     for (b = 0; b < 3; b++)
971         if (minsfb > cod_info->subblock_gain[b])
972             minsfb = cod_info->subblock_gain[b];
973 
974     if (minsfb > cod_info->global_gain/8)
975         minsfb = cod_info->global_gain/8;
976 
977     vbrmax                -= 8*minsfb;
978     cod_info->global_gain -= 8*minsfb;
979 
980     for (b = 0; b < 3; b++)
981         cod_info->subblock_gain[b] -= minsfb;
982 
983     *VBRmax = vbrmax;
984 }
985 
986 
987 
988 /******************************************************************
989  *
990  *  long block scalefacs
991  *
992  ******************************************************************/
993 
994 static void
long_block_scalefacs(lame_global_flags * gfp,gr_info * const cod_info,III_scalefac_t * const scalefac,III_scalefac_t * const vbrsf,int * const VBRmax)995 long_block_scalefacs (
996     lame_global_flags *gfp,
997           gr_info        * const cod_info,
998           III_scalefac_t * const scalefac,
999           III_scalefac_t * const vbrsf,
1000           int            * const VBRmax )
1001 {
1002     lame_internal_flags *gfc=gfp->internal_flags;
1003     const int * max_range;
1004     const int * max_rangep;
1005      int sfb;
1006     int maxover, maxover0, maxover1, maxover0p, maxover1p, mover;
1007     int v0, v1, v0p, v1p;
1008     int vbrmax = *VBRmax;
1009 
1010     max_range  = gfc->is_mpeg1 ? max_range_long : max_range_long_lsf;
1011     max_rangep = gfc->is_mpeg1 ? max_range_long : max_range_long_lsf_pretab;
1012 
1013     maxover0  = 0;
1014     maxover1  = 0;
1015     maxover0p = 0; /* pretab */
1016     maxover1p = 0; /* pretab */
1017 
1018     for ( sfb = 0; sfb < SBPSY_l; sfb++ ) {
1019         v0  = (vbrmax - vbrsf->l[sfb]) - 2*max_range[sfb];
1020         v1  = (vbrmax - vbrsf->l[sfb]) - 4*max_range[sfb];
1021         v0p = (vbrmax - vbrsf->l[sfb]) - 2*(max_rangep[sfb]+pretab[sfb]);
1022         v1p = (vbrmax - vbrsf->l[sfb]) - 4*(max_rangep[sfb]+pretab[sfb]);
1023         if (maxover0 < v0)
1024             maxover0 = v0;
1025         if (maxover1 < v1)
1026             maxover1 = v1;
1027         if (maxover0p < v0p)
1028             maxover0p = v0p;
1029         if (maxover1p < v1p)
1030             maxover1p = v1p;
1031     }
1032 
1033     mover = Min (maxover0, maxover0p);
1034     if (gfc->noise_shaping == 2) {
1035         /* allow scalefac_scale=1 */
1036         mover = Min (mover, maxover1);
1037         mover = Min (mover, maxover1p);
1038     }
1039 
1040     vbrmax    -= mover;
1041     maxover0  -= mover;
1042     maxover0p -= mover;
1043     maxover1  -= mover;
1044     maxover1p -= mover;
1045 
1046     if (maxover0 <= 0) {
1047         cod_info->scalefac_scale = 0;
1048         cod_info->preflag = 0;
1049         vbrmax -= maxover0;
1050     } else if (maxover0p <= 0) {
1051         cod_info->scalefac_scale = 0;
1052         cod_info->preflag = 1;
1053         vbrmax -= maxover0p;
1054     } else if (maxover1 == 0) {
1055         cod_info->scalefac_scale = 1;
1056         cod_info->preflag = 0;
1057     } else if (maxover1p == 0) {
1058         cod_info->scalefac_scale = 1;
1059         cod_info->preflag = 1;
1060     } else {
1061         assert(0); /* this should not happen */
1062     }
1063 
1064     /* sf =  (cod_info->global_gain-210.0) */
1065     cod_info->global_gain = vbrmax + 210;
1066     assert (cod_info->global_gain < 256);
1067 
1068     if (vbr_mtrh == gfp->VBR && cod_info->global_gain > 1) {
1069         /*  just to be safe, reduce global gain by one
1070          */
1071         cod_info->global_gain -= 1;
1072     }
1073 
1074     if (cod_info->global_gain > 255)
1075         cod_info->global_gain = 255;
1076 
1077     for (sfb = 0; sfb < SBPSY_l; sfb++)
1078         vbrsf->l[sfb] -= vbrmax;
1079 
1080     if ( gfc->is_mpeg1 == 1 )
1081         maxover = compute_scalefacs_long (vbrsf->l, cod_info, scalefac->l);
1082     else
1083         maxover = compute_scalefacs_long_lsf (vbrsf->l, cod_info, scalefac->l);
1084 
1085     assert (maxover <= 0);
1086 
1087     *VBRmax = vbrmax;
1088 }
1089 
1090 
1091 
1092 /***********************************************************************
1093  *
1094  *      calc_fac()
1095  *
1096  *  Mark Taylor 2000-??-??
1097  *  Robert Hegemann 2000-10-20 made functions of it
1098  *
1099  ***********************************************************************/
1100 
calc_fac(const int ifac)1101 static FLOAT8 calc_fac ( const int ifac )
1102 {
1103     if (ifac+210 < Q_MAX)
1104         return 1/IPOW20 (ifac+210);
1105     else
1106         return pow (2.0, 0.75*ifac/4.0);
1107 }
1108 
1109 
1110 
1111 /***********************************************************************
1112  *
1113  *  quantize xr34 based on scalefactors
1114  *
1115  *  calc_short_block_xr34
1116  *  calc_long_block_xr34
1117  *
1118  *  Mark Taylor 2000-??-??
1119  *  Robert Hegemann 2000-10-20 made functions of them
1120  *
1121  ***********************************************************************/
1122 
1123 static void
short_block_xr34(const lame_internal_flags * const gfc,const gr_info * const cod_info,const III_scalefac_t * const scalefac,const FLOAT8 xr34_orig[576],FLOAT8 xr34[576])1124 short_block_xr34 (
1125     const lame_internal_flags        * const gfc,
1126     const gr_info        * const cod_info,
1127     const III_scalefac_t * const scalefac,
1128     const FLOAT8                 xr34_orig[576],
1129           FLOAT8                 xr34     [576] )
1130 {
1131      int sfb, l, j, b;
1132     int    ifac, ifqstep, start, end;
1133     FLOAT8 fac;
1134 
1135     /* even though there is no scalefactor for sfb12
1136      * subblock gain affects upper frequencies too, that's why
1137      * we have to go up to SBMAX_s
1138      */
1139     ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
1140     for ( j = 0, sfb = 0; sfb < SBMAX_s; sfb++ ) {
1141         start = gfc->scalefac_band.s[ sfb ];
1142         end   = gfc->scalefac_band.s[ sfb+1 ];
1143         for (b = 0; b < 3; b++) {
1144             ifac = 8*cod_info->subblock_gain[b]+ifqstep*scalefac->s[sfb][b];
1145             fac = calc_fac( ifac );
1146             /*
1147              *  loop unrolled into "Duff's Device".  Robert Hegemann
1148              */
1149             l = (end-start+7) / 8;
1150             switch ((end-start) % 8) {
1151                 default:
1152                 case 0: do{ xr34[j] = xr34_orig[j]*fac; j++;
1153                 case 7:     xr34[j] = xr34_orig[j]*fac; j++;
1154                 case 6:     xr34[j] = xr34_orig[j]*fac; j++;
1155                 case 5:     xr34[j] = xr34_orig[j]*fac; j++;
1156                 case 4:     xr34[j] = xr34_orig[j]*fac; j++;
1157                 case 3:     xr34[j] = xr34_orig[j]*fac; j++;
1158                 case 2:     xr34[j] = xr34_orig[j]*fac; j++;
1159                 case 1:     xr34[j] = xr34_orig[j]*fac; j++; } while (--l);
1160             }
1161         }
1162     }
1163 }
1164 
1165 
1166 
1167 static void
long_block_xr34(const lame_internal_flags * const gfc,const gr_info * const cod_info,const III_scalefac_t * const scalefac,const FLOAT8 xr34_orig[576],FLOAT8 xr34[576])1168 long_block_xr34 (
1169     const lame_internal_flags        * const gfc,
1170     const gr_info        * const cod_info,
1171     const III_scalefac_t * const scalefac,
1172     const FLOAT8                 xr34_orig[576],
1173           FLOAT8                 xr34     [576] )
1174 {
1175      int sfb, l, j;
1176     int    ifac, ifqstep, start, end;
1177     FLOAT8 fac;
1178 
1179     ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;
1180     for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
1181 
1182         ifac = ifqstep*scalefac->l[sfb];
1183         if (cod_info->preflag)
1184             ifac += ifqstep*pretab[sfb];
1185 
1186         fac = calc_fac( ifac );
1187 
1188         start = gfc->scalefac_band.l[ sfb ];
1189         end   = gfc->scalefac_band.l[ sfb+1 ];
1190         /*
1191          *  loop unrolled into "Duff's Device".  Robert Hegemann
1192          */
1193         j = start;
1194         l = (end-start+7) / 8;
1195         switch ((end-start) % 8) {
1196             default:
1197             case 0: do{ xr34[j] = xr34_orig[j]*fac; j++;
1198             case 7:     xr34[j] = xr34_orig[j]*fac; j++;
1199             case 6:     xr34[j] = xr34_orig[j]*fac; j++;
1200             case 5:     xr34[j] = xr34_orig[j]*fac; j++;
1201             case 4:     xr34[j] = xr34_orig[j]*fac; j++;
1202             case 3:     xr34[j] = xr34_orig[j]*fac; j++;
1203             case 2:     xr34[j] = xr34_orig[j]*fac; j++;
1204             case 1:     xr34[j] = xr34_orig[j]*fac; j++; } while (--l);
1205         }
1206     }
1207 }
1208 
1209 
1210 
1211 
1212 
1213 
1214 
1215 
1216 
1217 /************************************************************************
1218  *
1219  * VBR_noise_shaping()
1220  *
1221  * compute scalefactors, l3_enc, and return number of bits needed to encode
1222  *
1223  * return code:    0   scalefactors were found with all noise < masking
1224  *
1225  *               n>0   scalefactors required too many bits.  global gain
1226  *                     was decreased by n
1227  *                     If n is large, we should probably recompute scalefacs
1228  *                     with a lower quality.
1229  *
1230  *               n<0   scalefactors used less than minbits.
1231  *                     global gain was increased by n.
1232  *                     If n is large, might want to recompute scalefacs
1233  *                     with a higher quality setting?
1234  *
1235  ************************************************************************/
1236 static int
VBR_noise_shaping(lame_global_flags * gfp,FLOAT8 xr[576],FLOAT8 xr34orig[576],III_psy_ratio * ratio,int l3_enc[576],int digital_silence,int minbits,int maxbits,III_scalefac_t * scalefac,III_psy_xmin * l3_xmin,int gr,int ch)1237 VBR_noise_shaping (
1238     lame_global_flags *gfp,
1239     FLOAT8             xr       [576],
1240     FLOAT8             xr34orig [576],
1241     III_psy_ratio     *ratio,
1242     int                l3_enc   [576],
1243     int                digital_silence,
1244     int                minbits,
1245     int                maxbits,
1246     III_scalefac_t    *scalefac,
1247     III_psy_xmin      *l3_xmin,
1248     int                gr,
1249     int                ch )
1250 {
1251     lame_internal_flags *gfc=gfp->internal_flags;
1252     III_scalefac_t save_sf;
1253     III_scalefac_t vbrsf;
1254     gr_info *cod_info;
1255     FLOAT8 xr34[576];
1256     int shortblock;
1257     int vbrmax;
1258     int global_gain_adjust = 0;
1259 
1260     cod_info   = &gfc->l3_side.gr[gr].ch[ch].tt;
1261     shortblock = (cod_info->block_type == SHORT_TYPE);
1262 
1263     if (shortblock)
1264         vbrmax = short_block_vbr_sf (gfc, l3_xmin, xr34orig, xr, &vbrsf);
1265     else
1266         vbrmax = long_block_vbr_sf (gfc, l3_xmin, xr34orig, xr, &vbrsf);
1267 
1268     /* save a copy of vbrsf, incase we have to recomptue scalefacs */
1269     memcpy (&save_sf, &vbrsf, sizeof(III_scalefac_t));
1270 
1271     do {
1272         memset (scalefac, 0, sizeof(III_scalefac_t));
1273 
1274         if (shortblock) {
1275             short_block_scalefacs (gfp, cod_info, scalefac, &vbrsf, &vbrmax);
1276             short_block_xr34      (gfc, cod_info, scalefac, xr34orig, xr34);
1277         } else {
1278             long_block_scalefacs (gfp, cod_info, scalefac, &vbrsf, &vbrmax);
1279             long_block_xr34      (gfc, cod_info, scalefac, xr34orig, xr34);
1280         }
1281         VBR_quantize_granule (gfp, xr34, l3_enc, ratio, scalefac, gr, ch);
1282 
1283 
1284         /* decrease noise until we use at least minbits
1285          */
1286         if (cod_info->part2_3_length < minbits) {
1287             if (digital_silence) break;
1288             //if (cod_info->part2_3_length == cod_info->part2_length) break;
1289             if (vbrmax+210 == 0) break;
1290 
1291             /* decrease global gain, recompute scale factors */
1292             --vbrmax;
1293             --global_gain_adjust;
1294             memcpy (&vbrsf, &save_sf, sizeof(III_scalefac_t));
1295         }
1296 
1297     } while (cod_info->part2_3_length < minbits);
1298 
1299     /* inject noise until we meet our bit limit
1300      */
1301     while (cod_info->part2_3_length > Min (maxbits, MAX_BITS)) {
1302         /* increase global gain, keep existing scale factors */
1303         ++cod_info->global_gain;
1304         if (cod_info->global_gain > 255)
1305             ERRORF (gfc,"%ld impossible to encode ??? frame! bits=%d\n",
1306                     //  gfp->frameNum, cod_info->part2_3_length);
1307                              -1,       cod_info->part2_3_length);
1308         VBR_quantize_granule (gfp, xr34, l3_enc, ratio, scalefac, gr, ch);
1309 
1310         ++global_gain_adjust;
1311     }
1312 
1313     return global_gain_adjust;
1314 }
1315 
1316 
1317 
1318 /************************************************************************
1319  *
1320  *  VBR_noise_shaping2()
1321  *
1322  *  may result in a need of too many bits, then do it CBR like
1323  *
1324  *  Robert Hegemann 2000-10-25
1325  *
1326  ***********************************************************************/
1327 
1328 int
VBR_noise_shaping2(lame_global_flags * gfp,FLOAT8 xr[576],FLOAT8 xr34orig[576],III_psy_ratio * const ratio,int l3_enc[576],int digital_silence,int minbits,int maxbits,III_scalefac_t * const scalefac,III_psy_xmin * const l3_xmin,int gr,int ch)1329 VBR_noise_shaping2 (
1330     lame_global_flags        *gfp,
1331     FLOAT8                 xr       [576],
1332     FLOAT8                 xr34orig [576],
1333     III_psy_ratio  * const ratio,
1334     int                    l3_enc   [576],
1335     int                    digital_silence,
1336     int                    minbits,
1337     int                    maxbits,
1338     III_scalefac_t * const scalefac,
1339     III_psy_xmin   * const l3_xmin,
1340     int                    gr,
1341     int                    ch )
1342 {
1343     lame_internal_flags *gfc=gfp->internal_flags;
1344     III_scalefac_t vbrsf;
1345     gr_info *cod_info;
1346     FLOAT8 xr34[576];
1347     int shortblock, ret, bits, huffbits;
1348     int vbrmax, best_huffman = gfc->use_best_huffman;
1349 
1350     cod_info   = &gfc->l3_side.gr[gr].ch[ch].tt;
1351     shortblock = (cod_info->block_type == SHORT_TYPE);
1352 
1353     if (shortblock) {
1354         vbrmax = short_block_sf (gfc, l3_xmin, xr34orig, xr, &vbrsf);
1355         short_block_scalefacs (gfp, cod_info, scalefac, &vbrsf, &vbrmax);
1356         short_block_xr34      (gfc, cod_info, scalefac, xr34orig, xr34);
1357     } else {
1358         vbrmax = long_block_sf (gfc, l3_xmin, xr34orig, xr, &vbrsf);
1359         long_block_scalefacs (gfp, cod_info, scalefac, &vbrsf, &vbrmax);
1360         long_block_xr34      (gfc, cod_info, scalefac, xr34orig, xr34);
1361     }
1362 
1363     gfc->use_best_huffman = 0; /* we will do it later */
1364 
1365     ret = VBR_quantize_granule (gfp, xr34, l3_enc, ratio, scalefac, gr, ch);
1366 
1367     gfc->use_best_huffman = best_huffman;
1368 
1369     if (ret == -1) /* Houston, we have a problem */
1370         return -1;
1371 
1372     if (cod_info->part2_3_length < minbits) {
1373         huffbits = minbits - cod_info->part2_length;
1374         bits = bin_search_StepSize (gfc, cod_info, huffbits,
1375                                     gfc->OldValue[ch], xr34, l3_enc);
1376         gfc->OldValue[ch] = cod_info->global_gain;
1377         cod_info->part2_3_length  = bits + cod_info->part2_length;
1378     }
1379     if (cod_info->part2_3_length > maxbits) {
1380         huffbits = maxbits - cod_info->part2_length;
1381         bits = bin_search_StepSize (gfc, cod_info, huffbits,
1382                                     gfc->OldValue[ch], xr34, l3_enc);
1383         gfc->OldValue[ch] = cod_info->global_gain;
1384         cod_info->part2_3_length = bits;
1385         if (bits > huffbits) {
1386             bits = inner_loop (gfc, cod_info, huffbits, xr34, l3_enc);
1387             cod_info->part2_3_length  = bits;
1388         }
1389         if (bits >= LARGE_BITS) /* Houston, we have a problem */
1390             return -2;
1391         cod_info->part2_3_length += cod_info->part2_length;
1392     }
1393 
1394     if (cod_info->part2_length >= LARGE_BITS) /* Houston, we have a problem */
1395         return -2;
1396 
1397     assert (cod_info->global_gain < 256);
1398 
1399     return 0;
1400 }
1401 
1402 
1403 
1404 
1405 void
VBR_quantize(lame_global_flags * gfp,FLOAT8 pe[2][2],FLOAT8 ms_ener_ratio[2],FLOAT8 xr[2][2][576],III_psy_ratio ratio[2][2],int l3_enc[2][2][576],III_scalefac_t scalefac[2][2])1406 VBR_quantize(lame_global_flags *gfp,
1407                 FLOAT8 pe[2][2], FLOAT8 ms_ener_ratio[2],
1408                 FLOAT8 xr[2][2][576], III_psy_ratio ratio[2][2],
1409                 int l3_enc[2][2][576],
1410                 III_scalefac_t scalefac[2][2])
1411 {
1412   lame_internal_flags *gfc=gfp->internal_flags;
1413   III_psy_xmin l3_xmin[2][2];
1414   int minbits,maxbits,max_frame_bits,totbits,gr,ch,i,bits_ok;
1415   int bitsPerFrame,mean_bits;
1416   int analog_silence;
1417   FLOAT8 qadjust;
1418   III_side_info_t * l3_side;
1419   gr_info *cod_info;
1420   int digital_silence[2][2];
1421   FLOAT8 masking_lower_db=0;
1422   FLOAT8 xr34[2][2][576];
1423   // static const FLOAT8 dbQ[10]={-6.0,-5.0,-4.0,-3.0, -2.0, -1.0, -.25, .5, 1.25, 2.0};
1424   /* from quantize.c VBR algorithm */
1425   /*static const FLOAT8 dbQ[10]=
1426    {-5.5,-4.25,-3.0,-2.50, -1.75, -.75, -.5, -.25, .25, .75};*/
1427   /* a third dbQ table ?!? */
1428   static const FLOAT8 dbQ[10]=
1429   {-6.06,-4.4,-2.9,-1.57, -0.4, 0.61, 1.45, 2.13, 2.65, 3.0};
1430 
1431   qadjust=0;   /* start with -1 db quality improvement over quantize.c VBR */
1432 
1433   l3_side = &gfc->l3_side;
1434   //gfc->ATHlower += (4-gfp->VBR_q)*4.0;
1435   //if (gfc->ATHlower < 0) gfc->ATHlower=0;
1436 
1437 
1438   /* now find out: if the frame can be considered analog silent
1439    *               if each granule can be considered digital silent
1440    * and calculate l3_xmin and the fresh xr34 array
1441    */
1442 
1443   assert( gfp->VBR_q <= 9 );
1444   assert( gfp->VBR_q >= 0 );
1445   analog_silence=1;
1446   for (gr = 0; gr < gfc->mode_gr; gr++) {
1447     /* copy data to be quantized into xr */
1448     if (gfc->mode_ext==MPG_MD_MS_LR) {
1449       ms_convert(xr[gr],xr[gr]);
1450     }
1451     for (ch = 0; ch < gfc->channels_out; ch++) {
1452       /* if in the following sections the quality would not be adjusted
1453        * then we would only have to call calc_xmin once here and
1454        * could drop subsequently calls (rh 2000/07/17)
1455        */
1456       int over_ath;
1457       cod_info = &l3_side->gr[gr].ch[ch].tt;
1458       cod_info->part2_3_length=LARGE_BITS;
1459 
1460       if (cod_info->block_type == SHORT_TYPE) {
1461           cod_info->sfb_lmax = 0; /* No sb*/
1462           cod_info->sfb_smin = 0;
1463       } else {
1464           /* MPEG 1 doesnt use last scalefactor band */
1465           cod_info->sfb_lmax = SBPSY_l;
1466           cod_info->sfb_smin = SBPSY_s;    /* No sb */
1467 	  if (cod_info->mixed_block_flag) {
1468 	    cod_info->sfb_lmax        = 8;
1469 	    cod_info->sfb_smin        = 3;
1470 	  }
1471       }
1472 
1473       /* quality setting */
1474       masking_lower_db = dbQ[gfp->VBR_q];
1475       if (pe[gr][ch]>750) {
1476         masking_lower_db -= Min(10,4*(pe[gr][ch]-750.)/750.);
1477       }
1478       gfc->masking_lower = pow(10.0,masking_lower_db/10);
1479 
1480       /* masking thresholds */
1481       over_ath = calc_xmin(gfp,xr[gr][ch],&ratio[gr][ch],cod_info,&l3_xmin[gr][ch]);
1482 
1483       /* if there are bands with more energy than the ATH
1484        * then we say the frame is not analog silent */
1485       if (over_ath) {
1486         analog_silence = 0;
1487       }
1488 
1489       /* if there is no line with more energy than 1e-20
1490        * then this granule is considered to be digital silent
1491        * plus calculation of xr34 */
1492       digital_silence[gr][ch] = 1;
1493       for(i=0;i<576;i++) {
1494         FLOAT8 temp=fabs(xr[gr][ch][i]);
1495         xr34[gr][ch][i]=sqrt(sqrt(temp)*temp);
1496         digital_silence[gr][ch] &= temp < 1E-20;
1497       }
1498     } /* ch */
1499   }  /* gr */
1500 
1501 
1502   /* compute minimum allowed bits from minimum allowed bitrate */
1503   if (analog_silence) {
1504     gfc->bitrate_index=1;
1505   } else {
1506     gfc->bitrate_index=gfc->VBR_min_bitrate;
1507   }
1508   getframebits(gfp, &bitsPerFrame, &mean_bits);
1509   minbits = (mean_bits/gfc->channels_out);
1510 
1511   /* compute maximum allowed bits from max allowed bitrate */
1512   gfc->bitrate_index=gfc->VBR_max_bitrate;
1513   getframebits(gfp, &bitsPerFrame, &mean_bits);
1514   max_frame_bits = ResvFrameBegin(gfp, l3_side, mean_bits, bitsPerFrame);
1515   maxbits=2.5*(mean_bits/gfc->channels_out);
1516 
1517   {
1518   /* compute a target  mean_bits based on compression ratio
1519    * which was set based on VBR_q
1520    */
1521   int bit_rate = gfp->out_samplerate*16*gfc->channels_out/(1000.0*gfp->compression_ratio);
1522   bitsPerFrame = (bit_rate*gfp->framesize*1000)/gfp->out_samplerate;
1523   mean_bits = (bitsPerFrame - 8*gfc->sideinfo_len) / gfc->mode_gr;
1524   }
1525 
1526 
1527   minbits = Max(minbits,125);
1528   minbits=Max(minbits,.40*(mean_bits/gfc->channels_out));
1529   maxbits=Min(maxbits,2.5*(mean_bits/gfc->channels_out));
1530 
1531 
1532 
1533 
1534 
1535 
1536 
1537   /*
1538    * loop over all ch,gr, encoding anything with bits > .5*(max_frame_bits/4)
1539    *
1540    * If a particular granule uses way too many bits, it will be re-encoded
1541    * on the next iteration of the loop (with a lower quality setting).
1542    * But granules which dont use
1543    * use too many bits will not be re-encoded.
1544    *
1545    * minbits:  minimum allowed bits for 1 granule 1 channel
1546    * maxbits:  maximum allowwed bits for 1 granule 1 channel
1547    * max_frame_bits:  maximum allowed bits for entire frame
1548    * (max_frame_bits/4)   estimate of average bits per granule per channel
1549    *
1550    */
1551 
1552   do {
1553 
1554     totbits=0;
1555     for (gr = 0; gr < gfc->mode_gr; gr++) {
1556       int minbits_lr[2];
1557       minbits_lr[0]=minbits;
1558       minbits_lr[1]=minbits;
1559 
1560 #if 0
1561       if (gfc->mode_ext==MPG_MD_MS_LR) {
1562         FLOAT8 fac;
1563         fac = .33*(.5-ms_ener_ratio[gr])/.5;
1564         if (fac<0) fac=0;
1565         if (fac>.5) fac=.5;
1566         minbits_lr[0] = (1+fac)*minbits;
1567         minbits_lr[1] = Max(125,(1-fac)*minbits);
1568       }
1569 #endif
1570 
1571 
1572       for (ch = 0; ch < gfc->channels_out; ch++) {
1573         int adjusted,shortblock;
1574         cod_info = &l3_side->gr[gr].ch[ch].tt;
1575 
1576         /* ENCODE this data first pass, and on future passes unless it uses
1577          * a very small percentage of the max_frame_bits  */
1578         if (cod_info->part2_3_length > (max_frame_bits/(2*gfc->channels_out*gfc->mode_gr))) {
1579 
1580           shortblock = (cod_info->block_type == SHORT_TYPE);
1581 
1582           /* Adjust allowed masking based on quality setting */
1583           if (qadjust!=0 /*|| shortblock*/) {
1584             masking_lower_db = dbQ[gfp->VBR_q] + qadjust;
1585 
1586             /*
1587             if (shortblock) masking_lower_db -= 4;
1588             */
1589 
1590             if (pe[gr][ch]>750)
1591               masking_lower_db -= Min(10,4*(pe[gr][ch]-750.)/750.);
1592             gfc->masking_lower = pow(10.0,masking_lower_db/10);
1593             calc_xmin( gfp, xr[gr][ch], ratio[gr]+ch, cod_info, l3_xmin[gr]+ch);
1594           }
1595 
1596           /* digital silent granules do not need the full round trip,
1597            * but this can be optimized later on
1598            */
1599           adjusted = VBR_noise_shaping (gfp,xr[gr][ch],xr34[gr][ch],
1600                                         ratio[gr]+ch,l3_enc[gr][ch],
1601                                         digital_silence[gr][ch],
1602                                         minbits_lr[ch],
1603                                         maxbits,scalefac[gr]+ch,
1604                                         l3_xmin[gr]+ch,gr,ch);
1605           if (adjusted>10) {
1606             /* global_gain was changed by a large amount to get bits < maxbits */
1607             /* quality is set to high.  we could set bits = LARGE_BITS
1608              * to force re-encoding.  But most likely the other channels/granules
1609              * will also use too many bits, and the entire frame will
1610              * be > max_frame_bits, forcing re-encoding below.
1611              */
1612             // cod_info->part2_3_bits = LARGE_BITS;
1613           }
1614         }
1615         totbits += cod_info->part2_3_length;
1616       }
1617     }
1618     bits_ok=1;
1619     if (totbits>max_frame_bits) {
1620       /* lower quality */
1621       qadjust += Max(.125,Min(1,(totbits-max_frame_bits)/300.0));
1622       /* adjusting minbits and maxbits is necessary too
1623        * cos lowering quality is not enough in rare cases
1624        * when each granule still needs almost maxbits, it wont fit */
1625       minbits = Max(125,minbits*0.975);
1626       maxbits = Max(minbits,maxbits*0.975);
1627       //      DEBUGF("%i totbits>max_frame_bits   totbits=%i  maxbits=%i \n",gfp->frameNum,totbits,max_frame_bits);
1628       //      DEBUGF("next masking_lower_db = %f \n",masking_lower_db + qadjust);
1629       bits_ok=0;
1630     }
1631 
1632   } while (!bits_ok);
1633 
1634 
1635 
1636   /* find optimal scalefac storage.  Cant be done above because
1637    * might enable scfsi which breaks the interation loops */
1638   totbits=0;
1639   for (gr = 0; gr < gfc->mode_gr; gr++) {
1640     for (ch = 0; ch < gfc->channels_out; ch++) {
1641       best_scalefac_store(gfc, gr, ch, l3_enc, l3_side, scalefac);
1642       totbits += l3_side->gr[gr].ch[ch].tt.part2_3_length;
1643     }
1644   }
1645 
1646 
1647 
1648 
1649   if (analog_silence && !gfp->VBR_hard_min) {
1650     gfc->bitrate_index = 1;
1651   } else {
1652     gfc->bitrate_index = gfc->VBR_min_bitrate;
1653   }
1654   for( ; gfc->bitrate_index < gfc->VBR_max_bitrate; gfc->bitrate_index++ ) {
1655 
1656     getframebits (gfp, &bitsPerFrame, &mean_bits);
1657     maxbits = ResvFrameBegin(gfp, l3_side, mean_bits, bitsPerFrame);
1658     if (totbits <= maxbits) break;
1659   }
1660   if (gfc->bitrate_index == gfc->VBR_max_bitrate) {
1661     getframebits (gfp, &bitsPerFrame, &mean_bits);
1662     maxbits = ResvFrameBegin(gfp, l3_side, mean_bits, bitsPerFrame);
1663   }
1664 
1665   //  DEBUGF("%i total_bits=%i max_frame_bits=%i index=%i  \n",gfp->frameNum,totbits,max_frame_bits,gfc->bitrate_index);
1666 
1667   for (gr = 0; gr < gfc->mode_gr; gr++) {
1668     for (ch = 0; ch < gfc->channels_out; ch++) {
1669       cod_info = &l3_side->gr[gr].ch[ch].tt;
1670 
1671 
1672       ResvAdjust (gfc, cod_info, l3_side, mean_bits);
1673 
1674       /*******************************************************************
1675        * set the sign of l3_enc from the sign of xr
1676        *******************************************************************/
1677       for ( i = 0; i < 576; i++) {
1678         if (xr[gr][ch][i] < 0) l3_enc[gr][ch][i] *= -1;
1679       }
1680     }
1681   }
1682   ResvFrameEnd (gfc, l3_side, mean_bits);
1683 
1684 
1685 
1686 }
1687 
1688 
1689 
1690