python 2.6 support
[swftools.git] / lib / lame / quantize_pvt.c
1 /*
2  *      quantize_pvt source file
3  *
4  *      Copyright (c) 1999 Takehiro TOMINAGA
5  *
6  * This library is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU Library General Public
8  * License as published by the Free Software Foundation; either
9  * version 2 of the License, or (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14  * Library General Public License for more details.
15  *
16  * You should have received a copy of the GNU Library General Public
17  * License along with this library; if not, write to the
18  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
19  * Boston, MA 02111-1307, USA.
20  */
21
22 /* $Id: quantize_pvt.c,v 1.2 2006/02/09 16:56:23 kramm Exp $ */
23 #include <stdlib.h>
24 #include "config_static.h"
25
26 #include <assert.h>
27 #include "util.h"
28 #include "lame-analysis.h"
29 #include "tables.h"
30 #include "reservoir.h"
31 #include "quantize_pvt.h"
32
33 #ifdef WITH_DMALLOC
34 #include <dmalloc.h>
35 #endif
36
37
38 #define NSATHSCALE 100 // Assuming dynamic range=96dB, this value should be 92
39
40 const char  slen1_tab [16] = { 0, 0, 0, 0, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4 };
41 const char  slen2_tab [16] = { 0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3, 1, 2, 3, 2, 3 };
42
43
44 /*
45   The following table is used to implement the scalefactor
46   partitioning for MPEG2 as described in section
47   2.4.3.2 of the IS. The indexing corresponds to the
48   way the tables are presented in the IS:
49
50   [table_number][row_in_table][column of nr_of_sfb]
51 */
52 const int  nr_of_sfb_block [6] [3] [4] =
53 {
54   {
55     {6, 5, 5, 5},
56     {9, 9, 9, 9},
57     {6, 9, 9, 9}
58   },
59   {
60     {6, 5, 7, 3},
61     {9, 9, 12, 6},
62     {6, 9, 12, 6}
63   },
64   {
65     {11, 10, 0, 0},
66     {18, 18, 0, 0},
67     {15,18,0,0}
68   },
69   {
70     {7, 7, 7, 0},
71     {12, 12, 12, 0},
72     {6, 15, 12, 0}
73   },
74   {
75     {6, 6, 6, 3},
76     {12, 9, 9, 6},
77     {6, 12, 9, 6}
78   },
79   {
80     {8, 8, 5, 0},
81     {15,12,9,0},
82     {6,18,9,0}
83   }
84 };
85
86
87 /* Table B.6: layer3 preemphasis */
88 const char  pretab [SBMAX_l] =
89 {
90     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
91     1, 1, 1, 1, 2, 2, 3, 3, 3, 2, 0
92 };
93
94 /*
95   Here are MPEG1 Table B.8 and MPEG2 Table B.1
96   -- Layer III scalefactor bands. 
97   Index into this using a method such as:
98     idx  = fr_ps->header->sampling_frequency
99            + (fr_ps->header->version * 3)
100 */
101
102
103 const scalefac_struct sfBandIndex[9] =
104 {
105   { /* Table B.2.b: 22.05 kHz */
106     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
107     {0,4,8,12,18,24,32,42,56,74,100,132,174,192}
108   },
109   { /* Table B.2.c: 24 kHz */                 /* docs: 332. mpg123(broken): 330 */
110     {0,6,12,18,24,30,36,44,54,66,80,96,114,136,162,194,232,278, 332, 394,464,540,576},
111     {0,4,8,12,18,26,36,48,62,80,104,136,180,192}
112   },
113   { /* Table B.2.a: 16 kHz */
114     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
115     {0,4,8,12,18,26,36,48,62,80,104,134,174,192}
116   },
117   { /* Table B.8.b: 44.1 kHz */
118     {0,4,8,12,16,20,24,30,36,44,52,62,74,90,110,134,162,196,238,288,342,418,576},
119     {0,4,8,12,16,22,30,40,52,66,84,106,136,192}
120   },
121   { /* Table B.8.c: 48 kHz */
122     {0,4,8,12,16,20,24,30,36,42,50,60,72,88,106,128,156,190,230,276,330,384,576},
123     {0,4,8,12,16,22,28,38,50,64,80,100,126,192}
124   },
125   { /* Table B.8.a: 32 kHz */
126     {0,4,8,12,16,20,24,30,36,44,54,66,82,102,126,156,194,240,296,364,448,550,576},
127     {0,4,8,12,16,22,30,42,58,78,104,138,180,192}
128   },
129   { /* MPEG-2.5 11.025 kHz */
130     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
131     {0/3,12/3,24/3,36/3,54/3,78/3,108/3,144/3,186/3,240/3,312/3,402/3,522/3,576/3}
132   },
133   { /* MPEG-2.5 12 kHz */
134     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
135     {0/3,12/3,24/3,36/3,54/3,78/3,108/3,144/3,186/3,240/3,312/3,402/3,522/3,576/3}
136   },
137   { /* MPEG-2.5 8 kHz */
138     {0,12,24,36,48,60,72,88,108,132,160,192,232,280,336,400,476,566,568,570,572,574,576},
139     {0/3,24/3,48/3,72/3,108/3,156/3,216/3,288/3,372/3,480/3,486/3,492/3,498/3,576/3}
140   }
141 };
142
143
144
145 FLOAT8 pow20[Q_MAX];
146 FLOAT8 ipow20[Q_MAX];
147 FLOAT8 iipow20[Q_MAX];
148 FLOAT8 *iipow20_;
149 FLOAT8 pow43[PRECALC_SIZE];
150 /* initialized in first call to iteration_init */
151 FLOAT8 adj43asm[PRECALC_SIZE];
152 FLOAT8 adj43[PRECALC_SIZE];
153
154 /************************************************************************/
155 /*  initialization for iteration_loop */
156 /************************************************************************/
157 void
158 iteration_init( lame_global_flags *gfp)
159 {
160   lame_internal_flags *gfc=gfp->internal_flags;
161   III_side_info_t * const l3_side = &gfc->l3_side;
162   int i;
163
164   if ( gfc->iteration_init_init==0 ) {
165     gfc->iteration_init_init=1;
166
167     l3_side->main_data_begin = 0;
168     compute_ath(gfp,gfc->ATH->l,gfc->ATH->s);
169
170     pow43[0] = 0.0;
171     for(i=1;i<PRECALC_SIZE;i++)
172         pow43[i] = pow((FLOAT8)i, 4.0/3.0);
173
174     adj43asm[0] = 0.0;
175     for (i = 1; i < PRECALC_SIZE; i++)
176       adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]),0.75);
177     for (i = 0; i < PRECALC_SIZE-1; i++)
178         adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
179     adj43[i] = 0.5;
180     iipow20_ = &iipow20[210];
181     for (i = 0; i < Q_MAX; i++) {
182         iipow20[i] = pow(2.0, (double)(i - 210) * 0.1875);
183         ipow20[i] = pow(2.0, (double)(i - 210) * -0.1875);
184         pow20[i] = pow(2.0, (double)(i - 210) * 0.25);
185     }
186     huffman_init(gfc);
187   }
188 }
189
190
191
192
193
194 /* 
195 compute the ATH for each scalefactor band 
196 cd range:  0..96db
197
198 Input:  3.3kHz signal  32767 amplitude  (3.3kHz is where ATH is smallest = -5db)
199 longblocks:  sfb=12   en0/bw=-11db    max_en0 = 1.3db
200 shortblocks: sfb=5           -9db              0db
201
202 Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
203 longblocks:  amp=1      sfb=12   en0/bw=-103 db      max_en0 = -92db
204             amp=32767   sfb=12           -12 db                 -1.4db 
205
206 Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
207 shortblocks: amp=1      sfb=5   en0/bw= -99                    -86 
208             amp=32767   sfb=5           -9  db                  4db 
209
210
211 MAX energy of largest wave at 3.3kHz = 1db
212 AVE energy of largest wave at 3.3kHz = -11db
213 Let's take AVE:  -11db = maximum signal in sfb=12.  
214 Dynamic range of CD: 96db.  Therefor energy of smallest audible wave 
215 in sfb=12  = -11  - 96 = -107db = ATH at 3.3kHz.  
216
217 ATH formula for this wave: -5db.  To adjust to LAME scaling, we need
218 ATH = ATH_formula  - 103  (db)
219 ATH = ATH * 2.5e-10      (ener)
220
221 */
222
223 FLOAT8 ATHmdct( lame_global_flags *gfp, FLOAT8 f )
224 {
225     lame_internal_flags *gfc = gfp->internal_flags;
226     FLOAT8 ath;
227   
228     ath = ATHformula( f , gfp );
229           
230     if (gfc->nsPsy.use) {
231         ath -= NSATHSCALE;
232     } else {
233         ath -= 114;
234     }
235     
236     /*  modify the MDCT scaling for the ATH  */
237     ath -= gfp->ATHlower;
238
239     /* convert to energy */
240     ath = pow( 10.0, ath/10.0 );
241     return ath;
242 }
243  
244
245 void compute_ath( lame_global_flags *gfp, FLOAT8 ATH_l[], FLOAT8 ATH_s[] )
246 {
247     lame_internal_flags *gfc = gfp->internal_flags;
248     int sfb, i, start, end;
249     FLOAT8 ATH_f;
250     FLOAT8 samp_freq = gfp->out_samplerate;
251
252     for (sfb = 0; sfb < SBMAX_l; sfb++) {
253         start = gfc->scalefac_band.l[ sfb ];
254         end   = gfc->scalefac_band.l[ sfb+1 ];
255         ATH_l[sfb]=FLOAT8_MAX;
256         for (i = start ; i < end; i++) {
257             FLOAT8 freq = i*samp_freq/(2*576);
258             ATH_f = ATHmdct( gfp, freq );  /* freq in kHz */
259             ATH_l[sfb] = Min( ATH_l[sfb], ATH_f );
260         }
261     }
262
263     for (sfb = 0; sfb < SBMAX_s; sfb++){
264         start = gfc->scalefac_band.s[ sfb ];
265         end   = gfc->scalefac_band.s[ sfb+1 ];
266         ATH_s[sfb] = FLOAT8_MAX;
267         for (i = start ; i < end; i++) {
268             FLOAT8 freq = i*samp_freq/(2*192);
269             ATH_f = ATHmdct( gfp, freq );    /* freq in kHz */
270             ATH_s[sfb] = Min( ATH_s[sfb], ATH_f );
271         }
272     }
273
274     /*  no-ATH mode:
275      *  reduce ATH to -200 dB
276      */
277     
278     if (gfp->noATH) {
279         for (sfb = 0; sfb < SBMAX_l; sfb++) {
280             ATH_l[sfb] = 1E-37;
281         }
282         for (sfb = 0; sfb < SBMAX_s; sfb++) {
283             ATH_s[sfb] = 1E-37;
284         }
285     }
286     
287     /*  work in progress, don't rely on it too much
288      */
289     gfc->ATH-> floor = 10. * log10( ATHmdct( gfp, -1. ) );
290     
291     /*
292     {   FLOAT8 g=10000, t=1e30, x;
293         for ( f = 100; f < 10000; f++ ) {
294             x = ATHmdct( gfp, f );
295             if ( t > x ) t = x, g = f;
296         }
297         printf("min=%g\n", g);
298     }*/
299 }
300
301
302
303
304
305 /* convert from L/R <-> Mid/Side, src == dst allowed */
306 void ms_convert(FLOAT8 dst[2][576], FLOAT8 src[2][576])
307 {
308     FLOAT8 l;
309     FLOAT8 r;
310     int i;
311     for (i = 0; i < 576; ++i) {
312         l = src[0][i];
313         r = src[1][i];
314         dst[0][i] = (l+r) * (FLOAT8)(SQRT2*0.5);
315         dst[1][i] = (l-r) * (FLOAT8)(SQRT2*0.5);
316     }
317 }
318
319
320
321 /************************************************************************
322  * allocate bits among 2 channels based on PE
323  * mt 6/99
324  * bugfixes rh 8/01: often allocated more than the allowed 4095 bits
325  ************************************************************************/
326 int on_pe( lame_global_flags *gfp, FLOAT8 pe[][2], III_side_info_t *l3_side,
327            int targ_bits[2], int mean_bits, int gr )
328 {
329     lame_internal_flags * gfc = gfp->internal_flags;
330     gr_info *   cod_info;
331     int     extra_bits, tbits, bits;
332     int     add_bits[2]; 
333     int     max_bits;  /* maximum allowed bits for this granule */
334     int     ch;
335
336     /* allocate targ_bits for granule */
337     ResvMaxBits( gfp, mean_bits, &tbits, &extra_bits );
338     max_bits = tbits + extra_bits;
339     if (max_bits > MAX_BITS) /* hard limit per granule */
340         max_bits = MAX_BITS;
341     
342     for ( bits = 0, ch = 0; ch < gfc->channels_out; ++ch ) {
343         /******************************************************************
344          * allocate bits for each channel 
345          ******************************************************************/
346         cod_info = &l3_side->gr[gr].ch[ch].tt;
347     
348         targ_bits[ch] = Min( MAX_BITS, tbits/gfc->channels_out );
349     
350         if (gfc->nsPsy.use) {
351             add_bits[ch] = targ_bits[ch] * pe[gr][ch] / 700.0 - targ_bits[ch];
352         } 
353         else {
354             add_bits[ch] = (pe[gr][ch]-750) / 1.4;
355             /* short blocks us a little extra, no matter what the pe */
356             if ( cod_info->block_type == SHORT_TYPE ) {
357                 if (add_bits[ch] < mean_bits/4) 
358                     add_bits[ch] = mean_bits/4;
359             }
360         }
361
362         /* at most increase bits by 1.5*average */
363         if (add_bits[ch] > mean_bits*3/4) 
364             add_bits[ch] = mean_bits*3/4;
365         if (add_bits[ch] < 0) 
366             add_bits[ch] = 0;
367
368         if (add_bits[ch]+targ_bits[ch] > MAX_BITS) 
369             add_bits[ch] = Max( 0, MAX_BITS-targ_bits[ch] );
370
371         bits += add_bits[ch];
372     }
373     if (bits > extra_bits) {
374         for ( ch = 0; ch < gfc->channels_out; ++ch ) {
375             add_bits[ch] = extra_bits * add_bits[ch] / bits;
376         }
377     }
378
379     for ( ch = 0; ch < gfc->channels_out; ++ch ) {
380         targ_bits[ch] += add_bits[ch];
381         extra_bits    -= add_bits[ch];
382         assert( targ_bits[ch] <= MAX_BITS );
383     }
384     assert( max_bits <= MAX_BITS );
385     return max_bits;
386 }
387
388
389
390
391 void reduce_side(int targ_bits[2],FLOAT8 ms_ener_ratio,int mean_bits,int max_bits)
392 {
393   int move_bits;
394   FLOAT fac;
395
396
397   /*  ms_ener_ratio = 0:  allocate 66/33  mid/side  fac=.33  
398    *  ms_ener_ratio =.5:  allocate 50/50 mid/side   fac= 0 */
399   /* 75/25 split is fac=.5 */
400   /* float fac = .50*(.5-ms_ener_ratio[gr])/.5;*/
401   fac = .33*(.5-ms_ener_ratio)/.5;
402   if (fac<0) fac=0;
403   if (fac>.5) fac=.5;
404   
405     /* number of bits to move from side channel to mid channel */
406     /*    move_bits = fac*targ_bits[1];  */
407     move_bits = fac*.5*(targ_bits[0]+targ_bits[1]);  
408
409     if (move_bits > MAX_BITS - targ_bits[0]) {
410         move_bits = MAX_BITS - targ_bits[0];
411     }
412     if (move_bits<0) move_bits=0;
413     
414     if (targ_bits[1] >= 125) {
415       /* dont reduce side channel below 125 bits */
416       if (targ_bits[1]-move_bits > 125) {
417
418         /* if mid channel already has 2x more than average, dont bother */
419         /* mean_bits = bits per granule (for both channels) */
420         if (targ_bits[0] < mean_bits)
421           targ_bits[0] += move_bits;
422         targ_bits[1] -= move_bits;
423       } else {
424         targ_bits[0] += targ_bits[1] - 125;
425         targ_bits[1] = 125;
426       }
427     }
428     
429     move_bits=targ_bits[0]+targ_bits[1];
430     if (move_bits > max_bits) {
431       targ_bits[0]=(max_bits*targ_bits[0])/move_bits;
432       targ_bits[1]=(max_bits*targ_bits[1])/move_bits;
433     }
434     assert (targ_bits[0] <= MAX_BITS);
435     assert (targ_bits[1] <= MAX_BITS);
436 }
437
438
439 /**
440  *  Robert Hegemann 2001-04-27:
441  *  this adjusts the ATH, keeping the original noise floor
442  *  affects the higher frequencies more than the lower ones
443  */
444
445 FLOAT8 athAdjust( FLOAT8 a, FLOAT8 x, FLOAT8 athFloor )
446 {
447     /*  work in progress
448      */
449     FLOAT8 const o = 90.30873362;
450     FLOAT8 const p = 94.82444863;
451     FLOAT8 u = 10. * log10(x); 
452     FLOAT8 v = a*a;
453     FLOAT8 w = 0.0;   
454     u -= athFloor;                                  // undo scaling
455     if ( v > 1E-20 ) w = 1. + 10. * log10(v) / o;
456     if ( w < 0  )    w = 0.; 
457     u *= w; 
458     u += athFloor + o-p;                            // redo scaling
459
460     return pow( 10., 0.1*u );
461 }
462
463
464 /*************************************************************************/
465 /*            calc_xmin                                                  */
466 /*************************************************************************/
467
468 /*
469   Calculate the allowed distortion for each scalefactor band,
470   as determined by the psychoacoustic model.
471   xmin(sb) = ratio(sb) * en(sb) / bw(sb)
472
473   returns number of sfb's with energy > ATH
474 */
475 int calc_xmin( 
476         lame_global_flags *gfp,
477         const FLOAT8                xr [576],
478         const III_psy_ratio * const ratio,
479         const gr_info       * const cod_info, 
480               III_psy_xmin  * const l3_xmin ) 
481 {
482     lame_internal_flags *gfc = gfp->internal_flags;
483     int sfb,j,start, end, bw,l, b, ath_over=0;
484     FLOAT8 en0, xmin, ener, tmpATH;
485     ATH_t * ATH = gfc->ATH;
486
487     if (cod_info->block_type == SHORT_TYPE) {
488
489         for ( j = 0, sfb = 0; sfb < SBMAX_s; sfb++ ) {
490             if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh )
491                 tmpATH = athAdjust( ATH->adjust, ATH->s[sfb], ATH->floor );
492             else
493                 tmpATH = ATH->adjust * ATH->s[sfb];
494             start = gfc->scalefac_band.s[ sfb ];
495             end   = gfc->scalefac_band.s[ sfb + 1 ];
496             bw = end - start;
497
498             for ( b = 0; b < 3; b++ ) {
499                 for (en0 = 0.0, l = start; l < end; l++) {
500                     ener = xr[j++];
501                     ener = ener * ener;
502                     en0 += ener;
503                 }
504                 en0 /= bw;
505       
506                 if (gfp->ATHonly || gfp->ATHshort) {
507                     xmin = tmpATH;
508                 } 
509                 else {
510                     xmin = ratio->en.s[sfb][b];
511                     if (xmin > 0.0)
512                         xmin = en0 * ratio->thm.s[sfb][b] * gfc->masking_lower / xmin;
513                     if (xmin < tmpATH) 
514                         xmin = tmpATH;
515                 }
516                 xmin *= bw;
517
518                 if (gfc->nsPsy.use) {
519                     if      (sfb <=  5) xmin *= gfc->nsPsy.bass;
520                     else if (sfb <= 10) xmin *= gfc->nsPsy.alto;
521                     else                xmin *= gfc->nsPsy.treble;
522                     if ((gfp->VBR == vbr_off || gfp->VBR == vbr_abr) && gfp->quality <= 1)
523                         xmin *= 0.001;
524                 }
525                 l3_xmin->s[sfb][b] = xmin;
526                 if (en0 > tmpATH) ath_over++;
527             }   /* b */
528         }   /* sfb */
529
530         if (gfp->useTemporal) {
531             for (sfb = 0; sfb < SBMAX_s; sfb++ ) {
532                 for ( b = 1; b < 3; b++ ) {
533                     xmin = l3_xmin->s[sfb][b] * (1.0 - gfc->decay)
534                          + l3_xmin->s[sfb][b-1] * gfc->decay;
535                     if (l3_xmin->s[sfb][b] < xmin)
536                         l3_xmin->s[sfb][b] = xmin;
537                 }
538             }
539         }
540
541     }   /* end of short block case */
542     else {
543         
544         for ( sfb = 0; sfb < SBMAX_l; sfb++ ){
545             if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh )
546                 tmpATH = athAdjust( ATH->adjust, ATH->l[sfb], ATH->floor );
547             else
548                 tmpATH = ATH->adjust * ATH->l[sfb];
549             start = gfc->scalefac_band.l[ sfb ];
550             end   = gfc->scalefac_band.l[ sfb+1 ];
551             bw = end - start;
552     
553             for (en0 = 0.0, l = start; l < end; l++ ) {
554                 ener = xr[l] * xr[l];
555                 en0 += ener;
556             }
557             /* why is it different from short blocks <?> */
558             if ( !gfc->nsPsy.use ) en0 /= bw;   
559     
560             if (gfp->ATHonly) {
561                 xmin = tmpATH;
562             }
563             else {
564                 xmin = ratio->en.l[sfb];
565                 if (xmin > 0.0)
566                     xmin = en0 * ratio->thm.l[sfb] * gfc->masking_lower / xmin;
567                 if (xmin < tmpATH) 
568                     xmin = tmpATH;
569             }
570             /* why is it different from short blocks <?> */
571             if ( !gfc->nsPsy.use ) {
572                 xmin *= bw;
573             }
574             else {
575                 if      (sfb <=  6) xmin *= gfc->nsPsy.bass;
576                 else if (sfb <= 13) xmin *= gfc->nsPsy.alto;
577                 else if (sfb <= 20) xmin *= gfc->nsPsy.treble;
578                 else                xmin *= gfc->nsPsy.sfb21;
579                 if ((gfp->VBR == vbr_off || gfp->VBR == vbr_abr) && gfp->quality <= 1)
580                     xmin *= 0.001;
581             }
582             l3_xmin->l[sfb] = xmin;
583             if (en0 > tmpATH) ath_over++;
584         }   /* sfb */
585         
586     }   /* end of long block case */
587
588     return ath_over;
589 }
590
591
592
593
594
595
596
597 /*************************************************************************/
598 /*            calc_noise                                                 */
599 /*************************************************************************/
600
601 // -oo dB  =>  -1.00
602 // - 6 dB  =>  -0.97
603 // - 3 dB  =>  -0.80
604 // - 2 dB  =>  -0.64
605 // - 1 dB  =>  -0.38
606 //   0 dB  =>   0.00
607 // + 1 dB  =>  +0.49
608 // + 2 dB  =>  +1.06
609 // + 3 dB  =>  +1.68
610 // + 6 dB  =>  +3.69
611 // +10 dB  =>  +6.45
612
613 double penalties ( double noise )
614 {
615     return log ( 0.368 + 0.632 * noise * noise * noise );
616 }
617
618 /*  mt 5/99:  Function: Improved calc_noise for a single channel   */
619
620 int  calc_noise( 
621         const lame_internal_flags           * const gfc,
622         const FLOAT8                    xr [576],
623         const int                       ix [576],
624         const gr_info           * const cod_info,
625         const III_psy_xmin      * const l3_xmin, 
626         const III_scalefac_t    * const scalefac,
627               III_psy_xmin      * xfsf,
628               calc_noise_result * const res )
629 {
630   int sfb,start, end, j,l, i, over=0;
631   FLOAT8 sum;
632   
633   int count=0;
634   FLOAT8 noise,noise_db;
635   FLOAT8 over_noise_db = 0;
636   FLOAT8 tot_noise_db  = 0;     /*    0 dB relative to masking */
637   FLOAT8 max_noise  = 1E-20; /* -200 dB relative to masking */
638   double klemm_noise = 1E-37;
639
640   if (cod_info->block_type == SHORT_TYPE) {
641     int max_index = gfc->sfb21_extra ? SBMAX_s : SBPSY_s;
642
643     for ( j=0, sfb = 0; sfb < max_index; sfb++ ) {
644          start = gfc->scalefac_band.s[ sfb ];
645          end   = gfc->scalefac_band.s[ sfb+1 ];
646          for ( i = 0; i < 3; i++ ) {
647             FLOAT8 step;
648             int s;
649
650             s = (scalefac->s[sfb][i] << (cod_info->scalefac_scale + 1))
651                 + cod_info->subblock_gain[i] * 8;
652             s = cod_info->global_gain - s;
653
654             assert(s<Q_MAX);
655             assert(s>=0);
656             step = POW20(s);
657             sum = 0.0;
658             l = start;
659             do {
660               FLOAT8 temp;
661               temp = pow43[ix[j]];
662               temp *= step;
663               temp -= fabs(xr[j]);
664               ++j;
665               sum += temp * temp;
666               l++;
667             } while (l < end);
668             noise = xfsf->s[sfb][i]  = sum / l3_xmin->s[sfb][i];
669
670             max_noise    = Max(max_noise,noise);
671             klemm_noise += penalties (noise);
672
673             noise_db=10*log10(Max(noise,1E-20));
674             tot_noise_db += noise_db;
675
676             if (noise > 1) {
677                 over++;
678                 over_noise_db += noise_db;
679             }
680             count++;        
681         }
682     }
683   }else{ /* cod_info->block_type == SHORT_TYPE */
684     int max_index = gfc->sfb21_extra ? SBMAX_l : SBPSY_l;
685     
686     for ( sfb = 0; sfb < max_index; sfb++ ) {
687         FLOAT8 step;
688         int s = scalefac->l[sfb];
689
690         if (cod_info->preflag)
691             s += pretab[sfb];
692
693         s = cod_info->global_gain - (s << (cod_info->scalefac_scale + 1));
694         assert(s<Q_MAX);
695         assert(s>=0);
696         step = POW20(s);
697
698         start = gfc->scalefac_band.l[ sfb ];
699         end   = gfc->scalefac_band.l[ sfb+1 ];
700
701         for ( sum = 0.0, l = start; l < end; l++ ) {
702           FLOAT8 temp;
703           temp = fabs(xr[l]) - pow43[ix[l]] * step;
704           sum += temp * temp;
705         }
706         noise = xfsf->l[sfb] = sum / l3_xmin->l[sfb];
707         max_noise=Max(max_noise,noise);
708         klemm_noise += penalties (noise);
709
710         noise_db=10*log10(Max(noise,1E-20));
711         /* multiplying here is adding in dB, but can overflow */
712         //tot_noise *= Max(noise, 1E-20);
713         tot_noise_db += noise_db;
714
715         if (noise > 1) {
716           over++;
717           /* multiplying here is adding in dB -but can overflow */
718           //over_noise *= noise;
719           over_noise_db += noise_db;
720         }
721
722         count++;
723     }
724   } /* cod_info->block_type == SHORT_TYPE */
725
726   /* normalization at this point by "count" is not necessary, since
727    * the values are only used to compare with previous values */
728   res->over_count = over;
729
730   /* convert to db. DO NOT CHANGE THESE */
731   /* tot_noise = is really the average over each sfb of: 
732      [noise(db) - allowed_noise(db)]
733
734      and over_noise is the same average, only over only the
735      bands with noise > allowed_noise.  
736
737   */
738
739   //res->tot_noise   = 10.*log10(Max(1e-20,tot_noise )); 
740   //res->over_noise  = 10.*log10(Max(1e+00,over_noise)); 
741   res->tot_noise   = tot_noise_db;
742   res->over_noise  = over_noise_db;
743   res->max_noise   = 10.*log10(Max(1e-20,max_noise ));
744   res->klemm_noise = 10.*log10(Max(1e-20,klemm_noise));
745
746   return over;
747 }
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762 /************************************************************************
763  *
764  *  set_pinfo()
765  *
766  *  updates plotting data    
767  *
768  *  Mark Taylor 2000-??-??                
769  *
770  *  Robert Hegemann: moved noise/distortion calc into it                          
771  *
772  ************************************************************************/
773
774 static
775 void set_pinfo (
776         lame_global_flags *gfp,
777         const gr_info        * const cod_info,
778         const III_psy_ratio  * const ratio, 
779         const III_scalefac_t * const scalefac,
780         const FLOAT8                 xr[576],
781         const int                    l3_enc[576],        
782         const int                    gr,
783         const int                    ch )
784 {
785     lame_internal_flags *gfc=gfp->internal_flags;
786     int sfb;
787     int j,i,l,start,end,bw;
788     FLOAT8 en0,en1;
789     FLOAT ifqstep = ( cod_info->scalefac_scale == 0 ) ? .5 : 1.0;
790
791
792     III_psy_xmin l3_xmin;
793     calc_noise_result noise;
794     III_psy_xmin xfsf;
795
796     calc_xmin (gfp,xr, ratio, cod_info, &l3_xmin);
797
798     calc_noise (gfc, xr, l3_enc, cod_info, &l3_xmin, scalefac, &xfsf, &noise);
799
800     if (cod_info->block_type == SHORT_TYPE) {
801         for (j=0, sfb = 0; sfb < SBMAX_s; sfb++ )  {
802             start = gfc->scalefac_band.s[ sfb ];
803             end   = gfc->scalefac_band.s[ sfb + 1 ];
804             bw = end - start;
805             for ( i = 0; i < 3; i++ ) {
806                 for ( en0 = 0.0, l = start; l < end; l++ ) {
807                     en0 += xr[j] * xr[j];
808                     ++j;
809                 }
810                 en0=Max(en0/bw,1e-20);
811
812
813 #if 0
814 {
815     double tot1,tot2;
816     if (sfb<SBMAX_s-1) {
817         if (sfb==0) {
818             tot1=0;
819             tot2=0;
820         }
821         tot1 += en0;
822         tot2 += ratio->en.s[sfb][i];
823
824         DEBUGF("%i %i sfb=%i mdct=%f fft=%f  fft-mdct=%f db \n",
825                 gr,ch,sfb,
826                 10*log10(Max(1e-25,ratio->en.s[sfb][i])),
827                 10*log10(Max(1e-25,en0)),
828                 10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.s[sfb][i]))));
829
830         if (sfb==SBMAX_s-2) {
831             DEBUGF("%i %i toti %f %f ratio=%f db \n",gr,ch,
832                     10*log10(Max(1e-25,tot2)),
833                     10*log(Max(1e-25,tot1)),
834                     10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
835
836         }
837     }
838     /*
839         masking: multiplied by en0/fft_energy
840         average seems to be about -135db.
841      */
842 }
843 #endif
844
845
846                 /* convert to MDCT units */
847                 en1=1e15;  /* scaling so it shows up on FFT plot */
848                 gfc->pinfo->xfsf_s[gr][ch][3*sfb+i] 
849                     = en1*xfsf.s[sfb][i]*l3_xmin.s[sfb][i]/bw;
850                 gfc->pinfo->en_s[gr][ch][3*sfb+i] = en1*en0;
851
852                 if (ratio->en.s[sfb][i]>0)
853                     en0 = en0/ratio->en.s[sfb][i];
854                 else
855                     en0=0;
856                 if (gfp->ATHonly || gfp->ATHshort)
857                     en0=0;
858
859                 gfc->pinfo->thr_s[gr][ch][3*sfb+i] = 
860                         en1*Max(en0*ratio->thm.s[sfb][i],gfc->ATH->s[sfb]);
861
862  
863                 /* there is no scalefactor bands >= SBPSY_s */
864                 if (sfb < SBPSY_s) {
865                     gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=
866                                             -ifqstep*scalefac->s[sfb][i];
867                 } else {
868                     gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=0;
869                 }
870                 gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i] -=
871                                              2*cod_info->subblock_gain[i];
872             }
873         }
874     } else {
875         for ( sfb = 0; sfb < SBMAX_l; sfb++ )   {
876             start = gfc->scalefac_band.l[ sfb ];
877             end   = gfc->scalefac_band.l[ sfb+1 ];
878             bw = end - start;
879             for ( en0 = 0.0, l = start; l < end; l++ ) 
880                 en0 += xr[l] * xr[l];
881             if (!gfc->nsPsy.use) en0/=bw;
882       /*
883     DEBUGF("diff  = %f \n",10*log10(Max(ratio[gr][ch].en.l[sfb],1e-20))
884                             -(10*log10(en0)+150));
885        */
886
887 #if 0
888  {
889     double tot1,tot2;
890     if (sfb==0) {
891         tot1=0;
892         tot2=0;
893     }
894     tot1 += en0;
895     tot2 += ratio->en.l[sfb];
896
897
898     DEBUGF("%i %i sfb=%i mdct=%f fft=%f  fft-mdct=%f db \n",
899             gr,ch,sfb,
900             10*log10(Max(1e-25,ratio->en.l[sfb])),
901             10*log10(Max(1e-25,en0)),
902             10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.l[sfb]))));
903
904     if (sfb==SBMAX_l-1) {
905         DEBUGF("%i %i toti %f %f ratio=%f db \n",
906             gr,ch,
907             10*log10(Max(1e-25,tot2)),
908             10*log(Max(1e-25,tot1)),
909             10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
910     }
911     /*
912         masking: multiplied by en0/fft_energy
913         average seems to be about -147db.
914      */
915 }
916 #endif
917
918
919             /* convert to MDCT units */
920             en1=1e15;  /* scaling so it shows up on FFT plot */
921             gfc->pinfo->xfsf[gr][ch][sfb] =  en1*xfsf.l[sfb]*l3_xmin.l[sfb]/bw;
922             gfc->pinfo->en[gr][ch][sfb] = en1*en0;
923             if (ratio->en.l[sfb]>0)
924                 en0 = en0/ratio->en.l[sfb];
925             else
926                 en0=0;
927             if (gfp->ATHonly)
928                 en0=0;
929             gfc->pinfo->thr[gr][ch][sfb] =
930                              en1*Max(en0*ratio->thm.l[sfb],gfc->ATH->l[sfb]);
931
932             /* there is no scalefactor bands >= SBPSY_l */
933             if (sfb<SBPSY_l) {
934                 if (scalefac->l[sfb]<0)  /* scfsi! */
935                     gfc->pinfo->LAMEsfb[gr][ch][sfb] =
936                                             gfc->pinfo->LAMEsfb[0][ch][sfb];
937                 else
938                     gfc->pinfo->LAMEsfb[gr][ch][sfb] = -ifqstep*scalefac->l[sfb];
939             }else{
940                 gfc->pinfo->LAMEsfb[gr][ch][sfb] = 0;
941             }
942
943             if (cod_info->preflag && sfb>=11) 
944                 gfc->pinfo->LAMEsfb[gr][ch][sfb] -= ifqstep*pretab[sfb];
945         } /* for sfb */
946     } /* block type long */
947     gfc->pinfo->LAMEqss     [gr][ch] = cod_info->global_gain;
948     gfc->pinfo->LAMEmainbits[gr][ch] = cod_info->part2_3_length;
949     gfc->pinfo->LAMEsfbits  [gr][ch] = cod_info->part2_length;
950
951     gfc->pinfo->over      [gr][ch] = noise.over_count;
952     gfc->pinfo->max_noise [gr][ch] = noise.max_noise;
953     gfc->pinfo->over_noise[gr][ch] = noise.over_noise;
954     gfc->pinfo->tot_noise [gr][ch] = noise.tot_noise;
955 }
956
957
958 /************************************************************************
959  *
960  *  set_frame_pinfo()
961  *
962  *  updates plotting data for a whole frame  
963  *
964  *  Robert Hegemann 2000-10-21                          
965  *
966  ************************************************************************/
967
968 void set_frame_pinfo( 
969         lame_global_flags *gfp,
970         FLOAT8          xr       [2][2][576],
971         III_psy_ratio   ratio    [2][2],  
972         int             l3_enc   [2][2][576],
973         III_scalefac_t  scalefac [2][2] )
974 {
975     lame_internal_flags *gfc=gfp->internal_flags;
976     unsigned int          sfb;
977         int                   ch;
978         int                   gr;
979     int                   act_l3enc[576];
980     III_scalefac_t        act_scalefac [2];
981     int scsfi[2] = {0,0};
982     
983     
984     gfc->masking_lower = 1.0;
985     
986     /* reconstruct the scalefactors in case SCSFI was used 
987      */
988     for (ch = 0; ch < gfc->channels_out; ch ++) {
989         for (sfb = 0; sfb < SBMAX_l; sfb ++) {
990             if (scalefac[1][ch].l[sfb] == -1) {/* scfsi */
991                 act_scalefac[ch].l[sfb] = scalefac[0][ch].l[sfb];
992                 scsfi[ch] = 1;
993             } else {
994                 act_scalefac[ch].l[sfb] = scalefac[1][ch].l[sfb];
995             }
996         }
997     }
998     
999     /* for every granule and channel patch l3_enc and set info
1000      */
1001     for (gr = 0; gr < gfc->mode_gr; gr ++) {
1002         for (ch = 0; ch < gfc->channels_out; ch ++) {
1003             int i;
1004             gr_info *cod_info = &gfc->l3_side.gr[gr].ch[ch].tt;
1005             
1006             /* revert back the sign of l3enc */
1007             for ( i = 0; i < 576; i++) {
1008                 if (xr[gr][ch][i] < 0) 
1009                     act_l3enc[i] = -l3_enc[gr][ch][i];
1010                 else
1011                     act_l3enc[i] = +l3_enc[gr][ch][i];
1012             }
1013             if (gr == 1 && scsfi[ch]) 
1014                 set_pinfo (gfp, cod_info, &ratio[gr][ch], &act_scalefac[ch],
1015                         xr[gr][ch], act_l3enc, gr, ch);                    
1016             else
1017                 set_pinfo (gfp, cod_info, &ratio[gr][ch], &scalefac[gr][ch],
1018                         xr[gr][ch], act_l3enc, gr, ch);                    
1019         } /* for ch */
1020     }    /* for gr */
1021 }
1022         
1023
1024
1025
1026 /*********************************************************************
1027  * nonlinear quantization of xr 
1028  * More accurate formula than the ISO formula.  Takes into account
1029  * the fact that we are quantizing xr -> ix, but we want ix^4/3 to be 
1030  * as close as possible to x^4/3.  (taking the nearest int would mean
1031  * ix is as close as possible to xr, which is different.)
1032  * From Segher Boessenkool <segher@eastsite.nl>  11/1999
1033  * ASM optimization from 
1034  *    Mathew Hendry <scampi@dial.pipex.com> 11/1999
1035  *    Acy Stapp <AStapp@austin.rr.com> 11/1999
1036  *    Takehiro Tominaga <tominaga@isoternet.org> 11/1999
1037  * 9/00: ASM code removed in favor of IEEE754 hack.  If you need
1038  * the ASM code, check CVS circa Aug 2000.  
1039  *********************************************************************/
1040
1041
1042 #ifdef TAKEHIRO_IEEE754_HACK
1043
1044 typedef union {
1045     float f;
1046     int i;
1047 } fi_union;
1048
1049 #define MAGIC_FLOAT (65536*(128))
1050 #define MAGIC_INT 0x4b000000
1051
1052 void quantize_xrpow(const FLOAT8 *xp, int *pi, FLOAT8 istep)
1053 {
1054     /* quantize on xr^(3/4) instead of xr */
1055     int j;
1056     fi_union *fi;
1057
1058     fi = (fi_union *)pi;
1059     for (j = 576 / 4 - 1; j >= 0; --j) {
1060         double x0 = istep * xp[0];
1061         double x1 = istep * xp[1];
1062         double x2 = istep * xp[2];
1063         double x3 = istep * xp[3];
1064
1065         x0 += MAGIC_FLOAT; fi[0].f = x0;
1066         x1 += MAGIC_FLOAT; fi[1].f = x1;
1067         x2 += MAGIC_FLOAT; fi[2].f = x2;
1068         x3 += MAGIC_FLOAT; fi[3].f = x3;
1069
1070         fi[0].f = x0 + (adj43asm - MAGIC_INT)[fi[0].i];
1071         fi[1].f = x1 + (adj43asm - MAGIC_INT)[fi[1].i];
1072         fi[2].f = x2 + (adj43asm - MAGIC_INT)[fi[2].i];
1073         fi[3].f = x3 + (adj43asm - MAGIC_INT)[fi[3].i];
1074
1075         fi[0].i -= MAGIC_INT;
1076         fi[1].i -= MAGIC_INT;
1077         fi[2].i -= MAGIC_INT;
1078         fi[3].i -= MAGIC_INT;
1079         fi += 4;
1080         xp += 4;
1081     }
1082 }
1083
1084 #  define ROUNDFAC -0.0946
1085 void quantize_xrpow_ISO(const FLOAT8 *xp, int *pi, FLOAT8 istep)
1086 {
1087     /* quantize on xr^(3/4) instead of xr */
1088     int j;
1089     fi_union *fi;
1090
1091     fi = (fi_union *)pi;
1092     for (j=576/4 - 1;j>=0;j--) {
1093         fi[0].f = istep * xp[0] + (ROUNDFAC + MAGIC_FLOAT);
1094         fi[1].f = istep * xp[1] + (ROUNDFAC + MAGIC_FLOAT);
1095         fi[2].f = istep * xp[2] + (ROUNDFAC + MAGIC_FLOAT);
1096         fi[3].f = istep * xp[3] + (ROUNDFAC + MAGIC_FLOAT);
1097
1098         fi[0].i -= MAGIC_INT;
1099         fi[1].i -= MAGIC_INT;
1100         fi[2].i -= MAGIC_INT;
1101         fi[3].i -= MAGIC_INT;
1102         fi+=4;
1103         xp+=4;
1104     }
1105 }
1106
1107 #else
1108
1109 /*********************************************************************
1110  * XRPOW_FTOI is a macro to convert floats to ints.  
1111  * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x]
1112  *                                         ROUNDFAC= -0.0946
1113  *
1114  * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]   
1115  *                                   ROUNDFAC=0.4054
1116  *
1117  * Note: using floor() or (int) is extermely slow. On machines where
1118  * the TAKEHIRO_IEEE754_HACK code above does not work, it is worthwile
1119  * to write some ASM for XRPOW_FTOI().  
1120  *********************************************************************/
1121 #define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
1122 #define QUANTFAC(rx)  adj43[rx]
1123 #define ROUNDFAC 0.4054
1124
1125
1126 void quantize_xrpow(const FLOAT8 *xr, int *ix, FLOAT8 istep) {
1127     /* quantize on xr^(3/4) instead of xr */
1128     /* from Wilfried.Behne@t-online.de.  Reported to be 2x faster than 
1129        the above code (when not using ASM) on PowerPC */
1130     int j;
1131
1132     for ( j = 576/8; j > 0; --j) {
1133         FLOAT8  x1, x2, x3, x4, x5, x6, x7, x8;
1134         int     rx1, rx2, rx3, rx4, rx5, rx6, rx7, rx8;
1135         x1 = *xr++ * istep;
1136         x2 = *xr++ * istep;
1137         XRPOW_FTOI(x1, rx1);
1138         x3 = *xr++ * istep;
1139         XRPOW_FTOI(x2, rx2);
1140         x4 = *xr++ * istep;
1141         XRPOW_FTOI(x3, rx3);
1142         x5 = *xr++ * istep;
1143         XRPOW_FTOI(x4, rx4);
1144         x6 = *xr++ * istep;
1145         XRPOW_FTOI(x5, rx5);
1146         x7 = *xr++ * istep;
1147         XRPOW_FTOI(x6, rx6);
1148         x8 = *xr++ * istep;
1149         XRPOW_FTOI(x7, rx7);
1150         x1 += QUANTFAC(rx1);
1151         XRPOW_FTOI(x8, rx8);
1152         x2 += QUANTFAC(rx2);
1153         XRPOW_FTOI(x1,*ix++);
1154         x3 += QUANTFAC(rx3);
1155         XRPOW_FTOI(x2,*ix++);
1156         x4 += QUANTFAC(rx4);
1157         XRPOW_FTOI(x3,*ix++);
1158         x5 += QUANTFAC(rx5);
1159         XRPOW_FTOI(x4,*ix++);
1160         x6 += QUANTFAC(rx6);
1161         XRPOW_FTOI(x5,*ix++);
1162         x7 += QUANTFAC(rx7);
1163         XRPOW_FTOI(x6,*ix++);
1164         x8 += QUANTFAC(rx8);
1165         XRPOW_FTOI(x7,*ix++);
1166         XRPOW_FTOI(x8,*ix++);
1167     }
1168 }
1169
1170
1171
1172
1173
1174
1175 void quantize_xrpow_ISO( const FLOAT8 *xr, int *ix, FLOAT8 istep )
1176 {
1177     /* quantize on xr^(3/4) instead of xr */
1178     const FLOAT8 compareval0 = (1.0 - 0.4054)/istep;
1179     int j;
1180     /* depending on architecture, it may be worth calculating a few more
1181        compareval's.
1182
1183        eg.  compareval1 = (2.0 - 0.4054/istep);
1184        .. and then after the first compare do this ...
1185        if compareval1>*xr then ix = 1;
1186
1187        On a pentium166, it's only worth doing the one compare (as done here),
1188        as the second compare becomes more expensive than just calculating
1189        the value. Architectures with slow FP operations may want to add some
1190        more comparevals. try it and send your diffs statistically speaking
1191
1192        73% of all xr*istep values give ix=0
1193        16% will give 1
1194        4%  will give 2
1195     */
1196     for (j=576;j>0;j--) {
1197         if (compareval0 > *xr) {
1198             *(ix++) = 0;
1199             xr++;
1200         } else {
1201             /*    *(ix++) = (int)( istep*(*(xr++))  + 0.4054); */
1202             XRPOW_FTOI(  istep*(*(xr++))  + ROUNDFAC , *(ix++) );
1203         }
1204     }
1205 }
1206
1207 #endif