开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多媒体 > 流媒体/Mpeg4/MP4 > 查看源码
quantize-pvt.c
资源名称：NETVIDEO.rar　[点击查看]
上传用户：sun1608
上传日期：2007-02-02
资源大小：6116k
文件大小：36k
源码类别：

流媒体/Mpeg4/MP4

开发平台：

Visual C++

quantize-pvt.c：源码内容
  1. #include <assert.h>
  2. #include "util.h"
  3. #include "tables.h"
  4. #include "reservoir.h"
  5. #include "quantize-pvt.h"
  7. FLOAT masking_lower=1;
  8. int convert_mdct, reduce_sidechannel;
  9. /*
  10. mt 5/99.  These global flags denote 4 possibilities:
  11.                                                                 mode    l3_xmin
  12. 1   MDCT input L/R, quantize L/R,   psy-model thresholds: L/R   -m s     either
  13. 2   MDCT input L/R, quantize M/S,   psy-model thresholds: L/R   -m j     orig
  14. 3   MDCT input M/S, quantize M/S,   psy-model thresholds: M/S   -m f     either
  15. 4   MDCT input L/R, quantize M/S,   psy-model thresholds: M/S   -m j -h  m/s
  17. 1:  convert_mdct = 0, convert_psy=0,  reduce_sidechannel=0          
  18. 2:  convert_mdct = 1, convert_psy=1,  reduce_sidechannel=1
  19. 3:  convert_mdct = 0, convert_psy=0,  reduce_sidechannel=1   (this mode no longer used)
  20. 4:  convert_mdct = 1, convert_psy=0,  reduce_sidechannel=1
  22. if (convert_mdct), then iteration_loop will quantize M/S data from
  23. the L/R input MDCT coefficients.
  25. if (convert_psy), then calc_noise will compute the noise for the L/R
  26. channels from M/S MDCT data and L/R psy-model threshold information.
  27. Distortion in ether L or R channel will be marked as distortion in
  28. both Mid and Side channels.  
  29. NOTE: 3/00: this mode has been removed.  
  31. if (reduce_sidechannel) then outer_loop will allocate less bits
  32. to the side channel and more bits to the mid channel based on relative 
  33. energies.
  34. */
  38. /*
  39.   The following table is used to implement the scalefactor
  40.   partitioning for MPEG2 as described in section
  41.   2.4.3.2 of the IS. The indexing corresponds to the
  42.   way the tables are presented in the IS:
  44.   [table_number][row_in_table][column of nr_of_sfb]
  45. */
  46. unsigned nr_of_sfb_block[6][3][4] =
  47. {
  48.   {
  49.     {6, 5, 5, 5},
  50.     {9, 9, 9, 9},
  51.     {6, 9, 9, 9}
  52.   },
  53.   {
  54.     {6, 5, 7, 3},
  55.     {9, 9, 12, 6},
  56.     {6, 9, 12, 6}
  57.   },
  58.   {
  59.     {11, 10, 0, 0},
  60.     {18, 18, 0, 0},
  61.     {15,18,0,0}
  62.   },
  63.   {
  64.     {7, 7, 7, 0},
  65.     {12, 12, 12, 0},
  66.     {6, 15, 12, 0}
  67.   },
  68.   {
  69.     {6, 6, 6, 3},
  70.     {12, 9, 9, 6},
  71.     {6, 12, 9, 6}
  72.   },
  73.   {
  74.     {8, 8, 5, 0},
  75.     {15,12,9,0},
  76.     {6,18,9,0}
  77.   }
  78. };
  81. /* Table B.6: layer3 preemphasis */
  82. int  pretab[21] =
  83. {
  84.     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  85.     1, 1, 1, 1, 2, 2, 3, 3, 3, 2
  86. };
  88. /*
  89.   Here are MPEG1 Table B.8 and MPEG2 Table B.1
  90.   -- Layer III scalefactor bands. 
  91.   Index into this using a method such as:
  92.     idx  = fr_ps->header->sampling_frequency
  93.            + (fr_ps->header->version * 3)
  94. */
  96. struct scalefac_struct sfBandIndex[6] =
  97. {
  98.   { /* Table B.2.b: 22.05 kHz */
  99.     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
  100.     {0,4,8,12,18,24,32,42,56,74,100,132,174,192}
  101.   },
  102.   { /* Table B.2.c: 24 kHz */                 /* docs: 332. mpg123: 330 */
  103.     {0,6,12,18,24,30,36,44,54,66,80,96,114,136,162,194,232,278, 332, 394,464,540,576},
  104.     {0,4,8,12,18,26,36,48,62,80,104,136,180,192}
  105.   },
  106.   { /* Table B.2.a: 16 kHz */
  107.     {0,6,12,18,24,30,36,44,54,66,80,96,116,140,168,200,238,284,336,396,464,522,576},
  108.     {0,4,8,12,18,26,36,48,62,80,104,134,174,192}
  109.   },
  110.   { /* Table B.8.b: 44.1 kHz */
  111.     {0,4,8,12,16,20,24,30,36,44,52,62,74,90,110,134,162,196,238,288,342,418,576},
  112.     {0,4,8,12,16,22,30,40,52,66,84,106,136,192}
  113.   },
  114.   { /* Table B.8.c: 48 kHz */
  115.     {0,4,8,12,16,20,24,30,36,42,50,60,72,88,106,128,156,190,230,276,330,384,576},
  116.     {0,4,8,12,16,22,28,38,50,64,80,100,126,192}
  117.   },
  118.   { /* Table B.8.a: 32 kHz */
  119.     {0,4,8,12,16,20,24,30,36,44,54,66,82,102,126,156,194,240,296,364,448,550,576},
  120.     {0,4,8,12,16,22,30,42,58,78,104,138,180,192}
  121.   }
  122. };
  124. struct scalefac_struct scalefac_band;
  127. FLOAT8 pow20[Q_MAX];
  128. FLOAT8 ipow20[Q_MAX];
  129. FLOAT8 pow43[PRECALC_SIZE];
  130. static FLOAT8 adj43[PRECALC_SIZE];
  131. static FLOAT8 adj43asm[PRECALC_SIZE];
  132. static FLOAT8 ATH_l[SBPSY_l];
  133. static FLOAT8 ATH_s[SBPSY_l];
  135. FLOAT8 ATH_mdct_long[576];
  136. FLOAT8 ATH_mdct_short[192];
  139. /************************************************************************/
  140. /*  initialization for iteration_loop */
  141. /************************************************************************/
  142. void
  143. iteration_init( lame_global_flags *gfp,III_side_info_t *l3_side, int l3_enc[2][2][576])
  144. {
  145.   gr_info *cod_info;
  146.   int ch, gr, i;
  148.   l3_side->resvDrain = 0;
  150.   if ( gfp->frameNum==0 ) {
  151.     for (i = 0; i < SBMAX_l + 1; i++) {
  152.       scalefac_band.l[i] =
  153. sfBandIndex[gfp->samplerate_index + (gfp->version * 3)].l[i];
  154.     }
  155.     for (i = 0; i < SBMAX_s + 1; i++) {
  156.       scalefac_band.s[i] =
  157. sfBandIndex[gfp->samplerate_index + (gfp->version * 3)].s[i];
  158.     }
  160.     l3_side->main_data_begin = 0;
  161.     compute_ath(gfp,ATH_l,ATH_s);
  163.     for(i=0;i<PRECALC_SIZE;i++)
  164.         pow43[i] = pow((FLOAT8)i, 4.0/3.0);
  166.     for (i = 0; i < PRECALC_SIZE-1; i++)
  167. adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
  168.     adj43[i] = 0.5;
  171.     adj43asm[0] = 0.0;
  172.     for (i = 1; i < PRECALC_SIZE; i++)
  173.       adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]),0.75);
  175.     for (i = 0; i < Q_MAX; i++) {
  176. ipow20[i] = pow(2.0, (double)(i - 210) * -0.1875);
  177. pow20[i] = pow(2.0, (double)(i - 210) * 0.25);
  178.     }
  179.   }
  182.   convert_mdct=0;
  183.   reduce_sidechannel=0;
  184.   if (gfp->mode_ext==MPG_MD_MS_LR) {
  185.     convert_mdct = 1;
  186.     reduce_sidechannel=1;
  187.   }
  188.   
  189.   /* some intializations. */
  190.   for ( gr = 0; gr < gfp->mode_gr; gr++ ){
  191.     for ( ch = 0; ch < gfp->stereo; ch++ ){
  192.       cod_info = (gr_info *) &(l3_side->gr[gr].ch[ch]);
  194.       if (cod_info->block_type == SHORT_TYPE)
  195.         {
  196.   cod_info->sfb_lmax = 0; /* No sb*/
  197.   cod_info->sfb_smax = 0;
  198.         }
  199.       else
  200. {
  201.   /* MPEG 1 doesnt use last scalefactor band */
  202.   cod_info->sfb_lmax = SBPSY_l;
  203.   cod_info->sfb_smax = SBPSY_s;    /* No sb */
  204. }
  206.     }
  207.   }
  210.   /* dont bother with scfsi. */
  211.   for ( ch = 0; ch < gfp->stereo; ch++ )
  212.     for ( i = 0; i < 4; i++ )
  213.       l3_side->scfsi[ch][i] = 0;
  214. }
  220. /* 
  221. compute the ATH for each scalefactor band 
  222. cd range:  0..96db
  224. Input:  3.3kHz signal  32767 amplitude  (3.3kHz is where ATH is smallest = -5db)
  225. longblocks:  sfb=12   en0/bw=-11db    max_en0 = 1.3db
  226. shortblocks: sfb=5           -9db              0db
  228. Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
  229. longblocks:  amp=1      sfb=12   en0/bw=-103 db      max_en0 = -92db
  230.             amp=32767   sfb=12           -12 db                 -1.4db 
  232. Input:  1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 (repeated)
  233. shortblocks: amp=1      sfb=5   en0/bw= -99                    -86 
  234.             amp=32767   sfb=5           -9  db                  4db 
  237. MAX energy of largest wave at 3.3kHz = 1db
  238. AVE energy of largest wave at 3.3kHz = -11db
  239. Let's take AVE:  -11db = maximum signal in sfb=12.  
  240. Dynamic range of CD: 96db.  Therefor energy of smallest audible wave 
  241. in sfb=12  = -11  - 96 = -107db = ATH at 3.3kHz.  
  243. ATH formula for this wave: -5db.  To adjust to LAME scaling, we need
  244. ATH = ATH_formula  - 103  (db)
  245. ATH = ATH * 2.5e-10      (ener)
  247. */
  248. FLOAT8 ATHformula(lame_global_flags *gfp,FLOAT8 f)
  249. {
  250.   FLOAT8 ath;
  251.   f  = Max(0.02, f);
  252.   /* from Painter & Spanias, 1997 */
  253.   /* minimum: (i=77) 3.3kHz = -5db */
  254.   ath=(3.640 * pow(f,-0.8)
  255.        -  6.500 * exp(-0.6*pow(f-3.3,2.0))
  256.        +  0.001 * pow(f,4.0));
  257.   /* convert to energy */
  258.   if (gfp->noATH)
  259.     ath -= 200; /* disables ATH */
  260.   else {
  261.     ath -= 114;    /* MDCT scaling.  From tests by macik and MUS420 code */
  262.     /* ath -= 109; */
  263.   }
  264. #ifdef RH_QUALITY_CONTROL 
  265.   /* purpose of RH_QUALITY_CONTROL:
  266.    * at higher quality lower ATH masking abilities   => needs more bits
  267.    * at lower quality increase ATH masking abilities => needs less bits
  268.    * works together with adjusted masking lowering of GPSYCHO thresholds
  269.    * (Robert.Hegemann@gmx.de 2000-01-30)
  270.    */
  271.   ath -= (4-gfp->VBR_q)*4.0; 
  272. #endif
  273.   ath = pow( 10.0, ath/10.0 );
  274.   return ath;
  275. }
  276.  
  278. void compute_ath(lame_global_flags *gfp,FLOAT8 ATH_l[SBPSY_l],FLOAT8 ATH_s[SBPSY_l])
  279. {
  280.   int sfb,i,start,end;
  281.   FLOAT8 ATH_f;
  282.   FLOAT8 samp_freq = gfp->out_samplerate/1000.0;
  283. #ifdef RH_ATH
  284.   /* going from average to peak level ATH masking
  285.    */
  286.   FLOAT8 adjust_mdct_scaling = 10.0; 
  287. #endif
  288.   
  290.   /* last sfb is not used */
  291.   for ( sfb = 0; sfb < SBPSY_l; sfb++ ) {
  292.     start = scalefac_band.l[ sfb ];
  293.     end   = scalefac_band.l[ sfb+1 ];
  294.     ATH_l[sfb]=1e99;
  295.     for (i=start ; i < end; i++) {
  296.       ATH_f = ATHformula(gfp,samp_freq*i/(2*576)); /* freq in kHz */
  297.       ATH_l[sfb]=Min(ATH_l[sfb],ATH_f);
  298. #ifdef RH_ATH
  299.       ATH_mdct_long[i] = ATH_f*adjust_mdct_scaling;
  300. #endif
  301.     }
  302.     /*
  303.     printf("sfb=%i %f  ATH=%f %f  %f   n",sfb,samp_freq*start/(2*576),
  304. 10*log10(ATH_l[sfb]),
  305. 10*log10( ATHformula(samp_freq*start/(2*576)))  ,
  306. 10*log10(ATHformula(samp_freq*end/(2*576))));
  307.     */
  308.   }
  310.   for ( sfb = 0; sfb < SBPSY_s; sfb++ ){
  311.     start = scalefac_band.s[ sfb ];
  312.     end   = scalefac_band.s[ sfb+1 ];
  313.     ATH_s[sfb]=1e99;
  314.     for (i=start ; i < end; i++) {
  315.       ATH_f = ATHformula(gfp,samp_freq*i/(2*192));     /* freq in kHz */
  316.       ATH_s[sfb]=Min(ATH_s[sfb],ATH_f);
  317. #ifdef RH_ATH
  318.       ATH_mdct_short[i] = ATH_f*adjust_mdct_scaling;
  319. #endif
  320.     }
  321.   }
  322. }
  328. /* convert from L/R <-> Mid/Side */
  329. void ms_convert(FLOAT8 xr[2][576],FLOAT8 xr_org[2][576])
  330. {
  331.   int i;
  332.   for ( i = 0; i < 576; i++ ) {
  333.     FLOAT8 l = xr_org[0][i];
  334.     FLOAT8 r = xr_org[1][i];
  335.     xr[0][i] = (l+r)*(SQRT2*0.5);
  336.     xr[1][i] = (l-r)*(SQRT2*0.5);
  337.   }
  338. }
  342. /************************************************************************
  343.  * allocate bits among 2 channels based on PE
  344.  * mt 6/99
  345.  ************************************************************************/
  346. void on_pe(lame_global_flags *gfp,FLOAT8 pe[2][2],III_side_info_t *l3_side,
  347. int targ_bits[2],int mean_bits, int gr)
  348. {
  349.   gr_info *cod_info;
  350.   int extra_bits,tbits,bits;
  351.   int add_bits[2]; 
  352.   int ch;
  354.   /* allocate targ_bits for granule */
  355.   ResvMaxBits( mean_bits, &tbits, &extra_bits, gr);
  356.     
  358.   for (ch=0 ; ch < gfp->stereo ; ch ++) {
  359.     /******************************************************************
  360.      * allocate bits for each channel 
  361.      ******************************************************************/
  362.     cod_info = &l3_side->gr[gr].ch[ch].tt;
  363.     
  364.     targ_bits[ch]=tbits/gfp->stereo;
  365.     
  366.     /* allocate extra bits from reservoir based on PE */
  367.     bits=0;
  368.     
  369.     /* extra bits based on PE > 700 */
  370.     add_bits[ch]=(pe[gr][ch]-750)/1.55;  /* 1.4; */
  371.     
  372.     /* short blocks need extra, no matter what the pe */
  373.     if (cod_info->block_type==SHORT_TYPE) 
  374.       if (add_bits[ch]<500) add_bits[ch]=500;
  375.     
  376.     if (add_bits[ch] < 0) add_bits[ch]=0;
  377.     bits += add_bits[ch];
  378.     
  379.     if (bits > extra_bits) add_bits[ch] = (extra_bits*add_bits[ch])/bits;
  380.     if ((targ_bits[ch]+add_bits[ch]) > 4095) 
  381.       add_bits[ch]=4095-targ_bits[ch];
  383.     targ_bits[ch] = targ_bits[ch] + add_bits[ch];
  384.     extra_bits -= add_bits[ch];
  385.   }
  386. }
  388. void reduce_side(int targ_bits[2],FLOAT8 ms_ener_ratio,int mean_bits)
  389. {
  390. int ch;
  391. int numchn=2;
  392.     /*  ms_ener_ratio = 0:  allocate 66/33  mid/side  fac=.33  
  393.      *  ms_ener_ratio =.5:  allocate 50/50 mid/side   fac= 0 */
  394.     /* 75/25 split is fac=.5 */
  395.     /* float fac = .50*(.5-ms_ener_ratio[gr])/.5;*/
  396.     float fac = .33*(.5-ms_ener_ratio)/.5;
  397.     if (fac<0) fac=0;
  398.     
  399.     if (targ_bits[1] >= 125) {
  400.       /* dont reduce side channel below 125 bits */
  401.       if (targ_bits[1]-targ_bits[1]*fac > 125) {
  402. targ_bits[0] += targ_bits[1]*fac;
  403. targ_bits[1] -= targ_bits[1]*fac;
  404.       } else {
  405. targ_bits[0] += targ_bits[1] - 125;
  406. targ_bits[1] = 125;
  407.       }
  408.     }
  409.     
  410.     /* dont allow to many bits per channel */  
  411.     for (ch=0; ch<numchn; ch++) {
  412.       int max_bits = Min(4095,mean_bits/2 + 1200);
  413.       if (targ_bits[ch] > max_bits) {
  414. targ_bits[ch] = max_bits;
  415.       }
  416.     }
  418. }
  420. /*************************************************************************** 
  421.  *         inner_loop                                                      * 
  422.  *************************************************************************** 
  423.  * The code selects the best global gain for a particular set of scalefacs */
  424.  
  425. int
  426. inner_loop( lame_global_flags *gfp,FLOAT8 xrpow[576],
  427.     int l3_enc[576], int max_bits,
  428.     gr_info *cod_info)
  429. {
  430.     int bits;
  431.     assert( max_bits >= 0 );
  432.     cod_info->global_gain--;
  433.     do
  434.     {
  435.       cod_info->global_gain++;
  436.       bits = count_bits(gfp,l3_enc, xrpow, cod_info);
  437.     }
  438.     while ( bits > max_bits );
  439.     return bits;
  440. }
  444. /*************************************************************************/
  445. /*            scale_bitcount                                             */
  446. /*************************************************************************/
  448. /* Also calculates the number of bits necessary to code the scalefactors. */
  450. int scale_bitcount( III_scalefac_t *scalefac, gr_info *cod_info)
  451. {
  452.     int i, k, sfb, max_slen1 = 0, max_slen2 = 0, /*a, b, */ ep = 2;
  454.     static int slen1[16] = { 1, 1, 1, 1, 8, 2, 2, 2, 4, 4, 4, 8, 8, 8,16,16 };
  455.     static int slen2[16] = { 1, 2, 4, 8, 1, 2, 4, 8, 2, 4, 8, 2, 4, 8, 4, 8 };
  457.     static int slen1_tab[16] = {0,
  458. 18, 36, 54, 54, 36, 54, 72, 54, 72, 90, 72, 90,108,108,126
  459.     };
  460.     static int slen2_tab[16] = {0,
  461. 10, 20, 30, 33, 21, 31, 41, 32, 42, 52, 43, 53, 63, 64, 74
  462.     };
  463.     int *tab;
  466.     if ( cod_info->block_type == SHORT_TYPE )
  467.     {
  468.             tab = slen1_tab;
  469.             /* a = 18; b = 18;  */
  470.             for ( i = 0; i < 3; i++ )
  471.             {
  472.                 for ( sfb = 0; sfb < 6; sfb++ )
  473.                     if (scalefac->s[sfb][i] > max_slen1 )
  474.                         max_slen1 = scalefac->s[sfb][i];
  475.                 for (sfb = 6; sfb < SBPSY_s; sfb++ )
  476.                     if ( scalefac->s[sfb][i] > max_slen2 )
  477.                         max_slen2 = scalefac->s[sfb][i];
  478.             }
  479.     }
  480.     else
  481.     { /* block_type == 1,2,or 3 */
  482.         tab = slen2_tab;
  483.         /* a = 11; b = 10;   */
  484.         for ( sfb = 0; sfb < 11; sfb++ )
  485.             if ( scalefac->l[sfb] > max_slen1 )
  486.                 max_slen1 = scalefac->l[sfb];
  488. if (!cod_info->preflag) {
  489.     for ( sfb = 11; sfb < SBPSY_l; sfb++ )
  490. if (scalefac->l[sfb] < pretab[sfb])
  491.     break;
  493.     if (sfb == SBPSY_l) {
  494. cod_info->preflag = 1;
  495. for ( sfb = 11; sfb < SBPSY_l; sfb++ )
  496.     scalefac->l[sfb] -= pretab[sfb];
  497.     }
  498. }
  500.         for ( sfb = 11; sfb < SBPSY_l; sfb++ )
  501.             if ( scalefac->l[sfb] > max_slen2 )
  502.                 max_slen2 = scalefac->l[sfb];
  503.     }
  507.     /* from Takehiro TOMINAGA <tominaga@isoternet.org> 10/99
  508.      * loop over *all* posible values of scalefac_compress to find the
  509.      * one which uses the smallest number of bits.  ISO would stop
  510.      * at first valid index */
  511.     cod_info->part2_length = LARGE_BITS;
  512.     for ( k = 0; k < 16; k++ )
  513.     {
  514.         if ( (max_slen1 < slen1[k]) && (max_slen2 < slen2[k]) &&
  515.              ((int)cod_info->part2_length > tab[k])) {
  516.   cod_info->part2_length=tab[k];
  517.   cod_info->scalefac_compress=k;
  518.   ep=0;  /* we found a suitable scalefac_compress */
  519. }
  520.     }
  521.     return ep;
  522. }
  526. /*
  527.   table of largest scalefactors (number of bits) for MPEG2
  528. */
  529. /*
  530. static unsigned max_sfac_tab[6][4] =
  531. {
  532.     {4, 4, 3, 3},
  533.     {4, 4, 3, 0},
  534.     {3, 2, 0, 0},
  535.     {4, 5, 5, 0},
  536.     {3, 3, 3, 0},
  537.     {2, 2, 0, 0}
  538. };
  539. */
  540. /*
  541.   table of largest scalefactor values for MPEG2
  542. */
  543. static unsigned max_range_sfac_tab[6][4] =
  544. {
  545.  { 15, 15, 7,  7},
  546.  { 15, 15, 7,  0},
  547.  { 7,  3,  0,  0},
  548.  { 15, 31, 31, 0},
  549.  { 7,  7,  7,  0},
  550.  { 3,  3,  0,  0}
  551. };
  557. /*************************************************************************/
  558. /*            scale_bitcount_lsf                                         */
  559. /*************************************************************************/
  561. /* Also counts the number of bits to encode the scalefacs but for MPEG 2 */ 
  562. /* Lower sampling frequencies  (24, 22.05 and 16 kHz.)                   */
  563.  
  564. /*  This is reverse-engineered from section 2.4.3.2 of the MPEG2 IS,     */
  565. /* "Audio Decoding Layer III"                                            */
  567. int scale_bitcount_lsf(III_scalefac_t *scalefac, gr_info *cod_info)
  568. {
  569.     int table_number, row_in_table, partition, nr_sfb, window, over;
  570.     int i, sfb, max_sfac[ 4 ];
  571.     unsigned *partition_table;
  573.     /*
  574.       Set partition table. Note that should try to use table one,
  575.       but do not yet...
  576.     */
  577.     if ( cod_info->preflag )
  578. table_number = 2;
  579.     else
  580. table_number = 0;
  582.     for ( i = 0; i < 4; i++ )
  583. max_sfac[i] = 0;
  585.     if ( cod_info->block_type == SHORT_TYPE )
  586.     {
  587.     row_in_table = 1;
  588.     partition_table = &nr_of_sfb_block[table_number][row_in_table][0];
  589.     for ( sfb = 0, partition = 0; partition < 4; partition++ )
  590.     {
  591. nr_sfb = partition_table[ partition ] / 3;
  592. for ( i = 0; i < nr_sfb; i++, sfb++ )
  593.     for ( window = 0; window < 3; window++ )
  594. if ( scalefac->s[sfb][window] > max_sfac[partition] )
  595.     max_sfac[partition] = scalefac->s[sfb][window];
  596.     }
  597.     }
  598.     else
  599.     {
  600. row_in_table = 0;
  601. partition_table = &nr_of_sfb_block[table_number][row_in_table][0];
  602. for ( sfb = 0, partition = 0; partition < 4; partition++ )
  603. {
  604.     nr_sfb = partition_table[ partition ];
  605.     for ( i = 0; i < nr_sfb; i++, sfb++ )
  606. if ( scalefac->l[sfb] > max_sfac[partition] )
  607.     max_sfac[partition] = scalefac->l[sfb];
  608. }
  609.     }
  611.     for ( over = 0, partition = 0; partition < 4; partition++ )
  612.     {
  613. if ( max_sfac[partition] > (int)max_range_sfac_tab[table_number][partition] )
  614.     over++;
  615.     }
  616.     if ( !over )
  617.     {
  618. /*
  619.   Since no bands have been over-amplified, we can set scalefac_compress
  620.   and slen[] for the formatter
  621. */
  622. static int log2tab[] = { 0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4 };
  624. unsigned slen1, slen2, slen3, slen4;
  626.         cod_info->sfb_partition_table = &nr_of_sfb_block[table_number][row_in_table][0];
  627. for ( partition = 0; partition < 4; partition++ )
  628.     cod_info->slen[partition] = log2tab[max_sfac[partition]];
  630. /* set scalefac_compress */
  631. slen1 = cod_info->slen[ 0 ];
  632. slen2 = cod_info->slen[ 1 ];
  633. slen3 = cod_info->slen[ 2 ];
  634. slen4 = cod_info->slen[ 3 ];
  636. switch ( table_number )
  637. {
  638.   case 0:
  639.     cod_info->scalefac_compress = (((slen1 * 5) + slen2) << 4)
  640. + (slen3 << 2)
  641. + slen4;
  642.     break;
  644.   case 1:
  645.     cod_info->scalefac_compress = 400
  646. + (((slen1 * 5) + slen2) << 2)
  647. + slen3;
  648.     break;
  650.   case 2:
  651.     cod_info->scalefac_compress = 500 + (slen1 * 3) + slen2;
  652.     break;
  654.   default:
  655.     fprintf( stderr, "intensity stereo not implemented yetn" );
  656.     exit( EXIT_FAILURE );
  657.     break;
  658. }
  659.     }
  660. #ifdef DEBUG
  661.     if ( over ) 
  662.         printf( "---WARNING !! Amplification of some bands over limitsn" );
  663. #endif
  664.     if (!over) {
  665.       assert( cod_info->sfb_partition_table );     
  666.       cod_info->part2_length=0;
  667.       for ( partition = 0; partition < 4; partition++ )
  668. cod_info->part2_length += cod_info->slen[partition] * cod_info->sfb_partition_table[partition];
  669.     }
  670.     return over;
  671. }
  677. /*************************************************************************/
  678. /*            calc_xmin                                                  */
  679. /*************************************************************************/
  681. /*
  682.   Calculate the allowed distortion for each scalefactor band,
  683.   as determined by the psychoacoustic model.
  684.   xmin(sb) = ratio(sb) * en(sb) / bw(sb)
  686.   returns number of sfb's with energy > ATH
  687. */
  688. int calc_xmin( lame_global_flags *gfp,FLOAT8 xr[576], III_psy_ratio *ratio,
  689.        gr_info *cod_info, III_psy_xmin *l3_xmin)
  690. {
  691.     int start, end, bw,l, b, ath_over=0;
  692. u_int sfb;
  693.     FLOAT8 en0, xmin, ener;
  695.     if (gfp->ATHonly) {    
  696.       for ( sfb = cod_info->sfb_smax; sfb < SBPSY_s; sfb++ )
  697.   for ( b = 0; b < 3; b++ )
  698.       l3_xmin->s[sfb][b]=ATH_s[sfb];
  699.       for ( sfb = 0; sfb < cod_info->sfb_lmax; sfb++ )
  700.   l3_xmin->l[sfb]=ATH_l[sfb];
  702.     }else{
  704.       for ( sfb = cod_info->sfb_smax; sfb < SBPSY_s; sfb++ ) {
  705. start = scalefac_band.s[ sfb ];
  706.         end   = scalefac_band.s[ sfb + 1 ];
  707. bw = end - start;
  708.         for ( b = 0; b < 3; b++ ) {
  709.   for (en0 = 0.0, l = start; l < end; l++) {
  710.     ener = xr[l * 3 + b];
  711.     ener = ener * ener;
  712.     en0 += ener;
  713.   }
  714.   en0 /= bw;
  716.   xmin = ratio->en.s[sfb][b];
  717.   if (xmin > 0.0)
  718.     xmin = en0 * ratio->thm.s[sfb][b] * masking_lower / xmin;
  720. #ifdef RH_ATH
  721.           /* do not mix up ATH masking with GPSYCHO thresholds
  722.    */
  723.   l3_xmin->s[sfb][b] = Max(1e-20, xmin);
  724. #else
  725.   l3_xmin->s[sfb][b] = Max(ATH_s[sfb], xmin);
  726. #endif
  727.   if (en0 > ATH_s[sfb]) ath_over++;
  728. }
  729.       }
  731.       for ( sfb = 0; sfb < cod_info->sfb_lmax; sfb++ ){
  732. start = scalefac_band.l[ sfb ];
  733. end   = scalefac_band.l[ sfb+1 ];
  734. bw = end - start;
  736.         for (en0 = 0.0, l = start; l < end; l++ ) {
  737.   ener = xr[l] * xr[l];
  738.   en0 += ener;
  739. }
  740. en0 /= bw;
  742. xmin = ratio->en.l[sfb];
  743. if (xmin > 0.0)
  744.   xmin = en0 * ratio->thm.l[sfb] * masking_lower / xmin;
  747. #ifdef RH_ATH
  748.         /* do not mix up ATH masking with GPSYCHO thresholds
  749.  */
  750. l3_xmin->l[sfb]=Max(1e-20, xmin);
  751. #else
  752. l3_xmin->l[sfb]=Max(ATH_l[sfb], xmin);
  753. #endif
  754. if (en0 > ATH_l[sfb]) ath_over++;
  755.       }
  756.     }
  757.     return ath_over;
  758. }
  762. /*************************************************************************/
  763. /*            loop_break                                                 */
  764. /*************************************************************************/
  766. /*  Function: Returns zero if there is a scalefac which has not been
  767.     amplified. Otherwise it returns one. 
  768. */
  770. int loop_break( III_scalefac_t *scalefac, gr_info *cod_info)
  771. {
  772.     int i;
  773. u_int sfb;
  775.     for ( sfb = 0; sfb < cod_info->sfb_lmax; sfb++ )
  776.         if ( scalefac->l[sfb] == 0 )
  777.     return 0;
  779.     for ( sfb = cod_info->sfb_smax; sfb < SBPSY_s; sfb++ )
  780.       for ( i = 0; i < 3; i++ ) 
  781.             if ( scalefac->s[sfb][i] == 0 )
  782. return 0;
  784.     return 1;
  785. }
  799. /*
  800.  ----------------------------------------------------------------------
  801.   if someone wants to try to find a faster step search function,
  802.   here is some code which gives a lower bound for the step size:
  803.   
  804.   for (max_xrspow = 0, i = 0; i < 576; ++i)
  805.   {
  806.     max_xrspow = Max(max_xrspow, xrspow[i]);
  807.   }
  808.   lowerbound = 210+log10(max_xrspow/IXMAX_VAL)/(0.1875*LOG2);
  809.  
  810.  
  811.                                                  Robert.Hegemann@gmx.de
  812.  ----------------------------------------------------------------------
  813. */
  816. typedef enum {
  817.     BINSEARCH_NONE,
  818.     BINSEARCH_UP, 
  819.     BINSEARCH_DOWN
  820. } binsearchDirection_t;
  822. /*-------------------------------------------------------------------------*/
  823. int 
  824. bin_search_StepSize2 (lame_global_flags *gfp,int desired_rate, int start, int *ix, 
  825.                       FLOAT8 xrspow[576], gr_info *cod_info)
  826. /*-------------------------------------------------------------------------*/
  827. {
  828.     static int CurrentStep = 4;
  829.     int nBits;
  830.     int flag_GoneOver = 0;
  831.     int StepSize = start;
  832.     binsearchDirection_t Direction = BINSEARCH_NONE;
  834.     do
  835.     {
  836. cod_info->global_gain = StepSize;
  837. nBits = count_bits(gfp,ix, xrspow, cod_info);  
  839. if (CurrentStep == 1 )
  840.         {
  841.     break; /* nothing to adjust anymore */
  842. }
  843. if (flag_GoneOver)
  844. {
  845.     CurrentStep /= 2;
  846. }
  847. if (nBits > desired_rate)  /* increase Quantize_StepSize */
  848. {
  849.     if (Direction == BINSEARCH_DOWN && !flag_GoneOver)
  850.     {
  851. flag_GoneOver = 1;
  852. CurrentStep /= 2; /* late adjust */
  853.     }
  854.     Direction = BINSEARCH_UP;
  855.     StepSize += CurrentStep;
  856.     if (StepSize > 255) break;
  857. }
  858. else if (nBits < desired_rate)
  859. {
  860.     if (Direction == BINSEARCH_UP && !flag_GoneOver)
  861.     {
  862. flag_GoneOver = 1;
  863. CurrentStep /= 2; /* late adjust */
  864.     }
  865.     Direction = BINSEARCH_DOWN;
  866.     StepSize -= CurrentStep;
  867.     if (StepSize < 0) break;
  868. }
  869. else break; /* nBits == desired_rate;; most unlikely to happen.*/
  870.     } while (1); /* For-ever, break is adjusted. */
  872.     CurrentStep = abs(start - StepSize);
  873.     
  874.     if (CurrentStep >= 4) {
  875. CurrentStep = 4;
  876.     } else {
  877. CurrentStep = 2;
  878.     }
  880.     return nBits;
  881. }
  890. #if (defined(__GNUC__) && defined(__i386__))
  891. #define USE_GNUC_ASM
  892. #endif
  893. #ifdef _MSC_VER
  894. #define USE_MSC_ASM
  895. #endif
  899. /*********************************************************************
  900.  * XRPOW_FTOI is a macro to convert floats to ints.  
  901.  * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x]
  902.  *                                         ROUNDFAC= -0.0946
  903.  *
  904.  * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]   
  905.  *                                   ROUNDFAC=0.4054
  906.  *********************************************************************/
  907. #ifdef USE_GNUC_ASM
  908. #  define QUANTFAC(rx)  adj43asm[rx]
  909. #  define ROUNDFAC -0.0946
  910. #  define XRPOW_FTOI(src, dest) 
  911.      asm ("fistpl %0 " : "=m"(dest) : "t"(src) : "st")
  912. #elif defined (USE_MSC_ASM)
  913. #  define QUANTFAC(rx)  adj43asm[rx]
  914. #  define ROUNDFAC -0.0946
  915. #  define XRPOW_FTOI(src, dest) do { 
  916.      FLOAT8 src_ = (src); 
  917.      int dest_; 
  918.      { 
  919.        __asm fld src_ 
  920.        __asm fistp dest_ 
  921.      } 
  922.      (dest) = dest_; 
  923.    } while (0)
  924. #else
  925. #  define QUANTFAC(rx)  adj43[rx]
  926. #  define ROUNDFAC 0.4054
  927. #  define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
  928. #endif
  930. #ifdef USE_MSC_ASM
  931. /* define F8type and F8size according to type of FLOAT8 */
  932. # if defined FLOAT8_is_double
  933. #  define F8type qword
  934. #  define F8size 8
  935. # elif defined FLOAT8_is_float
  936. #  define F8type dword
  937. #  define F8size 4
  938. # else
  939. /* only float and double supported */
  940. #  error invalid FLOAT8 type for USE_MSC_ASM
  941. # endif
  942. #endif
  944. #ifdef USE_GNUC_ASM
  945. /* define F8type and F8size according to type of FLOAT8 */
  946. # if defined FLOAT8_is_double
  947. #  define F8type "l"
  948. #  define F8size "8"
  949. # elif defined FLOAT8_is_float
  950. #  define F8type "s"
  951. #  define F8size "4"
  952. # else
  953. /* only float and double supported */
  954. #  error invalid FLOAT8 type for USE_GNUC_ASM
  955. # endif
  956. #endif
  958. /*********************************************************************
  959.  * nonlinear quantization of xr 
  960.  * More accurate formula than the ISO formula.  Takes into account
  961.  * the fact that we are quantizing xr -> ix, but we want ix^4/3 to be 
  962.  * as close as possible to x^4/3.  (taking the nearest int would mean
  963.  * ix is as close as possible to xr, which is different.)
  964.  * From Segher Boessenkool <segher@eastsite.nl>  11/1999
  965.  * ASM optimization from 
  966.  *    Mathew Hendry <scampi@dial.pipex.com> 11/1999
  967.  *    Acy Stapp <AStapp@austin.rr.com> 11/1999
  968.  *    Takehiro Tominaga <tominaga@isoternet.org> 11/1999
  969.  *********************************************************************/
  971. void quantize_xrpow(FLOAT8 xr[576], int ix[576], gr_info *cod_info) {
  972.   /* quantize on xr^(3/4) instead of xr */
  973.   const FLOAT8 istep = IPOW20(cod_info->global_gain);
  975. #if defined (USE_GNUC_ASM) 
  976.   {
  977.       int rx[4];
  978.       __asm__ __volatile__(
  979.         "nnloop1:nt"
  981.         "fld" F8type " 0*" F8size "(%1)nt"
  982.         "fld" F8type " 1*" F8size "(%1)nt"
  983.         "fld" F8type " 2*" F8size "(%1)nt"
  984.         "fld" F8type " 3*" F8size "(%1)nt"
  986.         "fxch %%st(3)nt"
  987.         "fmul %%st(4)nt"
  988.         "fxch %%st(2)nt"
  989.         "fmul %%st(4)nt"
  990.         "fxch %%st(1)nt"
  991.         "fmul %%st(4)nt"
  992.         "fxch %%st(3)nt"
  993.         "fmul %%st(4)nt"
  995.         "addl $4*" F8size ", %1nt"
  996.         "addl $16, %3nt"
  998.         "fxch %%st(2)nt"
  999.         "fistl %5nt"
  1000.         "fxch %%st(1)nt"
  1001.         "fistl 4+%5nt"
  1002.         "fxch %%st(3)nt"
  1003.         "fistl 8+%5nt"
  1004.         "fxch %%st(2)nt"
  1005.         "fistl 12+%5nt"
  1007.         "dec %4nt"
  1009.         "movl %5, %%eaxnt"
  1010.         "movl 4+%5, %%ebxnt"
  1011.         "fxch %%st(1)nt"
  1012.         "fadd" F8type " (%2,%%eax," F8size ")nt"
  1013.         "fxch %%st(3)nt"
  1014.         "fadd" F8type " (%2,%%ebx," F8size ")nt"
  1016.         "movl 8+%5, %%eaxnt"
  1017.         "movl 12+%5, %%ebxnt"
  1018.         "fxch %%st(2)nt"
  1019.         "fadd" F8type " (%2,%%eax," F8size ")nt"
  1020.         "fxch %%st(1)nt"
  1021.         "fadd" F8type " (%2,%%ebx," F8size ")nt"
  1023.         "fxch %%st(3)nt"
  1024.         "fistpl -16(%3)nt"
  1025.         "fxch %%st(1)nt"
  1026.         "fistpl -12(%3)nt"
  1027.         "fistpl -8(%3)nt"
  1028.         "fistpl -4(%3)nt"
  1030.         "jnz loop1nn"
  1031.         : /* no outputs */
  1032.         : "t" (istep), "r" (xr), "r" (adj43asm), "r" (ix), "r" (576 / 4), "m" (rx)
  1033.         : "%eax", "%ebx", "memory", "cc"
  1034.       );
  1035.   }
  1036. #elif defined (USE_MSC_ASM)
  1037.   {
  1038.       /* asm from Acy Stapp <AStapp@austin.rr.com> */
  1039.       int rx[4];
  1040.       _asm {
  1041.           fld F8type ptr [istep]
  1042.           mov esi, dword ptr [xr]
  1043.           lea edi, dword ptr [adj43asm]
  1044.           mov edx, dword ptr [ix]
  1045.           mov ecx, 576/4
  1046.       } loop1: _asm {
  1047.           fld F8type ptr [esi+(0*F8size)] // 0
  1048.           fld F8type ptr [esi+(1*F8size)] // 1 0
  1049.           fld F8type ptr [esi+(2*F8size)] // 2 1 0
  1050.           fld F8type ptr [esi+(3*F8size)] // 3 2 1 0
  1051.           fxch st(3)                  // 0 2 1 3
  1052.           fmul st(0), st(4)
  1053.           fxch st(2)                  // 1 2 0 3
  1054.           fmul st(0), st(4)
  1055.           fxch st(1)                  // 2 1 0 3
  1056.           fmul st(0), st(4)
  1057.           fxch st(3)                  // 3 1 0 2
  1058.           fmul st(0), st(4)
  1060.           add esi, 4*F8size
  1061.           add edx, 16
  1063.           fxch st(2)                  // 0 1 3 2
  1064.           fist dword ptr [rx]
  1065.           fxch st(1)                  // 1 0 3 2
  1066.           fist dword ptr [rx+4]
  1067.           fxch st(3)                  // 2 0 3 1
  1068.           fist dword ptr [rx+8]
  1069.           fxch st(2)                  // 3 0 2 1
  1070.           fist dword ptr [rx+12]
  1072.           dec ecx
  1074.           mov eax, dword ptr [rx]
  1075.           mov ebx, dword ptr [rx+4]
  1076.           fxch st(1)                  // 0 3 2 1
  1077.           fadd F8type ptr [edi+eax*F8size]
  1078.           fxch st(3)                  // 1 3 2 0
  1079.           fadd F8type ptr [edi+ebx*F8size]
  1081.           mov eax, dword ptr [rx+8]
  1082.           mov ebx, dword ptr [rx+12]
  1083.           fxch st(2)                  // 2 3 1 0
  1084.           fadd F8type ptr [edi+eax*F8size]
  1085.           fxch st(1)                  // 3 2 1 0
  1086.           fadd F8type ptr [edi+ebx*F8size]
  1087.           fxch st(3)                  // 0 2 1 3
  1088.           fistp dword ptr [edx-16]    // 2 1 3
  1089.           fxch st(1)                  // 1 2 3
  1090.           fistp dword ptr [edx-12]    // 2 3
  1091.           fistp dword ptr [edx-8]     // 3
  1092.           fistp dword ptr [edx-4]
  1094.           jnz loop1
  1096.           mov dword ptr [xr], esi
  1097.           mov dword ptr [ix], edx
  1098.           fstp st(0)
  1099.       }
  1100.   }
  1101. #else
  1102. #if 0
  1103.   {   /* generic code if you write ASM for XRPOW_FTOI() */
  1104.       FLOAT8 x;
  1105.       int j, rx;
  1106.       for (j = 576 / 4; j > 0; --j) {
  1107.           x = *xr++ * istep;
  1108.           XRPOW_FTOI(x, rx);
  1109.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1111.           x = *xr++ * istep;
  1112.           XRPOW_FTOI(x, rx);
  1113.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1115.           x = *xr++ * istep;
  1116.           XRPOW_FTOI(x, rx);
  1117.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1119.           x = *xr++ * istep;
  1120.           XRPOW_FTOI(x, rx);
  1121.           XRPOW_FTOI(x + QUANTFAC(rx), *ix++);
  1122.       }
  1123.   }
  1124. #endif
  1125.   {/* from Wilfried.Behne@t-online.de.  Reported to be 2x faster than 
  1126.       the above code (when not using ASM) on PowerPC */
  1127.       int j;
  1128.      
  1129.       for ( j = 576/8; j > 0; --j)
  1130.       {
  1131. FLOAT8 x1, x2, x3, x4, x5, x6, x7, x8;
  1132. int rx1, rx2, rx3, rx4, rx5, rx6, rx7, rx8;
  1133. x1 = *xr++ * istep;
  1134. x2 = *xr++ * istep;
  1135. XRPOW_FTOI(x1, rx1);
  1136. x3 = *xr++ * istep;
  1137. XRPOW_FTOI(x2, rx2);
  1138. x4 = *xr++ * istep;
  1139. XRPOW_FTOI(x3, rx3);
  1140. x5 = *xr++ * istep;
  1141. XRPOW_FTOI(x4, rx4);
  1142. x6 = *xr++ * istep;
  1143. XRPOW_FTOI(x5, rx5);
  1144. x7 = *xr++ * istep;
  1145. XRPOW_FTOI(x6, rx6);
  1146. x8 = *xr++ * istep;
  1147. XRPOW_FTOI(x7, rx7);
  1148. x1 += QUANTFAC(rx1);
  1149. XRPOW_FTOI(x8, rx8);
  1150. x2 += QUANTFAC(rx2);
  1151. XRPOW_FTOI(x1,*ix++);
  1152. x3 += QUANTFAC(rx3);
  1153. XRPOW_FTOI(x2,*ix++);
  1154. x4 += QUANTFAC(rx4);
  1155. XRPOW_FTOI(x3,*ix++);
  1156. x5 += QUANTFAC(rx5);
  1157. XRPOW_FTOI(x4,*ix++);
  1158. x6 += QUANTFAC(rx6);
  1159. XRPOW_FTOI(x5,*ix++);
  1160. x7 += QUANTFAC(rx7);
  1161. XRPOW_FTOI(x6,*ix++);
  1162. x8 += QUANTFAC(rx8);
  1163. XRPOW_FTOI(x7,*ix++);
  1164. XRPOW_FTOI(x8,*ix++);
  1165.       }
  1166. }
  1167. #endif
  1168. }
  1175. void quantize_xrpow_ISO( FLOAT8 xr[576], int ix[576], gr_info *cod_info )
  1176. {
  1177.   /* quantize on xr^(3/4) instead of xr */
  1178.   const FLOAT8 istep = IPOW20(cod_info->global_gain);
  1179.   
  1180. #if defined(USE_GNUC_ASM)
  1181.    {
  1182.       __asm__ __volatile__ (
  1183.         "nnloop0:nt"
  1185.         "fld" F8type " 0*" F8size "(%3)nt"
  1186.         "fld" F8type " 1*" F8size "(%3)nt"
  1187.         "fld" F8type " 2*" F8size "(%3)nt"
  1188.         "fld" F8type " 3*" F8size "(%3)nt"
  1190.         "addl $4*" F8size ", %3nt"
  1191.         "addl $16, %4nt"
  1193.         "fxch %%st(3)nt"
  1194.         "fmul %%st(4)nt"
  1195.         "fxch %%st(2)nt"
  1196.         "fmul %%st(4)nt"
  1197.         "fxch %%st(1)nt"
  1198.         "fmul %%st(4)nt"
  1199.         "fxch %%st(3)nt"
  1200.         "fmul %%st(4)nt"
  1202.         "dec %0nt"
  1204.         "fxch %%st(2)nt"
  1205.         "fadd %%st(5)nt"
  1206.         "fxch %%st(1)nt"
  1207.         "fadd %%st(5)nt"
  1208.         "fxch %%st(3)nt"
  1209.         "fadd %%st(5)nt"
  1210.         "fxch %%st(2)nt"
  1211.         "fadd %%st(5)nt"
  1213.         "fxch %%st(1)nt"
  1214.         "fistpl -16(%4)nt"
  1215.         "fxch %%st(2)nt"
  1216.         "fistpl -12(%4)nt"
  1217.         "fistpl -8(%4)nt"
  1218.         "fistpl -4(%4)nt"
  1220.         "jnz loop0nn"
  1222.         : /* no outputs */
  1223.         : "r" (576 / 4), "u" ((FLOAT8)(0.4054 - 0.5)), "t" (istep), "r" (xr), "r" (ix)
  1224.         : "memory", "cc"
  1225.       );
  1226.   }
  1227. #elif defined(USE_MSC_ASM)
  1228.   {
  1229.       /* asm from Acy Stapp <AStapp@austin.rr.com> */
  1230.       const FLOAT8 temp0 = 0.4054 - 0.5;
  1231.       _asm {
  1232.           mov ecx, 576/4;
  1233.           fld F8type ptr [temp0];
  1234.           fld F8type ptr [istep];
  1235.           mov eax, dword ptr [xr];
  1236.           mov edx, dword ptr [ix];
  1237.       } loop0: _asm {
  1238.           fld F8type ptr [eax+0*F8size]; // 0
  1239.           fld F8type ptr [eax+1*F8size]; // 1 0
  1240.           fld F8type ptr [eax+2*F8size]; // 2 1 0
  1241.           fld F8type ptr [eax+3*F8size]; // 3 2 1 0
  1243.           add eax, 4*F8size;
  1244.           add edx, 16;
  1246.           fxch st(3); // 0 2 1 3
  1247.           fmul st(0), st(4);
  1248.           fxch st(2); // 1 2 0 3
  1249.           fmul st(0), st(4);
  1250.           fxch st(1); // 2 1 0 3
  1251.           fmul st(0), st(4);
  1252.           fxch st(3); // 3 1 0 2
  1253.           fmul st(0), st(4);
  1255.           dec ecx;
  1257.           fxch st(2); // 0 1 3 2
  1258.           fadd st(0), st(5);
  1259.           fxch st(1); // 1 0 3 2
  1260.           fadd st(0), st(5);
  1261.           fxch st(3); // 2 0 3 1
  1262.           fadd st(0), st(5);
  1263.           fxch st(2); // 3 0 2 1
  1264.           fadd st(0), st(5);
  1266.           fxch st(1); // 0 3 2 1 
  1267.           fistp dword ptr [edx-16]; // 3 2 1
  1268.           fxch st(2); // 1 2 3
  1269.           fistp dword ptr [edx-12];
  1270.           fistp dword ptr [edx-8];
  1271.           fistp dword ptr [edx-4];
  1273.           jnz loop0;
  1275.           mov dword ptr [xr], eax;
  1276.           mov dword ptr [ix], edx;
  1277.           fstp st(0);
  1278.           fstp st(0);
  1279.       }
  1280.   }
  1281. #else
  1282. #if 0
  1283.    /* generic ASM */
  1284.       register int j;
  1285.       for (j=576/4;j>0;j--) {
  1286.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1287.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1288.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1289.          XRPOW_FTOI(istep * (*xr++) + ROUNDFAC, *ix++);
  1290.       }
  1291. #endif
  1292.   {
  1293.       register int j;
  1294.       const FLOAT8 compareval0 = (1.0 - 0.4054)/istep;
  1295.       /* depending on architecture, it may be worth calculating a few more compareval's.
  1296.          eg.  compareval1 = (2.0 - 0.4054/istep); 
  1297.               .. and then after the first compare do this ...
  1298.               if compareval1>*xr then ix = 1;
  1299.          On a pentium166, it's only worth doing the one compare (as done here), as the second
  1300.          compare becomes more expensive than just calculating the value. Architectures with 
  1301.          slow FP operations may want to add some more comparevals. try it and send your diffs 
  1302.          statistically speaking
  1303.          73% of all xr*istep values give ix=0
  1304.          16% will give 1
  1305.          4%  will give 2
  1306.       */
  1307.       for (j=576;j>0;j--) 
  1308.         {
  1309.           if (compareval0 > *xr) {
  1310.             *(ix++) = 0;
  1311.             xr++;
  1312.           } else
  1313.     /*    *(ix++) = (int)( istep*(*(xr++))  + 0.4054); */
  1314.             XRPOW_FTOI(  istep*(*(xr++))  + ROUNDFAC , *(ix++) );
  1315.         }
  1316.   }
  1317. #endif
  1318. }