home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 13 / AACD13.ISO / AACD / Sound / LAME / src / quantize-pvt.c < prev    next >
C/C++ Source or Header  |  2000-08-06  |  40KB  |  1,525 lines
  1. #include <assert.h>
  2. #include "util.h"
  3. #include "gtkanal.h"
  4. #include "tables.h"
  5. #include "reservoir.h"
  6. #include "quantize-pvt.h"
  7.  
  8. /* some problems found with -O2 and above, gcc 2.95 */
  9. #if (defined(__GNUC__) && defined(__i386__))
  10. #undef TAKEHIRO_IEEE754_HACK
  11. #endif
  12.  
  13. #define NSTHRE 4  /* tuned by hearing tests */
  14.  
  15. const int slen1_tab[16] = { 0, 0, 0, 0, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4 };
  16. const int slen2_tab[16] = { 0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3, 1, 2, 3, 2, 3 };
  17.  
  18.  
  19. /*
  20.   The following table is used to implement the scalefactor
  21.   partitioning for MPEG2 as described in section
  22.   2.4.3.2 of the IS. The indexing corresponds to the
  23.   way the tables are presented in the IS:
  24.  
  25.   [table_number][row_in_table][column of nr_of_sfb]
  26. */
  27. unsigned int nr_of_sfb_block[6][3][4] =
  28. {
  29.   {
  30.     {6, 5, 5, 5},
  31.     {9, 9, 9, 9},
  32.     {6, 9, 9, 9}
  33.   },
  34.   {
  35.     {6, 5, 7, 3},
  36.     {9, 9, 12, 6},
  37.     {6, 9, 12, 6}
  38.   },
  39.   {
  40.     {11, 10, 0, 0},
  41.     {18, 18, 0, 0},
  42.     {15,18,0,0}
  43.   },
  44.   {
  45.     {7, 7, 7, 0},
  46.     {12, 12, 12, 0},
  47.     {6, 15, 12, 0}
  48.   },
  49.   {
  50.     {6, 6, 6, 3},
  51.     {12, 9, 9, 6},
  52.     {6, 12, 9, 6}
  53.   },
  54.   {
  55.     {8, 8, 5, 0},
  56.     {15,12,9,0},
  57.     {6,18,9,0}
  58.   }
  59. };
  60.  
  61.  
  62. /* Table B.6: layer3 preemphasis */
  63. int  pretab[SBMAX_l] =
  64. {
  65.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  66.     1, 1, 1, 1, 2, 2, 3, 3, 3, 2, 0
  67. };
  68.  
  69. /*
  70.   Here are MPEG1 Table B.8 and MPEG2 Table B.1
  71.   -- Layer III scalefactor bands. 
  72.   Index into this using a method such as:
  73.     idx  = fr_ps->header->sampling_frequency
  74.            + (fr_ps->header->version * 3)
  75. */
  76.  
  77.  
  78. const scalefac_struct sfBandIndex[9] =
  79. {
  80.   { /* Table B.2.b: 22.05 kHz */
  81.     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
  82.     {0,4,8,12,18,24,32,42,56,74,100,132,174,192}
  83.   },
  84.   { /* Table B.2.c: 24 kHz */                 /* docs: 332. mpg123(broken): 330 */
  85.     {0,6,12,18,24,30,36,44,54,66,80,96,114,136,162,194,232,278, 332, 394,464,540,576},
  86.     {0,4,8,12,18,26,36,48,62,80,104,136,180,192}
  87.   },
  88.   { /* Table B.2.a: 16 kHz */
  89.     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
  90.     {0,4,8,12,18,26,36,48,62,80,104,134,174,192}
  91.   },
  92.   { /* Table B.8.b: 44.1 kHz */
  93.     {0,4,8,12,16,20,24,30,36,44,52,62,74,90,110,134,162,196,238,288,342,418,576},
  94.     {0,4,8,12,16,22,30,40,52,66,84,106,136,192}
  95.   },
  96.   { /* Table B.8.c: 48 kHz */
  97.     {0,4,8,12,16,20,24,30,36,42,50,60,72,88,106,128,156,190,230,276,330,384,576},
  98.     {0,4,8,12,16,22,28,38,50,64,80,100,126,192}
  99.   },
  100.   { /* Table B.8.a: 32 kHz */
  101.     {0,4,8,12,16,20,24,30,36,44,54,66,82,102,126,156,194,240,296,364,448,550,576},
  102.     {0,4,8,12,16,22,30,42,58,78,104,138,180,192}
  103.   },
  104.   { /* MPEG-2.5 11.025 kHz */
  105.     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
  106.     {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}
  107.   },
  108.   { /* MPEG-2.5 12 kHz */
  109.     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
  110.     {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}
  111.   },
  112.   { /* MPEG-2.5 8 kHz */
  113.     {0,12,24,36,48,60,72,88,108,132,160,192,232,280,336,400,476,566,568,570,572,574,576},
  114.     {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}
  115.   }
  116. };
  117.  
  118.  
  119.  
  120. FLOAT8 pow20[Q_MAX];
  121. FLOAT8 ipow20[Q_MAX];
  122. FLOAT8 pow43[PRECALC_SIZE];
  123. /* initialized in first call to iteration_init */
  124. FLOAT8 adj43[PRECALC_SIZE];
  125. FLOAT8 adj43asm[PRECALC_SIZE];
  126.  
  127.  
  128. /************************************************************************/
  129. /*  initialization for iteration_loop */
  130. /************************************************************************/
  131. void
  132. iteration_init( lame_global_flags *gfp,III_side_info_t *l3_side, int l3_enc[2][2][576])
  133. {
  134.   lame_internal_flags *gfc=gfp->internal_flags;
  135.   gr_info *cod_info;
  136.   int ch, gr, i;
  137.  
  138.   if ( gfc->iteration_init_init==0 ) {
  139.     gfc->iteration_init_init=1;
  140.  
  141.     l3_side->main_data_begin = 0;
  142.     compute_ath(gfp,gfc->ATH_l,gfc->ATH_s);
  143.  
  144.     for(i=0;i<PRECALC_SIZE;i++)
  145.         pow43[i] = pow((FLOAT8)i, 4.0/3.0);
  146.  
  147.     for (i = 0; i < PRECALC_SIZE-1; i++)
  148.     adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
  149.     adj43[i] = 0.5;
  150.  
  151.     adj43asm[0] = 0.0;
  152.     for (i = 1; i < PRECALC_SIZE; i++)
  153.       adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]),0.75);
  154.  
  155.     for (i = 0; i < Q_MAX; i++) {
  156.     ipow20[i] = pow(2.0, (double)(i - 210) * -0.1875);
  157.     pow20[i] = pow(2.0, (double)(i - 210) * 0.25);
  158.     }
  159.   }
  160.  
  161.  
  162.  
  163.   
  164.   /* some intializations. */
  165.   for ( gr = 0; gr < gfc->mode_gr; gr++ ){
  166.     for ( ch = 0; ch < gfc->stereo; ch++ ){
  167.       cod_info = (gr_info *) &(l3_side->gr[gr].ch[ch]);
  168.  
  169.       if (cod_info->block_type == SHORT_TYPE)
  170.         {
  171.       cod_info->sfb_lmax = 0; /* No sb*/
  172.       cod_info->sfb_smax = 0;
  173.         }
  174.       else
  175.     {
  176.       /* MPEG 1 doesnt use last scalefactor band */
  177.       cod_info->sfb_lmax = SBPSY_l;
  178.       cod_info->sfb_smax = SBPSY_s;    /* No sb */
  179.     }
  180.  
  181.     }
  182.   }
  183.  
  184.   huffman_init(gfp);
  185. }
  186.  
  187.  
  188.  
  189.  
  190.  
  191. /* 
  192. compute the ATH for each scalefactor band 
  193. cd range:  0..96db
  194.  
  195. Input:  3.3kHz signal  32767 amplitude  (3.3kHz is where ATH is smallest = -5db)
  196. longblocks:  sfb=12   en0/bw=-11db    max_en0 = 1.3db
  197. shortblocks: sfb=5           -9db              0db
  198.  
  199. Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
  200. longblocks:  amp=1      sfb=12   en0/bw=-103 db      max_en0 = -92db
  201.             amp=32767   sfb=12           -12 db                 -1.4db 
  202.  
  203. Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
  204. shortblocks: amp=1      sfb=5   en0/bw= -99                    -86 
  205.             amp=32767   sfb=5           -9  db                  4db 
  206.  
  207.  
  208. MAX energy of largest wave at 3.3kHz = 1db
  209. AVE energy of largest wave at 3.3kHz = -11db
  210. Let's take AVE:  -11db = maximum signal in sfb=12.  
  211. Dynamic range of CD: 96db.  Therefor energy of smallest audible wave 
  212. in sfb=12  = -11  - 96 = -107db = ATH at 3.3kHz.  
  213.  
  214. ATH formula for this wave: -5db.  To adjust to LAME scaling, we need
  215. ATH = ATH_formula  - 103  (db)
  216. ATH = ATH * 2.5e-10      (ener)
  217.  
  218. */
  219.  
  220. FLOAT8 ATHmdct(lame_global_flags *gfp,FLOAT8 f)
  221. {
  222.   lame_internal_flags *gfc=gfp->internal_flags;
  223.   FLOAT8 ath;
  224.   
  225.   ath = ATHformula(f);
  226.       
  227.   /* convert to energy */
  228.   ath -= 114;    /* MDCT scaling.  From tests by macik and MUS420 code */
  229.   ath -= gfp->ATHlower;
  230.  
  231.   /* purpose of RH_QUALITY_CONTROL:
  232.    * at higher quality lower ATH masking abilities   => needs more bits
  233.    * at lower quality increase ATH masking abilities => needs less bits
  234.    * works together with adjusted masking lowering of GPSYCHO thresholds
  235.    * (Robert.Hegemann@gmx.de 2000-01-30)
  236.    */
  237.   if (gfp->VBR!=vbr_off) 
  238.     {
  239.       ath -= gfc->ATH_vbrlower;
  240.       ath = Min(gfp->VBR_q-62,ath);
  241.     }
  242.     
  243.   ath = pow( 10.0, ath/10.0 );
  244.   return ath;
  245. }
  246.  
  247.  
  248. void compute_ath(lame_global_flags *gfp,FLOAT8 ATH_l[],FLOAT8 ATH_s[])
  249. {
  250.   lame_internal_flags *gfc=gfp->internal_flags;
  251.   int sfb,i,start,end;
  252.   FLOAT8 ATH_f;
  253.   FLOAT8 samp_freq = gfp->out_samplerate/1000.0;
  254.  
  255.   for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
  256.     start = gfc->scalefac_band.l[ sfb ];
  257.     end   = gfc->scalefac_band.l[ sfb+1 ];
  258.     ATH_l[sfb]=1e99;
  259.     for (i=start ; i < end; i++) {
  260.       FLOAT8 freq = samp_freq*i/(2*576);
  261.       assert( freq < 25 );
  262.       ATH_f = ATHmdct(gfp,freq);  /* freq in kHz */
  263.       ATH_l[sfb]=Min(ATH_l[sfb],ATH_f);
  264.     }
  265.  
  266.  
  267.     /*
  268.     DEBUGF("sfb=%2i freq(khz): %5.2f ..%5.2f  ATH=%6.2f %6.2f  %6.2f   \n",sfb,samp_freq*start/(2*576),
  269. samp_freq*end/(2*576),
  270. 10*log10(ATH_l[sfb]),
  271. 10*log10( ATHmdct(gfp,samp_freq*start/(2*576)))  ,
  272. 10*log10(ATHmdct(gfp,samp_freq*end/(2*576))));
  273.     */
  274.   }
  275.  
  276.   for ( sfb = 0; sfb < SBMAX_s; sfb++ ){
  277.     start = gfc->scalefac_band.s[ sfb ];
  278.     end   = gfc->scalefac_band.s[ sfb+1 ];
  279.     ATH_s[sfb]=1e99;
  280.     for (i=start ; i < end; i++) {
  281.       FLOAT8 freq = samp_freq*i/(2*192);
  282.       assert( freq < 25 );
  283.       ATH_f = ATHmdct(gfp,freq);    /* freq in kHz */
  284.       ATH_s[sfb]=Min(ATH_s[sfb],ATH_f);
  285.     }
  286.   }
  287.   /* in no ATH mode leave ATH for the last scalefactor band in 
  288.    * because VBR mode needs it
  289.    */
  290.   if (gfp->noATH) {
  291.     for ( sfb = 0; sfb < SBMAX_l-1; sfb++ ) {
  292.       ATH_l[sfb]=1E-20;
  293.     }
  294.     for ( sfb = 0; sfb < SBMAX_s-1; sfb++ ) {
  295.       ATH_s[sfb]=1E-20;
  296.     }
  297.   }
  298. }
  299.  
  300.  
  301.  
  302.  
  303.  
  304. /* convert from L/R <-> Mid/Side */
  305. void ms_convert(FLOAT8 xr[2][576],FLOAT8 xr_org[2][576])
  306. {
  307.   int i;
  308.   for ( i = 0; i < 576; i++ ) {
  309.     FLOAT8 l = xr_org[0][i];
  310.     FLOAT8 r = xr_org[1][i];
  311.     xr[0][i] = (l+r)*(SQRT2*0.5);
  312.     xr[1][i] = (l-r)*(SQRT2*0.5);
  313.   }
  314. }
  315.  
  316.  
  317.  
  318. /************************************************************************
  319.  * allocate bits among 2 channels based on PE
  320.  * mt 6/99
  321.  ************************************************************************/
  322. int on_pe(lame_global_flags *gfp,FLOAT8 pe[2][2],III_side_info_t *l3_side,
  323. int targ_bits[2],int mean_bits, int gr)
  324. {
  325.   lame_internal_flags *gfc=gfp->internal_flags;
  326.   gr_info *cod_info;
  327.   int extra_bits,tbits,bits;
  328.   int add_bits[2]; 
  329.   int ch;
  330.   int max_bits;  /* maximum allowed bits for this granule */
  331.  
  332.   /* allocate targ_bits for granule */
  333.   ResvMaxBits(gfp, mean_bits, &tbits, &extra_bits);
  334.   max_bits=tbits+extra_bits;
  335.  
  336.   bits=0;
  337.  
  338.   for (ch=0 ; ch < gfc->stereo ; ch ++) {
  339.     /******************************************************************
  340.      * allocate bits for each channel 
  341.      ******************************************************************/
  342.     cod_info = &l3_side->gr[gr].ch[ch].tt;
  343.     
  344.     targ_bits[ch]=Min(4095,tbits/gfc->stereo);
  345.     
  346.     add_bits[ch]=(pe[gr][ch]-750)/1.4;
  347.     /* short blocks us a little extra, no matter what the pe */
  348.     if (cod_info->block_type==SHORT_TYPE) {
  349.       if (add_bits[ch]<mean_bits/4) add_bits[ch]=mean_bits/4;
  350.     }
  351.  
  352.     /* at most increase bits by 1.5*average */
  353.     if (add_bits[ch] > .75*mean_bits) add_bits[ch]=mean_bits*.75;
  354.     if (add_bits[ch] < 0) add_bits[ch]=0;
  355.  
  356.     if ((targ_bits[ch]+add_bits[ch]) > 4095) 
  357.       add_bits[ch]=Max(0,4095-targ_bits[ch]);
  358.  
  359.     bits += add_bits[ch];
  360.   }
  361.   if (bits > extra_bits)
  362.     for (ch=0 ; ch < gfc->stereo ; ch ++) {
  363.       add_bits[ch] = (extra_bits*add_bits[ch])/bits;
  364.     }
  365.  
  366.   for (ch=0 ; ch < gfc->stereo ; ch ++) {
  367.     targ_bits[ch] = targ_bits[ch] + add_bits[ch];
  368.     extra_bits -= add_bits[ch];
  369.   }
  370.   return max_bits;
  371. }
  372.  
  373.  
  374.  
  375.  
  376. void reduce_side(int targ_bits[2],FLOAT8 ms_ener_ratio,int mean_bits,int max_bits)
  377. {
  378.   int move_bits;
  379.   FLOAT fac;
  380.  
  381.  
  382.   /*  ms_ener_ratio = 0:  allocate 66/33  mid/side  fac=.33  
  383.    *  ms_ener_ratio =.5:  allocate 50/50 mid/side   fac= 0 */
  384.   /* 75/25 split is fac=.5 */
  385.   /* float fac = .50*(.5-ms_ener_ratio[gr])/.5;*/
  386.   fac = .33*(.5-ms_ener_ratio)/.5;
  387.   if (fac<0) fac=0;
  388.   if (fac>.5) fac=.5;
  389.   
  390.     /* number of bits to move from side channel to mid channel */
  391.     /*    move_bits = fac*targ_bits[1];  */
  392.     move_bits = fac*.5*(targ_bits[0]+targ_bits[1]);  
  393.  
  394.     if ((move_bits + targ_bits[0]) > 4095) {
  395.       move_bits = 4095 - targ_bits[0];
  396.     }
  397.     if (move_bits<0) move_bits=0;
  398.     
  399.     if (targ_bits[1] >= 125) {
  400.       /* dont reduce side channel below 125 bits */
  401.       if (targ_bits[1]-move_bits > 125) {
  402.  
  403.     /* if mid channel already has 2x more than average, dont bother */
  404.     /* mean_bits = bits per granule (for both channels) */
  405.     if (targ_bits[0] < mean_bits)
  406.       targ_bits[0] += move_bits;
  407.     targ_bits[1] -= move_bits;
  408.       } else {
  409.     targ_bits[0] += targ_bits[1] - 125;
  410.     targ_bits[1] = 125;
  411.       }
  412.     }
  413.     
  414.     move_bits=targ_bits[0]+targ_bits[1];
  415.     if (move_bits > max_bits) {
  416.       targ_bits[0]=(max_bits*targ_bits[0])/move_bits;
  417.       targ_bits[1]=(max_bits*targ_bits[1])/move_bits;
  418.     }
  419. }
  420.  
  421. /*************************************************************************** 
  422.  *         inner_loop                                                      * 
  423.  *************************************************************************** 
  424.  * The code selects the best global gain for a particular set of scalefacs */
  425.  
  426. int
  427. inner_loop( lame_global_flags *gfp,FLOAT8 xrpow[576],
  428.         int l3_enc[576], int max_bits,
  429.         gr_info *cod_info)
  430. {
  431.     int bits;
  432.     assert( max_bits >= 0 );
  433.     cod_info->global_gain--;
  434.  
  435.     do
  436.     {
  437.       cod_info->global_gain++;
  438.       bits = count_bits(gfp,l3_enc, xrpow, cod_info);
  439.     }
  440.     while ( bits > max_bits );
  441.  
  442.     return bits;
  443. }
  444.  
  445.  
  446.  
  447. /*************************************************************************/
  448. /*            calc_xmin                                                  */
  449. /*************************************************************************/
  450.  
  451. /*
  452.   Calculate the allowed distortion for each scalefactor band,
  453.   as determined by the psychoacoustic model.
  454.   xmin(sb) = ratio(sb) * en(sb) / bw(sb)
  455.  
  456.   returns number of sfb's with energy > ATH
  457. */
  458. int calc_xmin( lame_global_flags *gfp,FLOAT8 xr[576], III_psy_ratio *ratio,
  459.            gr_info *cod_info, III_psy_xmin *l3_xmin)
  460. {
  461.   lame_internal_flags *gfc=gfp->internal_flags;
  462.   int j,start, end, bw,l, b, ath_over=0;
  463.   u_int    sfb;
  464.   FLOAT8 en0, xmin, ener;
  465.  
  466.   if (cod_info->block_type==SHORT_TYPE) {
  467.  
  468.   for ( j=0, sfb = 0; sfb < SBMAX_s; sfb++ ) {
  469.     start = gfc->scalefac_band.s[ sfb ];
  470.     end   = gfc->scalefac_band.s[ sfb + 1 ];
  471.     bw = end - start;
  472.     for ( b = 0; b < 3; b++ ) {
  473.       for (en0 = 0.0, l = start; l < end; l++) {
  474.         ener = xr[j++];
  475.         ener = ener * ener;
  476.         en0 += ener;
  477.       }
  478.       en0 /= bw;
  479.       
  480.       if (gfp->ATHonly || gfp->ATHshort) {
  481.         l3_xmin->s[sfb][b]=gfc->ATH_s[sfb];
  482.       } else {
  483.         xmin = ratio->en.s[sfb][b];
  484.         if (xmin > 0.0)
  485.           xmin = en0 * ratio->thm.s[sfb][b] * gfc->masking_lower / xmin;
  486.         l3_xmin->s[sfb][b] = Max(gfc->ATH_s[sfb], xmin);
  487.       }
  488.       
  489.       if (en0 > gfc->ATH_s[sfb]) ath_over++;
  490.     }
  491.   }
  492.  
  493.   }else{
  494.   
  495.   for ( sfb = 0; sfb < SBMAX_l; sfb++ ){
  496.     start = gfc->scalefac_band.l[ sfb ];
  497.     end   = gfc->scalefac_band.l[ sfb+1 ];
  498.     bw = end - start;
  499.     
  500.     for (en0 = 0.0, l = start; l < end; l++ ) {
  501.       ener = xr[l] * xr[l];
  502.       en0 += ener;
  503.     }
  504.     en0 /= bw;
  505.     
  506.     if (gfp->ATHonly) {
  507.       l3_xmin->l[sfb]=gfc->ATH_l[sfb];
  508.     } else {
  509.       xmin = ratio->en.l[sfb];
  510.       if (xmin > 0.0)
  511.         xmin = en0 * ratio->thm.l[sfb] * gfc->masking_lower / xmin;
  512.       l3_xmin->l[sfb]=Max(gfc->ATH_l[sfb], xmin);
  513.     }
  514.     if (en0 > gfc->ATH_l[sfb]) ath_over++;
  515.   }
  516.   }
  517.   return ath_over;
  518. }
  519.  
  520.  
  521.  
  522.  
  523.  
  524.  
  525.  
  526. /*************************************************************************/
  527. /*            calc_noise                                                 */
  528. /*************************************************************************/
  529. /*  mt 5/99:  Function: Improved calc_noise for a single channel   */
  530. int calc_noise( lame_global_flags *gfp,
  531.                  FLOAT8 xr[576], int ix[576], gr_info *cod_info,
  532.          FLOAT8 xfsf[4][SBMAX_l], FLOAT8 distort[4][SBMAX_l],
  533.          III_psy_xmin *l3_xmin, III_scalefac_t *scalefac,
  534.          calc_noise_result *res)
  535. {
  536.   int start, end, j,l, i, over=0;
  537.   u_int sfb;
  538.   FLOAT8 sum, osum, bw;
  539.   lame_internal_flags *gfc=gfp->internal_flags;
  540.   
  541.   int count=0;
  542.   FLOAT8 noise;
  543.   FLOAT8 over_noise = 1;     /*    0 dB relative to masking */
  544.   FLOAT8 tot_noise  = 1;     /*    0 dB relative to masking */
  545.   FLOAT8 max_noise  = 1E-20; /* -200 dB relative to masking */
  546.  
  547.   
  548.   if (cod_info->block_type == SHORT_TYPE) {
  549.     int max_index = SBPSY_s;
  550.     if (gfp->VBR==vbr_rh || gfp->VBR==vbr_mt)
  551.       {
  552.         max_index = SBMAX_s;
  553.       }
  554.     for ( j=0, sfb = 0; sfb < max_index; sfb++ ) {
  555.          start = gfc->scalefac_band.s[ sfb ];
  556.          end   = gfc->scalefac_band.s[ sfb+1 ];
  557.          bw = end - start;
  558.          for ( i = 0; i < 3; i++ ) {
  559.         FLOAT8 step;
  560.         int s;
  561.  
  562.         s = (scalefac->s[sfb][i] << (cod_info->scalefac_scale + 1))
  563.         + cod_info->subblock_gain[i] * 8;
  564.         s = cod_info->global_gain - s;
  565.  
  566.         assert(s<Q_MAX);
  567.         assert(s>=0);
  568.         step = POW20(s);
  569.  
  570.         if (gfp->exp_nspsytune) {
  571.           for ( osum = sum = 0.0, l = start; l < end; l++ ) {
  572.         FLOAT8 temp;
  573.         temp = fabs(xr[j]) - pow43[ix[j]] * step;
  574.         ++j;
  575.  
  576.         temp = temp*temp;
  577.         osum = Max(osum,bw*temp);
  578.         sum += temp;
  579.           }       
  580.  
  581.           if (osum > sum*NSTHRE) sum = osum;
  582.         } else {
  583.           for ( sum = 0.0, l = start; l < end; l++ ) {
  584.         FLOAT8 temp;
  585.         temp = fabs(xr[j]) - pow43[ix[j]] * step;
  586.         ++j;
  587.         sum += temp * temp;
  588.           }   
  589.         }    
  590.  
  591.         xfsf[i+1][sfb] = sum / bw;
  592.         noise = xfsf[i+1][sfb] / l3_xmin->s[sfb][i];
  593.         /* multiplying here is adding in dB */
  594.         tot_noise *= Max(noise, 1E-20);
  595.  
  596.             if (noise > 1) {
  597.         over++;
  598.             /* multiplying here is adding in dB */
  599.         over_noise *= noise;
  600.         }
  601.         max_noise=Max(max_noise,noise);
  602.             distort[i+1][sfb] = noise;
  603.  
  604.         count++;        
  605.         }
  606.     }
  607.   }else{
  608.     int max_index = SBPSY_l;
  609.     
  610.     if (gfp->VBR==vbr_rh || gfp->VBR==vbr_mt)
  611.       {
  612.         max_index = SBMAX_l;
  613.       }
  614.     for ( sfb = 0; sfb < max_index; sfb++ ) {
  615.     FLOAT8 step;
  616.     int s = scalefac->l[sfb];
  617.  
  618.     if (cod_info->preflag)
  619.         s += pretab[sfb];
  620.  
  621.     s = cod_info->global_gain - (s << (cod_info->scalefac_scale + 1));
  622.     assert(s<Q_MAX);
  623.     assert(s>=0);
  624.     step = POW20(s);
  625.  
  626.     start = gfc->scalefac_band.l[ sfb ];
  627.         end   = gfc->scalefac_band.l[ sfb+1 ];
  628.         bw = end - start;
  629.  
  630.     if (gfp->exp_nspsytune) {
  631.       for ( osum = sum = 0.0, l = start; l < end; l++ )
  632.         {
  633.           FLOAT8 temp;
  634.           temp = fabs(xr[l]) - pow43[ix[l]] * step;
  635.           temp = temp*temp;
  636.           osum = Max(osum,bw*temp);
  637.           sum += temp;
  638.         }
  639.  
  640.       if (osum > sum*NSTHRE) sum = osum;
  641.     } else {
  642.       for ( sum = 0.0, l = start; l < end; l++ )
  643.         {
  644.           FLOAT8 temp;
  645.           temp = fabs(xr[l]) - pow43[ix[l]] * step;
  646.           sum += temp * temp;
  647.         }
  648.     }
  649.  
  650.         xfsf[0][sfb] = sum / bw;
  651.  
  652.     noise = xfsf[0][sfb] / l3_xmin->l[sfb];
  653.     /* multiplying here is adding in dB */
  654.         tot_noise *= Max(noise, 1E-20);
  655.         if (noise>1) {
  656.       over++;
  657.       /* multiplying here is adding in dB */
  658.       over_noise *= noise;
  659.     }
  660.     max_noise=Max(max_noise,noise);
  661.         distort[0][sfb] = noise;
  662.     count++;
  663.  
  664.     }
  665.  
  666.   }
  667.  
  668.   /* normalization at this point by "count" is not necessary, since
  669.    * the values are only used to compare with previous values */
  670.   res->tot_count  = count;
  671.   res->over_count = over;
  672.  
  673. #ifndef RH_NOISE_CALC
  674.   /* convert to db. DO NOT CHANGE THESE */
  675.   /* tot_noise = is really the average over each sfb of: 
  676.      [noise(db) - allowed_noise(db)]
  677.  
  678.      and over_noise is the same average, only over only the
  679.      bands with noise > allowed_noise.  
  680.  
  681.   */
  682.  
  683.   res->tot_noise = 10*log10(Max(.00001,tot_noise)); 
  684.   res->over_noise = 10*log10(Max(1.0,over_noise)); 
  685.   res->max_noise = 10*log10(Max(.00001,max_noise));
  686. #else
  687.   /* no need for dB. DO NOT CHANGE THESE */
  688.   res->tot_noise  = tot_noise; 
  689.   res->over_noise = over_noise; 
  690.   res->max_noise  = max_noise;
  691. #endif  
  692.   return over;
  693. }
  694.  
  695.  
  696.  
  697.  
  698.  
  699.  
  700.  
  701.  
  702.  
  703.  
  704.  
  705.  
  706.  
  707.  
  708. /*************************************************************************/
  709. /*            loop_break                                                 */
  710. /*************************************************************************/
  711.  
  712. /*  Function: Returns zero if there is a scalefac which has not been
  713.     amplified. Otherwise it returns one. 
  714. */
  715.  
  716. int loop_break( III_scalefac_t *scalefac, gr_info *cod_info)
  717. {
  718.     int i;
  719.     u_int sfb;
  720.  
  721.     for ( sfb = 0; sfb < cod_info->sfb_lmax; sfb++ )
  722.         if ( scalefac->l[sfb] == 0 )
  723.         return 0;
  724.  
  725.     for ( sfb = cod_info->sfb_smax; sfb < SBPSY_s; sfb++ )
  726.       for ( i = 0; i < 3; i++ ) 
  727.             if ( scalefac->s[sfb][i]+cod_info->subblock_gain[i] == 0 )
  728.         return 0;
  729.  
  730.     return 1;
  731. }
  732.  
  733.  
  734.  
  735.  
  736.  
  737.  
  738.  
  739.  
  740.  
  741.  
  742.  
  743.  
  744.  
  745. /*
  746.  ----------------------------------------------------------------------
  747.   if someone wants to try to find a faster step search function,
  748.   here is some code which gives a lower bound for the step size:
  749.   
  750.   for (max_xrspow = 0, i = 0; i < 576; ++i)
  751.   {
  752.     max_xrspow = Max(max_xrspow, xrspow[i]);
  753.   }
  754.   lowerbound = 210+log10(max_xrspow/IXMAX_VAL)/(0.1875*LOG2);
  755.  
  756.  
  757.                                                  Robert.Hegemann@gmx.de
  758.  ----------------------------------------------------------------------
  759. */
  760.  
  761.  
  762. typedef enum {
  763.     BINSEARCH_NONE,
  764.     BINSEARCH_UP, 
  765.     BINSEARCH_DOWN
  766. } binsearchDirection_t;
  767.  
  768. /*-------------------------------------------------------------------------*/
  769. int 
  770. bin_search_StepSize2 (lame_global_flags *gfp,int desired_rate, int start, int *ix, 
  771.                       FLOAT8 xrspow[576], gr_info *cod_info)
  772. /*-------------------------------------------------------------------------*/
  773. {
  774.     int nBits;
  775.     int flag_GoneOver = 0;
  776.     int StepSize = start;
  777.     int CurrentStep;
  778.     lame_internal_flags *gfc=gfp->internal_flags;
  779.  
  780.     binsearchDirection_t Direction = BINSEARCH_NONE;
  781.     assert(gfc->CurrentStep);
  782.     CurrentStep = gfc->CurrentStep;
  783.  
  784.     do
  785.     {
  786.     cod_info->global_gain = StepSize;
  787.     nBits = count_bits(gfp,ix, xrspow, cod_info);  
  788.  
  789.     if (CurrentStep == 1 )
  790.         {
  791.         break; /* nothing to adjust anymore */
  792.     }
  793.     if (flag_GoneOver)
  794.     {
  795.         CurrentStep /= 2;
  796.     }
  797.     if (nBits > desired_rate)  /* increase Quantize_StepSize */
  798.     {
  799.         if (Direction == BINSEARCH_DOWN && !flag_GoneOver)
  800.         {
  801.         flag_GoneOver = 1;
  802.         CurrentStep /= 2; /* late adjust */
  803.         }
  804.         Direction = BINSEARCH_UP;
  805.         StepSize += CurrentStep;
  806.         if (StepSize > 255) break;
  807.     }
  808.     else if (nBits < desired_rate)
  809.     {
  810.         if (Direction == BINSEARCH_UP && !flag_GoneOver)
  811.         {
  812.         flag_GoneOver = 1;
  813.         CurrentStep /= 2; /* late adjust */
  814.         }
  815.         Direction = BINSEARCH_DOWN;
  816.         StepSize -= CurrentStep;
  817.         if (StepSize < 0) break;
  818.     }
  819.     else break; /* nBits == desired_rate;; most unlikely to happen.*/
  820.     } while (1); /* For-ever, break is adjusted. */
  821.  
  822.     CurrentStep = abs(start - StepSize);
  823.     
  824.     if (CurrentStep >= 4) {
  825.     CurrentStep = 4;
  826.     } else {
  827.     CurrentStep = 2;
  828.     }
  829.  
  830.     gfc->CurrentStep = CurrentStep;
  831.  
  832.     return nBits;
  833. }
  834.  
  835.  
  836.  
  837.  
  838. /************************************************************************/
  839. /*  updates plotting data                                               */
  840. /************************************************************************/
  841. void 
  842. set_pinfo (lame_global_flags *gfp,
  843.     gr_info *cod_info,
  844.     III_psy_ratio *ratio, 
  845.     III_scalefac_t *scalefac,
  846.     FLOAT8 xr[576],        
  847.     FLOAT8 xfsf[4][SBMAX_l],
  848.     FLOAT8 noise[4],
  849.     int gr,
  850.     int ch
  851. )
  852. {
  853.   lame_internal_flags *gfc=gfp->internal_flags;
  854.   int sfb;
  855.   int j,i,l,start,end,bw;
  856.   FLOAT8 en0,en1;
  857.   FLOAT ifqstep = ( cod_info->scalefac_scale == 0 ) ? .5 : 1.0;
  858.  
  859.  
  860.   if (cod_info->block_type == SHORT_TYPE) {
  861.     for (j=0, sfb = 0; sfb < SBMAX_s; sfb++ )  {
  862.       start = gfc->scalefac_band.s[ sfb ];
  863.       end   = gfc->scalefac_band.s[ sfb + 1 ];
  864.       bw = end - start;
  865.       for ( i = 0; i < 3; i++ ) {
  866.     for ( en0 = 0.0, l = start; l < end; l++ ) {
  867.       en0 += xr[j] * xr[j];
  868.       ++j;
  869.     }
  870.     en0=Max(en0/bw,1e-20);
  871.  
  872.  
  873. #if 0
  874.     {
  875. double tot1,tot2;
  876.       if (sfb<SBMAX_s-1) {
  877.       if (sfb==0) {
  878.     tot1=0;
  879.     tot2=0;
  880.       }
  881.       tot1 += en0;
  882.       tot2 += ratio->en.s[sfb][i];
  883.  
  884.  
  885.       DEBUGF("%i %i sfb=%i mdct=%f fft=%f  fft-mdct=%f db \n",gr,ch,sfb,
  886. 10*log10(Max(1e-25,ratio->en.s[sfb][i])),
  887. 10*log10(Max(1e-25,en0)),
  888. 10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.s[sfb][i]))));
  889.  
  890.       if (sfb==SBMAX_s-2) {
  891.       DEBUGF("%i %i toti %f %f ratio=%f db \n",gr,ch,
  892. 10*log10(Max(1e-25,tot2)),
  893. 10*log(Max(1e-25,tot1)),
  894. 10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
  895.  
  896.       }
  897.       }
  898.       /*
  899. masking: multiplied by en0/fft_energy
  900. average seems to be about -135db.
  901.        */
  902.     }
  903. #endif
  904.  
  905.  
  906.  
  907.  
  908.     /* convert to MDCT units */
  909.     en1=1e15;  /* scaling so it shows up on FFT plot */
  910.     gfc->pinfo->xfsf_s[gr][ch][3*sfb+i] =  en1*xfsf[i+1][sfb];
  911.     gfc->pinfo->en_s[gr][ch][3*sfb+i] = en1*en0;
  912.     
  913.     if (ratio->en.s[sfb][i]>0)
  914.       en0 = en0/ratio->en.s[sfb][i];
  915.     else
  916.       en0=0;
  917.     if (gfp->ATHonly || gfp->ATHshort)
  918.       en0=0;
  919.  
  920.     gfc->pinfo->thr_s[gr][ch][3*sfb+i] = 
  921. en1*Max(en0*ratio->thm.s[sfb][i],gfc->ATH_s[sfb]);
  922.  
  923.     
  924.     /* there is no scalefactor bands >= SBPSY_s */
  925.     if (sfb < SBPSY_s) {
  926.       gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=
  927.         -ifqstep*scalefac->s[sfb][i];
  928.     } else {
  929.       gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=0;
  930.     }
  931.     gfc->pinfo->LAMEsfb_s[gr][ch][3*sfb+i] -= 2*cod_info->subblock_gain[i];
  932.       }
  933.     }
  934.   }else{
  935.     for ( sfb = 0; sfb < SBMAX_l; sfb++ )   {
  936.  
  937.       start = gfc->scalefac_band.l[ sfb ];
  938.       end   = gfc->scalefac_band.l[ sfb+1 ];
  939.       bw = end - start;
  940.       for ( en0 = 0.0, l = start; l < end; l++ ) 
  941.     en0 += xr[l] * xr[l];
  942.       en0/=bw;
  943.       /*
  944.     DEBUGF("diff  = %f \n",10*log10(Max(ratio[gr][ch].en.l[sfb],1e-20))
  945.     -(10*log10(en0)+150));
  946.       */
  947.  
  948. #if 0
  949.     {
  950. double tot1,tot2;
  951.       if (sfb==0) {
  952.     tot1=0;
  953.     tot2=0;
  954.       }
  955.       tot1 += en0;
  956.       tot2 += ratio->en.l[sfb];
  957.  
  958.  
  959.       DEBUGF("%i %i sfb=%i mdct=%f fft=%f  fft-mdct=%f db \n",gr,ch,sfb,
  960. 10*log10(Max(1e-25,ratio->en.l[sfb])),
  961. 10*log10(Max(1e-25,en0)),
  962. 10*log10((Max(1e-25,en0)/Max(1e-25,ratio->en.l[sfb]))));
  963.  
  964.       if (sfb==SBMAX_l-1) {
  965.       DEBUGF("%i %i toti %f %f ratio=%f db \n",gr,ch,
  966. 10*log10(Max(1e-25,tot2)),
  967. 10*log(Max(1e-25,tot1)),
  968. 10*log10(Max(1e-25,tot1)/(Max(1e-25,tot2))));
  969.  
  970.       }
  971.       /*
  972. masking: multiplied by en0/fft_energy
  973. average seems to be about -147db.
  974.  
  975.        */
  976.     }
  977. #endif
  978.  
  979.  
  980.  
  981.       /* convert to MDCT units */
  982.       en1=1e15;  /* scaling so it shows up on FFT plot */
  983.       gfc->pinfo->xfsf[gr][ch][sfb] =  en1*xfsf[0][sfb];
  984.       gfc->pinfo->en[gr][ch][sfb] = en1*en0;
  985.       if (ratio->en.l[sfb]>0)
  986.     en0 = en0/ratio->en.l[sfb];
  987.       else
  988.     en0=0;
  989.       if (gfp->ATHonly)
  990.     en0=0;
  991.       gfc->pinfo->thr[gr][ch][sfb] = en1*Max(en0*ratio->thm.l[sfb],gfc->ATH_l[sfb]);
  992.  
  993.       /* there is no scalefactor bands >= SBPSY_l */
  994.       if (sfb<SBPSY_l) {
  995.     if (scalefac->l[sfb]<0)  /* scfsi! */
  996.       gfc->pinfo->LAMEsfb[gr][ch][sfb]=gfc->pinfo->LAMEsfb[0][ch][sfb];
  997.     else
  998.       gfc->pinfo->LAMEsfb[gr][ch][sfb]=-ifqstep*scalefac->l[sfb];
  999.       }else{
  1000.     gfc->pinfo->LAMEsfb[gr][ch][sfb]=0;
  1001.       }
  1002.  
  1003.       if (cod_info->preflag && sfb>=11) 
  1004.     gfc->pinfo->LAMEsfb[gr][ch][sfb]-=ifqstep*pretab[sfb];
  1005.     }
  1006.   }
  1007.   gfc->pinfo->LAMEqss     [gr][ch] = cod_info->global_gain;
  1008.   gfc->pinfo->LAMEmainbits[gr][ch] = cod_info->part2_3_length;
  1009.   gfc->pinfo->LAMEsfbits  [gr][ch] = cod_info->part2_length;
  1010.  
  1011.   gfc->pinfo->over      [gr][ch] = noise[0];
  1012. #ifndef RH_NOISE_CALC
  1013.   gfc->pinfo->max_noise [gr][ch] = noise[1];
  1014.   gfc->pinfo->over_noise[gr][ch] = noise[2];
  1015.   gfc->pinfo->tot_noise [gr][ch] = noise[3];
  1016. #else
  1017.   gfc->pinfo->max_noise [gr][ch] = 10*log10(Max(1E-20,noise[1]));
  1018.   gfc->pinfo->over_noise[gr][ch] = 10*log10(Max(1.0,  noise[2]));
  1019.   gfc->pinfo->tot_noise [gr][ch] = 10*log10(Max(1E-20,noise[3]));
  1020. #endif
  1021. }
  1022.  
  1023.  
  1024.  
  1025.  
  1026.  
  1027.  
  1028. /*********************************************************************
  1029.  * nonlinear quantization of xr 
  1030.  * More accurate formula than the ISO formula.  Takes into account
  1031.  * the fact that we are quantizing xr -> ix, but we want ix^4/3 to be 
  1032.  * as close as possible to x^4/3.  (taking the nearest int would mean
  1033.  * ix is as close as possible to xr, which is different.)
  1034.  * From Segher Boessenkool <segher@eastsite.nl>  11/1999
  1035.  * ASM optimization from 
  1036.  *    Mathew Hendry <scampi@dial.pipex.com> 11/1999
  1037.  *    Acy Stapp <AStapp@austin.rr.com> 11/1999
  1038.  *    Takehiro Tominaga <tominaga@isoternet.org> 11/1999
  1039.  *********************************************************************/
  1040.  
  1041. #ifdef TAKEHIRO_IEEE754_HACK
  1042.  
  1043. #define MAGIC_FLOAT (65536*(128))
  1044. #define MAGIC_INT 0x4b000000
  1045.  
  1046. void quantize_xrpow(FLOAT8 xp[576], int pi[576], gr_info *cod_info)
  1047. {
  1048.     /* quantize on xr^(3/4) instead of xr */
  1049.     const FLOAT8 istep = IPOW20(cod_info->global_gain);
  1050.  
  1051.     int j;
  1052.     for (j = 576 / 4; j > 0; --j) {
  1053.     double x0 = istep * xp[0] + MAGIC_FLOAT;
  1054.     double x1 = istep * xp[1] + MAGIC_FLOAT;
  1055.     double x2 = istep * xp[2] + MAGIC_FLOAT;
  1056.     double x3 = istep * xp[3] + MAGIC_FLOAT;
  1057.  
  1058.     ((float*)pi)[0] = x0;
  1059.     ((float*)pi)[1] = x1;
  1060.     ((float*)pi)[2] = x2;
  1061.     ((float*)pi)[3] = x3;
  1062.  
  1063.     ((float *)pi)[0] = (x0 + adj43asm[pi[0] - MAGIC_INT]);
  1064.     ((float *)pi)[1] = (x1 + adj43asm[pi[1] - MAGIC_INT]);
  1065.     ((float *)pi)[2] = (x2 + adj43asm[pi[2] - MAGIC_INT]);
  1066.     ((float *)pi)[3] = (x3 + adj43asm[pi[3] - MAGIC_INT]);
  1067.  
  1068.     pi[0] -= MAGIC_INT;
  1069.     pi[1] -= MAGIC_INT;
  1070.     pi[2] -= MAGIC_INT;
  1071.     pi[3] -= MAGIC_INT;
  1072.     pi += 4;
  1073.     xp += 4;
  1074.     }
  1075. }
  1076.  
  1077. #  define ROUNDFAC -0.0946
  1078. void quantize_xrpow_ISO(FLOAT8 xp[576], int pi[576], gr_info *cod_info)
  1079. {
  1080.     /* quantize on xr^(3/4) instead of xr */
  1081.     const FLOAT8 istep = IPOW20(cod_info->global_gain);
  1082.  
  1083.     register int j;
  1084.     for (j=576/4;j>0;j--) {
  1085.     ((float *)pi)[0] = (istep * xp[0]) + (ROUNDFAC + MAGIC_FLOAT);
  1086.     ((float *)pi)[1] = (istep * xp[1]) + (ROUNDFAC + MAGIC_FLOAT);
  1087.     ((float *)pi)[2] = (istep * xp[2]) + (ROUNDFAC + MAGIC_FLOAT);
  1088.     ((float *)pi)[3] = (istep * xp[3]) + (ROUNDFAC + MAGIC_FLOAT);
  1089.  
  1090.     pi[0] -= MAGIC_INT;
  1091.     pi[1] -= MAGIC_INT;
  1092.     pi[2] -= MAGIC_INT;
  1093.     pi[3] -= MAGIC_INT;
  1094.     pi+=4;
  1095.     xp+=4;
  1096.     }
  1097. }
  1098.  
  1099. #else
  1100.  
  1101.  
  1102.  
  1103.  
  1104.  
  1105.  
  1106. #if (defined(__GNUC__) && defined(__i386__))
  1107. #define USE_GNUC_ASM
  1108. #endif
  1109. #ifdef _MSC_VER
  1110. #define USE_MSC_ASM
  1111. #endif
  1112.  
  1113.  
  1114.  
  1115. /*********************************************************************
  1116.  * XRPOW_FTOI is a macro to convert floats to ints.  
  1117.  * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x]
  1118.  *                                         ROUNDFAC= -0.0946
  1119.  *
  1120.  * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]   
  1121.  *                                   ROUNDFAC=0.4054
  1122.  *********************************************************************/
  1123. #ifdef USE_GNUC_ASM
  1124. #  define QUANTFAC(rx)  adj43asm[rx]
  1125. #  define ROUNDFAC -0.0946
  1126. #  define XRPOW_FTOI(src, dest) \
  1127.      asm ("fistpl %0 " : "=m"(dest) : "t"(src) : "st")
  1128. #elif defined (USE_MSC_ASM)
  1129. #  define QUANTFAC(rx)  adj43asm[rx]
  1130. #  define ROUNDFAC -0.0946
  1131. #  define XRPOW_FTOI(src, dest) do { \
  1132.      FLOAT8 src_ = (src); \
  1133.      int dest_; \
  1134.      { \
  1135.        __asm fld src_ \
  1136.        __asm fistp dest_ \
  1137.      } \
  1138.      (dest) = dest_; \
  1139.    } while (0)
  1140. #else
  1141. #  define QUANTFAC(rx)  adj43[rx]
  1142. #  define ROUNDFAC 0.4054
  1143. #  define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
  1144. #endif
  1145.  
  1146. #ifdef USE_MSC_ASM
  1147. /* define F8type and F8size according to type of FLOAT8 */
  1148. # if defined FLOAT8_is_double
  1149. #  define F8type qword
  1150. #  define F8size 8
  1151. # elif defined FLOAT8_is_float
  1152. #  define F8type dword
  1153. #  define F8size 4
  1154. # else
  1155. /* only float and double supported */
  1156. #  error invalid FLOAT8 type for USE_MSC_ASM
  1157. # endif
  1158. #endif
  1159.  
  1160. #ifdef USE_GNUC_ASM
  1161. /* define F8type and F8size according to type of FLOAT8 */
  1162. # if defined FLOAT8_is_double
  1163. #  define F8type "l"
  1164. #  define F8size "8"
  1165. # elif defined FLOAT8_is_float
  1166. #  define F8type "s"
  1167. #  define F8size "4"
  1168. # else
  1169. /* only float and double supported */
  1170. #  error invalid FLOAT8 type for USE_GNUC_ASM
  1171. # endif
  1172. #endif
  1173.  
  1174. void quantize_xrpow(FLOAT8 xr[576], int ix[576], gr_info *cod_info) {
  1175.   /* quantize on xr^(3/4) instead of xr */
  1176.   const FLOAT8 istep = IPOW20(cod_info->global_gain);
  1177.  
  1178. #if defined (USE_GNUC_ASM) 
  1179.   {
  1180.       int rx[4];
  1181.       __asm__ __volatile__(
  1182.         "\n\nloop1:\n\t"
  1183.  
  1184.         "fld" F8type " 0*" F8size "(%1)\n\t"
  1185.         "fld" F8type " 1*" F8size "(%1)\n\t"
  1186.         "fld" F8type " 2*" F8size "(%1)\n\t"
  1187.         "fld" F8type " 3*" F8size "(%1)\n\t"
  1188.  
  1189.         "fxch %%st(3)\n\t"
  1190.         "fmul %%st(4)\n\t"
  1191.         "fxch %%st(2)\n\t"
  1192.         "fmul %%st(4)\n\t"
  1193.         "fxch %%st(1)\n\t"
  1194.         "fmul %%st(4)\n\t"
  1195.         "fxch %%st(3)\n\t"
  1196.         "fmul %%st(4)\n\t"
  1197.  
  1198.         "addl $4*" F8size ", %1\n\t"
  1199.         "addl $16, %3\n\t"
  1200.  
  1201.         "fxch %%st(2)\n\t"
  1202.         "fistl %5\n\t"
  1203.         "fxch %%st(1)\n\t"
  1204.         "fistl 4+%5\n\t"
  1205.         "fxch %%st(3)\n\t"
  1206.         "fistl 8+%5\n\t"
  1207.         "fxch %%st(2)\n\t"
  1208.         "fistl 12+%5\n\t"
  1209.  
  1210.         "dec %4\n\t"
  1211.  
  1212.         "movl %5, %%eax\n\t"
  1213.         "movl 4+%5, %%ebx\n\t"
  1214.         "fxch %%st(1)\n\t"
  1215.         "fadd" F8type " (%2,%%eax," F8size ")\n\t"
  1216.         "fxch %%st(3)\n\t"
  1217.         "fadd" F8type " (%2,%%ebx," F8size ")\n\t"
  1218.  
  1219.         "movl 8+%5, %%eax\n\t"
  1220.         "movl 12+%5, %%ebx\n\t"
  1221.         "fxch %%st(2)\n\t"
  1222.         "fadd" F8type " (%2,%%eax," F8size ")\n\t"
  1223.         "fxch %%st(1)\n\t"
  1224.         "fadd" F8type " (%2,%%ebx," F8size ")\n\t"
  1225.  
  1226.         "fxch %%st(3)\n\t"
  1227.         "fistpl -16(%3)\n\t"
  1228.         "fxch %%st(1)\n\t"
  1229.         "fistpl -12(%3)\n\t"
  1230.         "fistpl -8(%3)\n\t"
  1231.         "fistpl -4(%3)\n\t"
  1232.  
  1233.         "jnz loop1\n\n"
  1234.         : /* no outputs */
  1235.         : "t" (istep), "r" (xr), "r" (adj43asm), "r" (ix), "r" (576 / 4), "m" (rx)
  1236.         : "%eax", "%ebx", "memory", "cc"
  1237.       );
  1238.   }
  1239. #elif defined (USE_MSC_ASM)
  1240.   {
  1241.       /* asm from Acy Stapp <AStapp@austin.rr.com> */
  1242.       int rx[4];
  1243.       _asm {
  1244.           fld F8type ptr [istep]
  1245.           mov esi, dword ptr [xr]
  1246.           lea edi, dword ptr [adj43asm]
  1247.           mov edx, dword ptr [ix]
  1248.           mov ecx, 576/4
  1249.       } loop1: _asm {
  1250.           fld F8type ptr [esi+(0*F8size)] // 0
  1251.           fld F8type ptr [esi+(1*F8size)] // 1 0
  1252.           fld F8type ptr [esi+(2*F8size)] // 2 1 0
  1253.           fld F8type ptr [esi+(3*F8size)] // 3 2 1 0
  1254.           fxch st(3)                  // 0 2 1 3
  1255.           fmul st(0), st(4)
  1256.           fxch st(2)                  // 1 2 0 3
  1257.           fmul st(0), st(4)
  1258.           fxch st(1)                  // 2 1 0 3
  1259.           fmul st(0), st(4)
  1260.           fxch st(3)                  // 3 1 0 2
  1261.           fmul st(0), st(4)
  1262.  
  1263.           add esi, 4*F8size
  1264.           add edx, 16
  1265.  
  1266.           fxch st(2)                  // 0 1 3 2
  1267.           fist dword ptr [rx]
  1268.           fxch st(1)                  // 1 0 3 2
  1269.           fist dword ptr [rx+4]
  1270.           fxch st(3)                  // 2 0 3 1
  1271.           fist dword ptr [rx+8]
  1272.           fxch st(2)                  // 3 0 2 1
  1273.           fist dword ptr [rx+12]
  1274.  
  1275.           dec ecx
  1276.  
  1277.           mov eax, dword ptr [rx]
  1278.           mov ebx, dword ptr [rx+4]
  1279.           fxch st(1)                  // 0 3 2 1
  1280.           fadd F8type ptr [edi+eax*F8size]
  1281.           fxch st(3)                  // 1 3 2 0
  1282.           fadd F8type ptr [edi+ebx*F8size]
  1283.  
  1284.           mov eax, dword ptr [rx+8]
  1285.           mov ebx, dword ptr [rx+12]
  1286.           fxch st(2)                  // 2 3 1 0
  1287.           fadd F8type ptr [edi+eax*F8size]
  1288.           fxch st(1)                  // 3 2 1 0
  1289.           fadd F8type ptr [edi+ebx*F8size]
  1290.           fxch st(3)                  // 0 2 1 3
  1291.           fistp dword ptr [edx-16]    // 2 1 3
  1292.           fxch st(1)                  // 1 2 3
  1293.           fistp dword ptr [edx-12]    // 2 3
  1294.           fistp dword ptr [edx-8]     // 3
  1295.           fistp dword ptr [edx-4]
  1296.  
  1297.           jnz loop1
  1298.  
  1299.           mov dword ptr [xr], esi
  1300.           mov dword ptr [ix], edx
  1301.           fstp st(0)
  1302.       }
  1303.   }
  1304. #else
  1305. #if 0
  1306.   {   /* generic code if you write ASM for XRPOW_FTOI() */
  1307.       FLOAT8 x;
  1308.       int j, rx;
  1309.       for (j = 576 / 4; j > 0; --j) {
  1310.           x = *xr++ * istep;
  1311.           XRPOW_FTOI(x, rx);
  1312.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1313.  
  1314.           x = *xr++ * istep;
  1315.           XRPOW_FTOI(x, rx);
  1316.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1317.  
  1318.           x = *xr++ * istep;
  1319.           XRPOW_FTOI(x, rx);
  1320.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1321.  
  1322.           x = *xr++ * istep;
  1323.           XRPOW_FTOI(x, rx);
  1324.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1325.       }
  1326.   }
  1327. #endif
  1328.   {/* from Wilfried.Behne@t-online.de.  Reported to be 2x faster than 
  1329.       the above code (when not using ASM) on PowerPC */
  1330.          int j;
  1331.          
  1332.          for ( j = 576/8; j > 0; --j)
  1333.          {
  1334.             FLOAT8    x1, x2, x3, x4, x5, x6, x7, x8;
  1335.             int        rx1, rx2, rx3, rx4, rx5, rx6, rx7, rx8;
  1336.             x1 = *xr++ * istep;
  1337.             x2 = *xr++ * istep;
  1338.             XRPOW_FTOI(x1, rx1);
  1339.             x3 = *xr++ * istep;
  1340.             XRPOW_FTOI(x2, rx2);
  1341.             x4 = *xr++ * istep;
  1342.             XRPOW_FTOI(x3, rx3);
  1343.             x5 = *xr++ * istep;
  1344.             XRPOW_FTOI(x4, rx4);
  1345.             x6 = *xr++ * istep;
  1346.             XRPOW_FTOI(x5, rx5);
  1347.             x7 = *xr++ * istep;
  1348.             XRPOW_FTOI(x6, rx6);
  1349.             x8 = *xr++ * istep;
  1350.             XRPOW_FTOI(x7, rx7);
  1351.             x1 += QUANTFAC(rx1);
  1352.             XRPOW_FTOI(x8, rx8);
  1353.             x2 += QUANTFAC(rx2);
  1354.             XRPOW_FTOI(x1,*ix++);
  1355.             x3 += QUANTFAC(rx3);
  1356.             XRPOW_FTOI(x2,*ix++);
  1357.             x4 += QUANTFAC(rx4);        
  1358.             XRPOW_FTOI(x3,*ix++);
  1359.             x5 += QUANTFAC(rx5);
  1360.             XRPOW_FTOI(x4,*ix++);
  1361.             x6 += QUANTFAC(rx6);
  1362.             XRPOW_FTOI(x5,*ix++);
  1363.             x7 += QUANTFAC(rx7);
  1364.             XRPOW_FTOI(x6,*ix++);
  1365.             x8 += QUANTFAC(rx8);        
  1366.             XRPOW_FTOI(x7,*ix++);
  1367.             XRPOW_FTOI(x8,*ix++);
  1368.          }
  1369.     }
  1370. #endif
  1371. }
  1372.  
  1373.  
  1374.  
  1375.  
  1376.  
  1377.  
  1378. void quantize_xrpow_ISO( FLOAT8 xr[576], int ix[576], gr_info *cod_info )
  1379. {
  1380.   /* quantize on xr^(3/4) instead of xr */
  1381.   const FLOAT8 istep = IPOW20(cod_info->global_gain);
  1382.   
  1383. #if defined(USE_GNUC_ASM)
  1384.  
  1385.    {
  1386.       __asm__ __volatile__ (
  1387.         "\n\nloop0:\n\t"
  1388.  
  1389.         "fld" F8type " 0*" F8size "(%3)\n\t"
  1390.         "fld" F8type " 1*" F8size "(%3)\n\t"
  1391.         "fld" F8type " 2*" F8size "(%3)\n\t"
  1392.         "fld" F8type " 3*" F8size "(%3)\n\t"
  1393.  
  1394.         "addl $4*" F8size ", %3\n\t"
  1395.         "addl $16, %4\n\t"
  1396.  
  1397.         "fxch %%st(3)\n\t"
  1398.         "fmul %%st(4)\n\t"
  1399.         "fxch %%st(2)\n\t"
  1400.         "fmul %%st(4)\n\t"
  1401.         "fxch %%st(1)\n\t"
  1402.         "fmul %%st(4)\n\t"
  1403.         "fxch %%st(3)\n\t"
  1404.         "fmul %%st(4)\n\t"
  1405.  
  1406.         "dec %0\n\t"
  1407.  
  1408.         "fxch %%st(2)\n\t"
  1409.         "fadd %%st(5)\n\t"
  1410.         "fxch %%st(1)\n\t"
  1411.         "fadd %%st(5)\n\t"
  1412.         "fxch %%st(3)\n\t"
  1413.         "fadd %%st(5)\n\t"
  1414.         "fxch %%st(2)\n\t"
  1415.         "fadd %%st(5)\n\t"
  1416.  
  1417.         "fxch %%st(1)\n\t"
  1418.         "fistpl -16(%4)\n\t"
  1419.         "fxch %%st(2)\n\t"
  1420.         "fistpl -12(%4)\n\t"
  1421.         "fistpl -8(%4)\n\t"
  1422.         "fistpl -4(%4)\n\t"
  1423.  
  1424.         "jnz loop0\n\n"
  1425.  
  1426.         : /* no outputs */
  1427.         : "r" (576 / 4), "u" ((FLOAT8)(0.4054 - 0.5)), "t" (istep), "r" (xr), "r" (ix)
  1428.         : "memory", "cc"
  1429.       );
  1430.   }
  1431. #elif defined(USE_MSC_ASM)
  1432.   {
  1433.       /* asm from Acy Stapp <AStapp@austin.rr.com> */
  1434.       const FLOAT8 temp0 = 0.4054 - 0.5;
  1435.       _asm {
  1436.           mov ecx, 576/4;
  1437.           fld F8type ptr [temp0];
  1438.           fld F8type ptr [istep];
  1439.           mov eax, dword ptr [xr];
  1440.           mov edx, dword ptr [ix];
  1441.       } loop0: _asm {
  1442.           fld F8type ptr [eax+0*F8size]; // 0
  1443.           fld F8type ptr [eax+1*F8size]; // 1 0
  1444.           fld F8type ptr [eax+2*F8size]; // 2 1 0
  1445.           fld F8type ptr [eax+3*F8size]; // 3 2 1 0
  1446.  
  1447.           add eax, 4*F8size;
  1448.           add edx, 16;
  1449.  
  1450.           fxch st(3); // 0 2 1 3
  1451.           fmul st(0), st(4);
  1452.           fxch st(2); // 1 2 0 3
  1453.           fmul st(0), st(4);
  1454.           fxch st(1); // 2 1 0 3
  1455.           fmul st(0), st(4);
  1456.           fxch st(3); // 3 1 0 2
  1457.           fmul st(0), st(4);
  1458.  
  1459.           dec ecx;
  1460.  
  1461.           fxch st(2); // 0 1 3 2
  1462.           fadd st(0), st(5);
  1463.           fxch st(1); // 1 0 3 2
  1464.           fadd st(0), st(5);
  1465.           fxch st(3); // 2 0 3 1
  1466.           fadd st(0), st(5);
  1467.           fxch st(2); // 3 0 2 1
  1468.           fadd st(0), st(5);
  1469.  
  1470.           fxch st(1); // 0 3 2 1 
  1471.           fistp dword ptr [edx-16]; // 3 2 1
  1472.           fxch st(2); // 1 2 3
  1473.           fistp dword ptr [edx-12];
  1474.           fistp dword ptr [edx-8];
  1475.           fistp dword ptr [edx-4];
  1476.  
  1477.           jnz loop0;
  1478.  
  1479.           mov dword ptr [xr], eax;
  1480.           mov dword ptr [ix], edx;
  1481.           fstp st(0);
  1482.           fstp st(0);
  1483.       }
  1484.   }
  1485. #else
  1486. #if 0
  1487.    /* generic ASM */
  1488.       register int j;
  1489.       for (j=576/4;j>0;j--) {
  1490.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1491.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1492.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1493.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1494.       }
  1495. #endif
  1496.   {
  1497.       register int j;
  1498.       const FLOAT8 compareval0 = (1.0 - 0.4054)/istep;
  1499.       /* depending on architecture, it may be worth calculating a few more compareval's.
  1500.          eg.  compareval1 = (2.0 - 0.4054/istep); 
  1501.               .. and then after the first compare do this ...
  1502.               if compareval1>*xr then ix = 1;
  1503.          On a pentium166, it's only worth doing the one compare (as done here), as the second
  1504.          compare becomes more expensive than just calculating the value. Architectures with 
  1505.          slow FP operations may want to add some more comparevals. try it and send your diffs 
  1506.          statistically speaking
  1507.          73% of all xr*istep values give ix=0
  1508.          16% will give 1
  1509.          4%  will give 2
  1510.       */
  1511.       for (j=576;j>0;j--) 
  1512.         {
  1513.           if (compareval0 > *xr) {
  1514.             *(ix++) = 0;
  1515.             xr++;
  1516.           } else
  1517.         /*    *(ix++) = (int)( istep*(*(xr++))  + 0.4054); */
  1518.             XRPOW_FTOI(  istep*(*(xr++))  + ROUNDFAC , *(ix++) );
  1519.         }
  1520.   }
  1521. #endif
  1522. }
  1523.  
  1524. #endif
  1525.