quantize.c
上传用户:sun1608
上传日期:2007-02-02
资源大小:6116k
文件大小:36k
源码类别:

流媒体/Mpeg4/MP4

开发平台:

Visual C++

  1. #define MAXNOISEXX
  2. /*
  3.  * MP3 quantization
  4.  *
  5.  * Copyright (c) 1999 Mark Taylor
  6.  *
  7.  * This program is free software; you can redistribute it and/or modify
  8.  * it under the terms of the GNU General Public License as published by
  9.  * the Free Software Foundation; either version 2, or (at your option)
  10.  * any later version.
  11.  *
  12.  * This program is distributed in the hope that it will be useful,
  13.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15.  * GNU General Public License for more details.
  16.  *
  17.  * You should have received a copy of the GNU General Public License
  18.  * along with this program; see the file COPYING.  If not, write to
  19.  * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  20.  */
  21. #include <assert.h>
  22. #include "util.h"
  23. #include "l3side.h"
  24. #include "quantize.h"
  25. #include "l3bitstream.h"
  26. #include "reservoir.h"
  27. #include "quantize-pvt.h"
  28. #ifdef HAVEGTK
  29. #include "gtkanal.h"
  30. #endif
  31. #ifdef HAVEGTK
  32. /************************************************************************/
  33. /*  updates plotting data                                               */
  34. /************************************************************************/
  35. void 
  36. set_pinfo (
  37.     gr_info *cod_info,
  38.     III_psy_ratio *ratio, 
  39.     III_scalefac_t *scalefac,
  40.     FLOAT8 xr[576],        
  41.     FLOAT8 xfsf[4][SBPSY_l],
  42.     FLOAT8 noise[4],
  43.     int gr,
  44.     int ch
  45. )
  46. {
  47.   int sfb;
  48.   FLOAT ifqstep;
  49.   int i,l,start,end,bw;
  50.   FLOAT8 en0;
  51.   D192_3 *xr_s = (D192_3 *)xr;
  52.   ifqstep = ( cod_info->scalefac_scale == 0 ) ? .5 : 1.0;
  53.   
  54.   if (cod_info->block_type == SHORT_TYPE) {
  55.     for ( i = 0; i < 3; i++ ) {
  56.       for ( sfb = 0; sfb < SBPSY_s; sfb++ )  {
  57. start = scalefac_band.s[ sfb ];
  58. end   = scalefac_band.s[ sfb + 1 ];
  59. bw = end - start;
  60. for ( en0 = 0.0, l = start; l < end; l++ ) 
  61.   en0 += (*xr_s)[l][i] * (*xr_s)[l][i];
  62. en0=Max(en0/bw,1e-20);
  63. /* conversion to FFT units */
  64. en0 = ratio->en.s[sfb][i]/en0;
  65. pinfo->xfsf_s[gr][ch][3*sfb+i] =  xfsf[i+1][sfb]*en0;
  66. pinfo->thr_s[gr][ch][3*sfb+i] = ratio->thm.s[sfb][i];
  67. pinfo->en_s[gr][ch][3*sfb+i] = ratio->en.s[sfb][i]; 
  68. pinfo->LAMEsfb_s[gr][ch][3*sfb+i]=
  69.   -2*cod_info->subblock_gain[i]-ifqstep*scalefac->s[sfb][i];
  70.       }
  71.     }
  72.   }else{
  73.     for ( sfb = 0; sfb < SBPSY_l; sfb++ )   {
  74.       start = scalefac_band.l[ sfb ];
  75.       end   = scalefac_band.l[ sfb+1 ];
  76.       bw = end - start;
  77.       for ( en0 = 0.0, l = start; l < end; l++ ) 
  78. en0 += xr[l] * xr[l];
  79.       en0=Max(en0/bw,1e-20);
  80.       /*
  81. printf("diff  = %f n",10*log10(Max(ratio[gr][ch].en.l[sfb],1e-20))
  82. -(10*log10(en0)+150));
  83.       */
  84.       
  85.       /* convert to FFT units */
  86.       en0 =   ratio->en.l[sfb]/en0;
  87.       
  88.       pinfo->xfsf[gr][ch][sfb] =  xfsf[0][sfb]*en0;
  89.       pinfo->thr[gr][ch][sfb] = ratio->thm.l[sfb];
  90.       pinfo->en[gr][ch][sfb] = ratio->en.l[sfb];
  91.       
  92.       pinfo->LAMEsfb[gr][ch][sfb]=-ifqstep*scalefac->l[sfb];
  93.       if (cod_info->preflag && sfb>=11) 
  94. pinfo->LAMEsfb[gr][ch][sfb]-=ifqstep*pretab[sfb];
  95.     }
  96.   }
  97.   pinfo->LAMEqss[gr][ch] = cod_info->global_gain;
  98.   pinfo->LAMEmainbits[gr][ch] = cod_info->part2_3_length;
  99.   pinfo->over      [gr][ch] = noise[0];
  100.   pinfo->max_noise [gr][ch] = noise[1];
  101.   pinfo->over_noise[gr][ch] = noise[2];
  102.   pinfo->tot_noise [gr][ch] = noise[3];
  103. }
  104. #endif
  105. /************************************************************************/
  106. /*  iteration_loop()                                                    */
  107. /************************************************************************/
  108. void
  109. iteration_loop( lame_global_flags *gfp,
  110.                 FLOAT8 pe[2][2], FLOAT8 ms_ener_ratio[2],
  111. FLOAT8 xr[2][2][576], III_psy_ratio ratio[2][2],
  112. III_side_info_t *l3_side, int l3_enc[2][2][576],
  113. III_scalefac_t scalefac[2][2])
  114. {
  115.   FLOAT8 xfsf[4][SBPSY_l];
  116.   FLOAT8 noise[4]; /* over,max_noise,over_noise,tot_noise; */
  117.   III_psy_xmin l3_xmin[2];
  118.   gr_info *cod_info;
  119.   int bitsPerFrame;
  120.   int mean_bits;
  121.   int ch, gr, i, bit_rate;
  122.   iteration_init(gfp,l3_side,l3_enc);
  123.   bit_rate = bitrate_table[gfp->version][gfp->bitrate_index];
  124.   getframebits(gfp,&bitsPerFrame, &mean_bits);
  125.   ResvFrameBegin(gfp, l3_side, mean_bits, bitsPerFrame );
  126.   /* quantize! */
  127.   for ( gr = 0; gr < gfp->mode_gr; gr++ ) {
  128.     int targ_bits[2];
  129.     if (convert_mdct) 
  130.       ms_convert(xr[gr], xr[gr]);
  131.     
  132.     on_pe(gfp,pe,l3_side,targ_bits,mean_bits, gr);
  133. #ifdef RH_SIDE_CBR
  134. #else
  135.     if (reduce_sidechannel) 
  136.       reduce_side(targ_bits,ms_ener_ratio[gr],mean_bits);
  137. #endif      
  138.     
  139.     for (ch=0 ; ch < gfp->stereo ; ch ++) {
  140.       cod_info = &l3_side->gr[gr].ch[ch].tt;
  141.       if (!init_outer_loop(gfp,xr[gr][ch], cod_info))
  142.         {
  143.           /* xr contains no energy 
  144.            * cod_info was set in init_outer_loop above
  145.    */
  146.           memset(&scalefac[gr][ch],0,sizeof(III_scalefac_t));
  147.           memset(l3_enc[gr][ch],0,576*sizeof(int));
  148.   noise[0]=noise[1]=noise[2]=noise[3]=0;
  149.         }
  150.       else
  151. {
  152.           calc_xmin(gfp,xr[gr][ch], &ratio[gr][ch], cod_info, &l3_xmin[ch]);
  153.   outer_loop( gfp,xr[gr][ch], targ_bits[ch], noise,
  154.       &l3_xmin[ch], l3_enc[gr][ch], 
  155.       &scalefac[gr][ch], cod_info, xfsf, ch);
  156.         }
  157.       best_scalefac_store(gfp,gr, ch, l3_enc, l3_side, scalefac);
  158.       if (gfp->use_best_huffman==1 && cod_info->block_type == NORM_TYPE) {
  159. best_huffman_divide(gr, ch, cod_info, l3_enc[gr][ch]);
  160.       }
  161. #ifdef HAVEGTK
  162.       if (gfp->gtkflag)
  163. set_pinfo (cod_info, &ratio[gr][ch], &scalefac[gr][ch], xr[gr][ch], xfsf, noise, gr, ch);
  164. #endif
  165. /*#define NORES_TEST */
  166. #ifndef NORES_TEST
  167.       ResvAdjust(gfp,cod_info, l3_side, mean_bits );
  168. #endif
  169.       /* set the sign of l3_enc */
  170.       for ( i = 0; i < 576; i++) {
  171. if (xr[gr][ch][i] < 0)
  172.   l3_enc[gr][ch][i] *= -1;
  173.       }
  174.     }
  175.   } /* loop over gr */
  176. #ifdef NORES_TEST
  177.   /* replace ResvAdjust above with this code if you do not want
  178.      the second granule to use bits saved by the first granule.
  179.      when combined with --nores, this is usefull for testing only */
  180.   for ( gr = 0; gr < gfp->mode_gr; gr++ ) {
  181.     for ( ch =  0; ch < gfp->stereo; ch++ ) {
  182. cod_info = &l3_side->gr[gr].ch[ch].tt;
  183. ResvAdjust(gfp, cod_info, l3_side, mean_bits );
  184.     }
  185.   }
  186. #endif
  187.   ResvFrameEnd(gfp,l3_side, mean_bits );
  188. }
  189. void 
  190. set_masking_lower (int VBR_q,int nbits)
  191. {
  192. FLOAT masking_lower_db, adjust;
  193. /* quality setting */
  194. /* Adjust allowed masking based on quality setting */
  195. #ifdef  RH_QUALITY_CONTROL
  196. /* - lower masking depending on Quality setting
  197.  * - quality control together with adjusted ATH MDCT scaling
  198.  *   on lower quality setting allocate more noise from
  199.  *   ATH masking, and on higher quality setting allocate
  200.  *   less noise from ATH masking.
  201.  * - experiments show that going more than 2dB over GPSYCHO's
  202.  *   limits ends up in very annoying artefacts
  203.  */
  204. static FLOAT dbQ[10]={-6.0,-4.5,-3.0,-1.5,0,0.3,0.6,1.0,1.5,2.0};
  205. assert( VBR_q <= 9 );
  206. assert( VBR_q >= 0 );
  207. masking_lower_db = dbQ[VBR_q];
  208. adjust = 0;
  209. #else
  210. /* masking_lower varies from -8 to +10 db */
  211. masking_lower_db = -6 + 2*VBR_q;
  212. /* adjust by -6(min)..0(max) depending on bitrate */
  213. adjust = (nbits-125)/(2500.0-125.0);
  214. adjust = 4*(adjust-1);
  215. #endif
  216. masking_lower_db += adjust;
  217. masking_lower = pow(10.0,masking_lower_db/10);
  218. }
  219. /************************************************************************
  220.  *
  221.  * VBR_iteration_loop()   
  222.  *
  223.  * tries to find out how many bits are needed for each granule and channel
  224.  * to get an acceptable quantization. An appropriate bitrate will then be
  225.  * choosed for quantization.  rh 8/99                                                
  226.  *
  227.  ************************************************************************/
  228. void
  229. VBR_iteration_loop (lame_global_flags *gfp,
  230.                 FLOAT8 pe[2][2], FLOAT8 ms_ener_ratio[2],
  231.                 FLOAT8 xr[2][2][576], III_psy_ratio ratio[2][2],
  232.                 III_side_info_t * l3_side, int l3_enc[2][2][576],
  233.                 III_scalefac_t scalefac[2][2])
  234. {
  235. #ifdef HAVEGTK
  236.   plotting_data bst_pinfo;
  237. #endif
  238.   gr_info         bst_cod_info, clean_cod_info;
  239.   III_scalefac_t  bst_scalefac;
  240.   int             bst_l3_enc[576]; 
  241.   
  242.   III_psy_xmin l3_xmin;
  243.   gr_info  *cod_info = NULL;
  244.   int       save_bits[2][2];
  245.   FLOAT8    noise[4];      /* over,max_noise,over_noise,tot_noise; */
  246.   FLOAT8    targ_noise[4]; /* over,max_noise,over_noise,tot_noise; */
  247.   FLOAT8    xfsf[4][SBPSY_l];
  248.   int       this_bits, dbits;
  249.   int       used_bits=0;
  250.   int       min_bits,max_bits,min_mean_bits=0;
  251.   int       frameBits[15];
  252.   int       bitsPerFrame;
  253.   int       bits;
  254.   int       mean_bits;
  255.   int       i,ch, gr, analog_silence;
  256.   int     reparted = 0;
  257.   iteration_init(gfp,l3_side,l3_enc);
  258. #ifdef RH_QUALITY_CONTROL
  259.   /* with RH_QUALITY_CONTROL we have to set masking_lower only once */
  260.   set_masking_lower(gfp->VBR_q, 0 );
  261. #endif      
  262.   /*******************************************************************
  263.    * how many bits are available for each bitrate?
  264.    *******************************************************************/
  265.   for( gfp->bitrate_index = 1;
  266.        gfp->bitrate_index <= gfp->VBR_max_bitrate;
  267.        gfp->bitrate_index++    ) {
  268.     getframebits (gfp,&bitsPerFrame, &mean_bits);
  269.     if (gfp->bitrate_index == gfp->VBR_min_bitrate) {
  270.       /* always use at least this many bits per granule per channel */
  271.       /* unless we detect analog silence, see below */
  272.       min_mean_bits=mean_bits/gfp->stereo;
  273.     }
  274.     frameBits[gfp->bitrate_index]=
  275.       ResvFrameBegin (gfp,l3_side, mean_bits, bitsPerFrame);
  276.   }
  277.   gfp->bitrate_index=gfp->VBR_max_bitrate;
  278.   
  279.   /*******************************************************************
  280.    * how many bits would we use of it?
  281.    *******************************************************************/
  282.   analog_silence=0;
  283.   for (gr = 0; gr < gfp->mode_gr; gr++) {
  284.     int num_chan=gfp->stereo;
  285. #ifdef  RH_SIDE_VBR
  286.     /* my experiences are, that side channel reduction  
  287.      * does more harm than good when VBR encoding
  288.      * (Robert.Hegemann@gmx.de 2000-02-18)
  289.      */
  290. #else
  291.     /* determine quality based on mid channel only */
  292.     if (reduce_sidechannel) num_chan=1;  
  293. #endif
  294.     /* copy data to be quantized into xr */
  295.     if (convert_mdct)
  296. ms_convert(xr[gr],xr[gr]);
  297.     for (ch = 0; ch < num_chan; ch++) { 
  298.       int real_bits;
  299.       
  300.       /******************************************************************
  301.        * find smallest number of bits for an allowable quantization
  302.        ******************************************************************/
  303.       cod_info = &l3_side->gr[gr].ch[ch].tt;
  304.       min_bits = Max(125,min_mean_bits);
  305.       if (!init_outer_loop(gfp,xr[gr][ch], cod_info))
  306.       {
  307.         /* xr contains no energy 
  308.          * cod_info was set in init_outer_loop above
  309.  */
  310.         memset(&scalefac[gr][ch],0,sizeof(III_scalefac_t));
  311.         memset(l3_enc[gr][ch],0,576*sizeof(int));
  312.         save_bits[gr][ch] = 0;
  313. #ifdef HAVEGTK
  314. if (gfp->gtkflag)
  315.   set_pinfo(cod_info, &ratio[gr][ch], &scalefac[gr][ch], xr[gr][ch], xfsf, noise, gr, ch);
  316. #endif
  317. analog_silence=1;
  318. continue; /* with next channel */
  319.       }
  320.       
  321.       memcpy( &clean_cod_info, cod_info, sizeof(gr_info) );
  322.       
  323. #ifdef RH_QUALITY_CONTROL
  324.       /*
  325.        * masking lower already set in the beginning
  326.        */
  327. #else
  328.       /*
  329.        * has to be set before calculating l3_xmin
  330.        */
  331.       set_masking_lower( gfp->VBR_q,2500 );
  332. #endif      
  333.       /* check for analolg silence */
  334.       /* if energy < ATH, set min_bits = 125 */
  335.       if (0==calc_xmin(gfp,xr[gr][ch], &ratio[gr][ch], cod_info, &l3_xmin)) {
  336.   analog_silence=1;
  337.   min_bits=125;
  338.       }
  339.       if (cod_info->block_type==SHORT_TYPE) {
  340.   min_bits += Max(1100,pe[gr][ch]);
  341.   min_bits=Min(min_bits,1800);
  342.       }
  343.       max_bits = 1200 + frameBits[gfp->VBR_max_bitrate]/(gfp->stereo*gfp->mode_gr);
  344.       max_bits=Min(max_bits,2500);
  345.       max_bits=Max(max_bits,min_bits);
  346.       dbits = (max_bits-min_bits)/4;
  347.       this_bits = (max_bits+min_bits)/2;
  348.       real_bits = max_bits+1;
  349.       /* bin search to within +/- 10 bits of optimal */
  350.       do {
  351.   int better;
  352.   assert(this_bits>=min_bits);
  353.   assert(this_bits<=max_bits);
  354.   if( this_bits >= real_bits ){
  355.       /* 
  356.        * we already found a quantization with fewer bits
  357.        * so we can skip this try
  358.        */
  359.       this_bits -= dbits;
  360.       dbits /= 2;
  361.       continue; /* skips the rest of this do-while loop */
  362.   }
  363.   /* VBR will look for a quantization which has better values
  364.    * then those specified below.*/
  365.   targ_noise[0]=0;          /* over */
  366.   targ_noise[1]=0;          /* max_noise */
  367.   targ_noise[2]=0;          /* over_noise */
  368.   targ_noise[3]=0;          /* tot_noise */
  369.   targ_noise[0]=Max(0,targ_noise[0]);
  370.   targ_noise[2]=Max(0,targ_noise[2]);
  371.   /*
  372.    *  OK, start with a fresh setting
  373.    *  - scalefac  will be set up by outer_loop
  374.    *  - l3_enc    will be set up by outer_loop
  375.    *  + cod_info  we will restore our initialized one, see below
  376.    */
  377.   memcpy( cod_info, &clean_cod_info, sizeof(gr_info) );
  378. #ifdef RH_QUALITY_CONTROL
  379.           /*
  380.    * there is no need for a recalculation of l3_xmin,
  381.    * because masking_lower did not change within this do-while
  382.    */
  383. #else
  384.   /* quality setting */
  385.   set_masking_lower( gfp->VBR_q,this_bits );
  386.           /* 
  387.    * compute max allowed distortion, masking lower has changed
  388.    */
  389.           calc_xmin(gfp,xr[gr][ch], &ratio[gr][ch], cod_info, &l3_xmin);
  390. #endif
  391.   outer_loop( gfp,xr[gr][ch], this_bits, noise, 
  392.       &l3_xmin, l3_enc[gr][ch],
  393.       &scalefac[gr][ch], cod_info, xfsf,
  394.       ch);
  395.   /* is quantization as good as we are looking for ? */
  396.   better=VBR_compare((int)targ_noise[0],targ_noise[3],targ_noise[2],
  397.      targ_noise[1],(int)noise[0],noise[3],noise[2],
  398.      noise[1]);
  399. #ifdef HAVEGTK
  400.   if (gfp->gtkflag)
  401.     set_pinfo(cod_info, &ratio[gr][ch], &scalefac[gr][ch], xr[gr][ch], xfsf, noise, gr, ch);
  402. #endif
  403.   if (better) {
  404.       /* 
  405.        * we now know it can be done with "real_bits"
  406.        * and maybe we can skip some iterations
  407.        */
  408.       real_bits = cod_info->part2_3_length;
  409.       /*
  410.        * save best quantization so far
  411.        */
  412.               memcpy( &bst_scalefac, &scalefac[gr][ch], sizeof(III_scalefac_t)  );
  413.               memcpy(  bst_l3_enc,    l3_enc  [gr][ch], sizeof(int)*576         );
  414.               memcpy( &bst_cod_info,  cod_info,         sizeof(gr_info)         );
  415. #ifdef HAVEGTK
  416.               if (gfp->gtkflag) 
  417.                 memcpy( &bst_pinfo, pinfo, sizeof(plotting_data) );
  418. #endif
  419.       /*
  420.        * try with fewer bits
  421.        */
  422.       this_bits -= dbits;
  423.   } else {
  424.       /*
  425.        * try with more bits
  426.        */
  427.       this_bits += dbits;
  428.   }
  429.   dbits /= 2;
  430.       } while (dbits>10) ;
  431.       
  432.       if (real_bits <= max_bits)
  433.       {
  434.         /* restore best quantization found */
  435.         memcpy(  cod_info,         &bst_cod_info, sizeof(gr_info)        );
  436.         memcpy( &scalefac[gr][ch], &bst_scalefac, sizeof(III_scalefac_t) );
  437.         memcpy(  l3_enc  [gr][ch],  bst_l3_enc,   sizeof(int)*576        );
  438. #ifdef HAVEGTK
  439.         if (gfp->gtkflag) 
  440.           memcpy( pinfo, &bst_pinfo, sizeof(plotting_data) );
  441. #endif
  442.       }
  443.       assert((int)cod_info->part2_3_length <= max_bits);
  444.       save_bits[gr][ch] = cod_info->part2_3_length;
  445.       used_bits += save_bits[gr][ch];
  446.       
  447.     } /* for ch */
  448.   } /* for gr */
  449. #ifdef  RH_SIDE_VBR
  450.   /* my experiences are, that side channel reduction  
  451.    * does more harm than good when VBR encoding
  452.    * (Robert.Hegemann@gmx.de 2000-02-18)
  453.    */
  454. #else
  455.   if (reduce_sidechannel) {
  456.     /* number of bits needed was found for MID channel above.  Use formula
  457.      * (fixed bitrate code) to set the side channel bits */
  458.     for (gr = 0; gr < gfp->mode_gr; gr++) {
  459.       FLOAT8 fac = .33*(.5-ms_ener_ratio[gr])/.5;
  460.       save_bits[gr][1]=((1-fac)/(1+fac))*save_bits[gr][0];
  461.       save_bits[gr][1]=Max(125,save_bits[gr][1]);
  462.       used_bits += save_bits[gr][1];
  463.     }
  464.   }
  465. #endif
  466.   /******************************************************************
  467.    * find lowest bitrate able to hold used bits
  468.    ******************************************************************/
  469.   for( gfp->bitrate_index =   (analog_silence ? 1 : gfp->VBR_min_bitrate );
  470.        gfp->bitrate_index < gfp->VBR_max_bitrate;
  471.        gfp->bitrate_index++    )
  472.     if( used_bits <= frameBits[gfp->bitrate_index] ) break;
  473.   /*******************************************************************
  474.    * calculate quantization for this bitrate
  475.    *******************************************************************/  
  476.   getframebits (gfp,&bitsPerFrame, &mean_bits);
  477.   bits=ResvFrameBegin (gfp,l3_side, mean_bits, bitsPerFrame);
  478.   /* repartion available bits in same proportion */
  479.   if (used_bits > bits ) {
  480.     reparted = 1;
  481.     for( gr = 0; gr < gfp->mode_gr; gr++) {
  482.       for(ch = 0; ch < gfp->stereo; ch++) {
  483. save_bits[gr][ch]=(save_bits[gr][ch]*frameBits[gfp->bitrate_index])/used_bits;
  484.       }
  485.     }
  486.     used_bits=0;
  487.     for( gr = 0; gr < gfp->mode_gr; gr++) {
  488.       for(ch = 0; ch < gfp->stereo; ch++) {
  489. used_bits += save_bits[gr][ch];
  490.       }
  491.     }
  492.   }
  493.   assert(used_bits <= bits);
  494.   for(gr = 0; gr < gfp->mode_gr; gr++) {
  495.     for(ch = 0; ch < gfp->stereo; ch++) {
  496. #ifdef RH_SIDE_VBR
  497.       if (reparted)
  498. #else
  499.       if (reparted || (reduce_sidechannel && ch == 1))
  500. #endif
  501.       {
  502.         cod_info = &l3_side->gr[gr].ch[ch].tt;
  503.        
  504. if (!init_outer_loop(gfp,xr[gr][ch], cod_info))
  505.         {
  506.           /* xr contains no energy 
  507.            * cod_info was set in init_outer_loop above
  508.    */
  509.           memset(&scalefac[gr][ch],0,sizeof(III_scalefac_t));
  510.           memset(l3_enc[gr][ch],0,576*sizeof(int));
  511.   noise[0]=noise[1]=noise[2]=noise[3]=0;
  512.         }
  513. else
  514. {
  515. #ifdef RH_QUALITY_CONTROL
  516.           /*
  517.            * masking lower already set in the beginning
  518.            */
  519. #else
  520.           /* quality setting */
  521.           set_masking_lower( gfp->VBR_q,save_bits[gr][ch] );
  522. #endif
  523.           calc_xmin(gfp,xr[gr][ch], &ratio[gr][ch], cod_info, &l3_xmin);
  524.           outer_loop( gfp,xr[gr][ch], save_bits[gr][ch], noise,
  525.         &l3_xmin, l3_enc[gr][ch], 
  526.       &scalefac[gr][ch], cod_info, xfsf, ch);
  527. }
  528. #ifdef HAVEGTK
  529. if (gfp->gtkflag)
  530.   set_pinfo(cod_info, &ratio[gr][ch], &scalefac[gr][ch], xr[gr][ch], xfsf, noise, gr, ch);
  531. #endif
  532.       }
  533.     }
  534.   }
  535.   /*******************************************************************
  536.    * update reservoir status after FINAL quantization/bitrate 
  537.    *******************************************************************/
  538.   for (gr = 0; gr < gfp->mode_gr; gr++)
  539.     for (ch = 0; ch < gfp->stereo; ch++) {
  540.       cod_info = &l3_side->gr[gr].ch[ch].tt;
  541.       best_scalefac_store(gfp,gr, ch, l3_enc, l3_side, scalefac);
  542.       if (cod_info->block_type == NORM_TYPE) {
  543. best_huffman_divide(gr, ch, cod_info, l3_enc[gr][ch]);
  544.       }
  545. #ifdef HAVEGTK
  546.       if (gfp->gtkflag)
  547. pinfo->LAMEmainbits[gr][ch]=cod_info->part2_3_length;
  548. #endif
  549.       ResvAdjust (gfp,cod_info, l3_side, mean_bits);
  550.     }
  551.   /*******************************************************************
  552.    * set the sign of l3_enc 
  553.    *******************************************************************/
  554.   for (gr = 0; gr < gfp->mode_gr; gr++)
  555.     for (ch = 0; ch < gfp->stereo; ch++) {
  556. /*
  557.  * is the following code correct?
  558.  *
  559.       int      *pi = &l3_enc[gr][ch][0];
  560.       for (i = 0; i < 576; i++) {
  561.         FLOAT8    pr = xr[gr][ch][i];
  562.         if ((pr < 0) && (pi[i] > 0))
  563.           pi[i] *= -1;
  564.       }
  565.  *
  566.  * or is the code used for CBR correct?
  567.  */
  568.       for ( i = 0; i < 576; i++) {
  569.         if (xr[gr][ch][i] < 0) l3_enc[gr][ch][i] *= -1;
  570.       }
  571.     }
  572.   ResvFrameEnd (gfp,l3_side, mean_bits);
  573. }
  574. /************************************************************************/
  575. /*  init_outer_loop  mt 6/99                                            */
  576. /*  returns 0 if all energies in xr are zero, else 1                    */
  577. /************************************************************************/
  578. int init_outer_loop(lame_global_flags *gfp,
  579.     FLOAT8 xr[576],        /*  could be L/R OR MID/SIDE */
  580.     gr_info *cod_info)
  581. {
  582.   int i;
  583.   for ( i = 0; i < 4; i++ )
  584.     cod_info->slen[i] = 0;
  585.   cod_info->sfb_partition_table = &nr_of_sfb_block[0][0][0];
  586.   cod_info->part2_3_length    = 0;
  587.   cod_info->big_values        = 0;
  588.   cod_info->count1            = 0;
  589.   cod_info->scalefac_compress = 0;
  590.   cod_info->table_select[0]   = 0;
  591.   cod_info->table_select[1]   = 0;
  592.   cod_info->table_select[2]   = 0;
  593.   cod_info->subblock_gain[0]  = 0;
  594.   cod_info->subblock_gain[1]  = 0;
  595.   cod_info->subblock_gain[2]  = 0;
  596.   cod_info->region0_count     = 0;
  597.   cod_info->region1_count     = 0;
  598.   cod_info->part2_length      = 0;
  599.   cod_info->preflag           = 0;
  600.   cod_info->scalefac_scale    = 0;
  601.   cod_info->global_gain       = 210;
  602.   cod_info->count1table_select= 0;
  603.   cod_info->count1bits        = 0;
  604.   
  605.   
  606.   if (gfp->experimentalZ) {
  607.     /* compute subblock gains */
  608.     int j,b;  FLOAT8 en[3],mx;
  609.     if ((cod_info->block_type==SHORT_TYPE) ) {
  610.       /* estimate energy within each subblock */
  611.       for (b=0; b<3; b++) en[b]=0;
  612.       for ( i=0,j = 0; j < 192; j++ ) {
  613. for (b=0; b<3; b++) {
  614.   en[b]+=xr[i] * xr[i];
  615.   i++;
  616. }
  617.       }
  618.       mx = 1e-12;
  619.       for (b=0; b<3; b++) mx=Max(mx,en[b]);
  620.       for (b=0; b<3; b++) en[b] = Max(en[b],1e-12)/mx;
  621.       /*printf("ener = %4.2f  %4.2f  %4.2f  n",en[0],en[1],en[2]);*/
  622.       /* pick gain so that 2^(2gain)*en[0] = 1  */
  623.       /* gain = .5* log( 1/en[0] )/LOG2 = -.5*log(en[])/LOG2 */
  624.       for (b=0; b<3; b++) {
  625. cod_info->subblock_gain[b] = (int)(-.5*log(en[b])/LOG2 + 0.5);
  626. if (cod_info->subblock_gain[b] > 2) 
  627.   cod_info->subblock_gain[b]=2;
  628. if (cod_info->subblock_gain[b] < 0) 
  629.   cod_info->subblock_gain[b]=0;
  630.       }
  631.       /*
  632.        *  check if there is some energy we have to quantize
  633.        *  if so, then return 1 else 0
  634.        */
  635.       if (1e-99 < en[0]+en[1]+en[2])
  636.         return 1;
  637.       else
  638.         return 0;
  639.     }
  640.   }
  641.   /*
  642.    *  check if there is some energy we have to quantize
  643.    *  if so, then return 1 else 0
  644.    */
  645.   for (i=0; i<576; i++) 
  646.     if ( 1e-99 < fabs (xr[i]) )
  647.       return 1;
  648.   
  649.   return 0;
  650. }
  651. /************************************************************************/
  652. /*  outer_loop                                                         */
  653. /************************************************************************/
  654. /*  Function: The outer iteration loop controls the masking conditions  */
  655. /*  of all scalefactorbands. It computes the best scalefac and          */
  656. /*  global gain. This module calls the inner iteration loop             
  657.  * 
  658.  *  mt 5/99 completely rewritten to allow for bit reservoir control,   
  659.  *  mid/side channels with L/R or mid/side masking thresholds, 
  660.  *  and chooses best quantization instead of last quantization when 
  661.  *  no distortion free quantization can be found.  
  662.  *  
  663.  *  added VBR support mt 5/99
  664.  ************************************************************************/
  665. void outer_loop(
  666.     lame_global_flags *gfp,
  667.     FLOAT8 xr[576],        
  668.     int targ_bits,
  669.     FLOAT8 best_noise[4],
  670.     III_psy_xmin *l3_xmin,   /* the allowed distortion of the scalefactor */
  671.     int l3_enc[576],         /* vector of quantized values ix(0..575) */
  672.     III_scalefac_t *scalefac, /* scalefactors */
  673.     gr_info *cod_info,
  674.     FLOAT8 xfsf[4][SBPSY_l],
  675.     int ch)
  676. {
  677.   III_scalefac_t scalefac_w;
  678.   gr_info save_cod_info;
  679.   int l3_enc_w[576]; 
  680.   int i, iteration;
  681.   int status,bits_found=0;
  682.   int huff_bits;
  683.   FLOAT8 xrpow[576],temp;
  684.   int better;
  685.   int over=0;
  686.   FLOAT8 max_noise;
  687.   FLOAT8 over_noise;
  688.   FLOAT8 tot_noise;
  689.   int best_over=100;
  690.   FLOAT8 best_max_noise=0;
  691.   FLOAT8 best_over_noise=0;
  692.   FLOAT8 best_tot_noise=0;
  693.   FLOAT8 xfsf_w[4][SBPSY_l];
  694.   FLOAT8 distort[4][SBPSY_l];
  695.   int compute_stepsize=1;
  696.   int notdone=1;
  697.   /* BEGIN MAIN LOOP */
  698.   iteration = 0;
  699.   while ( notdone  ) {
  700.     static int OldValue[2] = {180, 180};
  701.     int try_scale=0;
  702.     iteration ++;
  703.     if (compute_stepsize) {
  704.       /* init and compute initial quantization step */
  705.       compute_stepsize=0;
  706.       /* reset of iteration variables */
  707.       memset(&scalefac_w, 0, sizeof(III_scalefac_t));
  708.       for (i=0;i<576;i++) {
  709. temp=fabs(xr[i]);
  710. xrpow[i]=sqrt(sqrt(temp)*temp);
  711.       }
  712.       bits_found=bin_search_StepSize2(gfp,targ_bits,OldValue[ch],
  713.       l3_enc_w,xrpow,cod_info);
  714.       OldValue[ch] = cod_info->global_gain;
  715.     }
  716.     /* inner_loop starts with the initial quantization step computed above
  717.      * and slowly increases until the bits < huff_bits.
  718.      * Thus it is important not to start with too large of an inital
  719.      * quantization step.  Too small is ok, but inner_loop will take longer 
  720.      */
  721.     huff_bits = targ_bits - cod_info->part2_length;
  722.     if (huff_bits < 0) {
  723.       assert(iteration != 1);
  724.       /* scale factors too large, not enough bits. use previous quantizaton */
  725.       notdone=0;
  726.     } else {
  727.       /* if this is the first iteration, see if we can reuse the quantization
  728.        * computed in bin_search_StepSize above */
  729.       int real_bits;
  730.       if (iteration==1) {
  731. if(bits_found>huff_bits) {
  732.   cod_info->global_gain++;
  733.   real_bits = inner_loop(gfp,xrpow, l3_enc_w, huff_bits, cod_info);
  734. } else real_bits=bits_found;
  735.       }
  736.       else 
  737. real_bits=inner_loop(gfp,xrpow, l3_enc_w, huff_bits, cod_info);
  738.       cod_info->part2_3_length = real_bits;
  739.       /* compute the distortion in this quantization */
  740.       if (gfp->noise_shaping==0) {
  741.        over=0;
  742.       }else{
  743. /* coefficients and thresholds both l/r (or both mid/side) */
  744. over=calc_noise1( xr, l3_enc_w, cod_info, 
  745.   xfsf_w,distort, l3_xmin, &scalefac_w, &over_noise, 
  746.   &tot_noise, &max_noise);
  747.       }
  748.       /* check if this quantization is better the our saved quantization */
  749.       if (iteration == 1) better=1;
  750.       else 
  751. better=quant_compare(gfp->experimentalX,
  752.      best_over,best_tot_noise,best_over_noise,best_max_noise,
  753.                   over,     tot_noise,     over_noise,     max_noise);
  754.       /* save data so we can restore this quantization later */    
  755.       if (better) {
  756. best_over=over;
  757. best_max_noise=max_noise;
  758. best_over_noise=over_noise;
  759. best_tot_noise=tot_noise;
  760. memcpy(scalefac, &scalefac_w, sizeof(III_scalefac_t));
  761. memcpy(l3_enc,l3_enc_w,sizeof(int)*576);
  762. memcpy(&save_cod_info,cod_info,sizeof(save_cod_info));
  763. #ifdef HAVEGTK
  764. if (gfp->gtkflag) {
  765.   memcpy(xfsf, xfsf_w, sizeof(xfsf_w));
  766. }
  767. #endif
  768.       }
  769.     }
  770.     
  771.     /* if no bands with distortion, we are done */
  772.     if (gfp->noise_shaping_stop==0)
  773.       if (over==0) notdone=0;
  774.     if (notdone) {
  775. amp_scalefac_bands( xrpow, cod_info, &scalefac_w, distort);
  776. /* check to make sure we have not amplified too much */
  777. /* loop_break returns 0 if there is an unamplified scalefac */
  778. /* scale_bitcount returns 0 if no scalefactors are too large */
  779. if ( (status = loop_break(&scalefac_w, cod_info)) == 0 ) {
  780.     if ( gfp->version == 1 ) {
  781. status = scale_bitcount(&scalefac_w, cod_info);
  782.     }else{
  783. status = scale_bitcount_lsf(&scalefac_w, cod_info);
  784.     }
  785.     if (status && (cod_info->scalefac_scale==0)) try_scale=1; 
  786. }
  787. notdone = !status;
  788.     }
  789.     if (try_scale && gfp->experimentalY) {
  790.       init_outer_loop(gfp,xr, cod_info);
  791.       compute_stepsize=1;  /* compute a new global gain */
  792.       notdone=1;
  793.       cod_info->scalefac_scale=1;
  794.     }
  795.   }    /* done with main iteration */
  796.   memcpy(cod_info,&save_cod_info,sizeof(save_cod_info));
  797.   cod_info->part2_3_length += cod_info->part2_length;
  798.   /* finish up */
  799.   assert( cod_info->global_gain < 256 );
  800.   best_noise[0]=best_over;
  801.   best_noise[1]=best_max_noise;
  802.   best_noise[2]=best_over_noise;
  803.   best_noise[3]=best_tot_noise;
  804. }
  805.   
  806. /*************************************************************************/
  807. /*            calc_noise                                                 */
  808. /*************************************************************************/
  809. /*  mt 5/99:  Function: Improved calc_noise for a single channel   */
  810. int calc_noise1( FLOAT8 xr[576], int ix[576], gr_info *cod_info,
  811.  FLOAT8 xfsf[4][SBPSY_l], FLOAT8 distort[4][SBPSY_l],
  812.  III_psy_xmin *l3_xmin, III_scalefac_t *scalefac,
  813.  FLOAT8 *over_noise,
  814.  FLOAT8 *tot_noise, FLOAT8 *max_noise)
  815. {
  816.     int start, end, l, i, over=0;
  817. u_int sfb;
  818.     FLOAT8 sum,step,bw;
  819. #ifdef RH_ATH
  820.     FLOAT8 ath_max;
  821. #endif
  822.     int count=0;
  823.     FLOAT8 noise;
  824.     *over_noise=0;
  825.     *tot_noise=0;
  826.     *max_noise = -999;
  827.     for ( sfb = 0; sfb < cod_info->sfb_lmax; sfb++ ) {
  828. FLOAT8 step;
  829. int s = scalefac->l[sfb];
  830. if (cod_info->preflag)
  831.     s += pretab[sfb];
  832. s = cod_info->global_gain - (s << (cod_info->scalefac_scale + 1));
  833. assert(s<Q_MAX);
  834. assert(s>=0);
  835. step = POW20(s);
  836. start = scalefac_band.l[ sfb ];
  837.         end   = scalefac_band.l[ sfb+1 ];
  838.         bw = end - start;
  839. #ifdef RH_ATH
  840.         ath_max = 0;
  841. #endif
  842.         for ( sum = 0.0, l = start; l < end; l++ )
  843.         {
  844.             FLOAT8 temp;
  845.             temp = fabs(xr[l]) - pow43[ix[l]] * step;
  846. #ifdef MAXNOISE
  847.     temp = bw*temp*temp;
  848.     sum = Max(sum,temp);
  849. #elif RH_ATH
  850.     temp = temp*temp;
  851.             sum += temp;
  852.     ath_max = Max( ath_max, temp/ATH_mdct_long[l] );
  853. #else
  854.             sum += temp * temp;
  855. #endif
  856.     
  857.         }
  858.         xfsf[0][sfb] = sum / bw;
  859. /* max -30db noise below threshold */
  860. #ifdef RH_ATH
  861. noise = 10*log10(Max(.001,Min(ath_max,xfsf[0][sfb]/l3_xmin->l[sfb])));
  862. #else
  863. noise = 10*log10(Max(.001,xfsf[0][sfb] / l3_xmin->l[sfb]));
  864. #endif
  865.         distort[0][sfb] = noise;
  866.         if (noise>0) {
  867.   over++;
  868.   *over_noise += noise;
  869. }
  870. *tot_noise += noise;
  871. *max_noise=Max(*max_noise,noise);
  872. count++;
  873.     }
  874.     for ( i = 0; i < 3; i++ ) {
  875.         for ( sfb = cod_info->sfb_smax; sfb < SBPSY_s; sfb++ ) {
  876.     int s;
  877.     s = (scalefac->s[sfb][i] << (cod_info->scalefac_scale + 1))
  878. + cod_info->subblock_gain[i] * 8;
  879.     s = cod_info->global_gain - s;
  880.     assert(s<Q_MAX);
  881.     assert(s>=0);
  882.     step = POW20(s);
  883.     start = scalefac_band.s[ sfb ];
  884.     end   = scalefac_band.s[ sfb+1 ];
  885.             bw = end - start;
  886. #ifdef RH_ATH
  887.     ath_max = 0;
  888. #endif
  889.     for ( sum = 0.0, l = start; l < end; l++ ) {
  890. FLOAT8 temp;
  891. temp = fabs(xr[l * 3 + i]) - pow43[ix[l * 3 + i]] * step;
  892. #ifdef MAXNOISE
  893. temp = bw*temp*temp;
  894. sum = Max(sum,temp);
  895. #elif RH_ATH
  896. temp = temp*temp;
  897. sum += temp;
  898. ath_max = Max( ath_max, temp/ATH_mdct_short[l] );
  899. #else
  900. sum += temp * temp;
  901. #endif
  902.             }       
  903.     xfsf[i+1][sfb] = sum / bw;
  904.     /* max -30db noise below threshold */
  905. #ifdef RH_ATH
  906.     noise = 10*log10(Max(.001,Min(ath_max,xfsf[i+1][sfb]/l3_xmin->s[sfb][i])));
  907. #else
  908.     noise = 10*log10(Max(.001,xfsf[i+1][sfb] / l3_xmin->s[sfb][i] ));
  909. #endif
  910.             distort[i+1][sfb] = noise;
  911.             if (noise > 0) {
  912. over++;
  913. *over_noise += noise;
  914.     }
  915.     *tot_noise += noise;
  916.     *max_noise=Max(*max_noise,noise);
  917.     count++;     
  918.         }
  919.     }
  920.     if (count>1) *tot_noise /= count;
  921.     if (over>1) *over_noise /= over;
  922.     return over;
  923. }
  924. /*************************************************************************/
  925. /*            amp_scalefac_bands                                         */
  926. /*************************************************************************/
  927. /* 
  928.   Amplify the scalefactor bands that violate the masking threshold.
  929.   See ISO 11172-3 Section C.1.5.4.3.5
  930. */
  931. void amp_scalefac_bands(FLOAT8 xrpow[576], 
  932. gr_info *cod_info,
  933. III_scalefac_t *scalefac,
  934. FLOAT8 distort[4][SBPSY_l])
  935. {
  936.     int start, end, l,i;
  937. u_int sfb;
  938.     FLOAT8 ifqstep34;
  939.     FLOAT8 distort_thresh;
  940.     if ( cod_info->scalefac_scale == 0 )
  941. ifqstep34 = 1.29683955465100964055;
  942.     else
  943. ifqstep34 = 1.68179283050742922612;
  944.     /* distort_thresh = 0, unless all bands have distortion 
  945.      * less than masking.  In that case, just amplify bands with distortion
  946.      * within 95% of largest distortion/masking ratio */
  947.     distort_thresh = -900;
  948.     for ( sfb = 0; sfb < cod_info->sfb_lmax; sfb++ ) {
  949. distort_thresh = Max(distort[0][sfb],distort_thresh);
  950.     }
  951.     for ( sfb = cod_info->sfb_smax; sfb < 12; sfb++ ) {
  952. for ( i = 0; i < 3; i++ ) {
  953.     distort_thresh = Max(distort[i+1][sfb],distort_thresh);
  954. }
  955.     }
  956.     distort_thresh=Min(distort_thresh * 1.05, 0.0);
  957.     for ( sfb = 0; sfb < cod_info->sfb_lmax; sfb++ ) {
  958. if ( distort[0][sfb]>distort_thresh  ) {
  959.     scalefac->l[sfb]++;
  960.     start = scalefac_band.l[sfb];
  961.     end   = scalefac_band.l[sfb+1];
  962.     for ( l = start; l < end; l++ )
  963. xrpow[l] *= ifqstep34;
  964. }
  965.     }
  966.     for ( i = 0; i < 3; i++ ) {
  967. for ( sfb = cod_info->sfb_smax; sfb < 12; sfb++ ) {
  968.             if ( distort[i+1][sfb]>distort_thresh) {
  969.                 scalefac->s[sfb][i]++;
  970.                 start = scalefac_band.s[sfb];
  971.                 end   = scalefac_band.s[sfb+1];
  972. for (l = start; l < end; l++)
  973.     xrpow[l * 3 + i] *= ifqstep34;
  974.             }
  975. }
  976.     }
  977. }
  978. int quant_compare(int experimentalX,
  979. int best_over,FLOAT8 best_tot_noise,FLOAT8 best_over_noise,FLOAT8 best_max_noise,
  980. int over,FLOAT8 tot_noise, FLOAT8 over_noise, FLOAT8 max_noise)
  981. {
  982.   /*
  983.     noise is given in decibals (db) relative to masking thesholds.
  984.     over_noise:  sum of quantization noise > masking
  985.     tot_noise:   sum of all quantization noise
  986.     max_noise:   max quantization noise 
  987.    */
  988.   int better=0;
  989.   if (experimentalX==0) {
  990.     better = ((over < best_over) ||
  991.       ((over==best_over) && (over_noise<=best_over_noise)) ) ;
  992.   }
  993.   if (experimentalX==1) 
  994.     better = max_noise < best_max_noise;
  995.   if (experimentalX==2) {
  996.     better = tot_noise < best_tot_noise;
  997.   }
  998.   if (experimentalX==3) {
  999.     better = (tot_noise < best_tot_noise) &&
  1000.       (max_noise < best_max_noise + 2);
  1001.   }
  1002.   if (experimentalX==4) {
  1003.     better = ( ( (0>=max_noise) && (best_max_noise>2)) ||
  1004.      ( (0>=max_noise) && (best_max_noise<0) && ((best_max_noise+2)>max_noise) && (tot_noise<best_tot_noise) ) ||
  1005.      ( (0>=max_noise) && (best_max_noise>0) && ((best_max_noise+2)>max_noise) && (tot_noise<(best_tot_noise+best_over_noise)) ) ||
  1006.      ( (0<max_noise) && (best_max_noise>-0.5) && ((best_max_noise+1)>max_noise) && ((tot_noise+over_noise)<(best_tot_noise+best_over_noise)) ) ||
  1007.      ( (0<max_noise) && (best_max_noise>-1) && ((best_max_noise+1.5)>max_noise) && ((tot_noise+over_noise+over_noise)<(best_tot_noise+best_over_noise+best_over_noise)) ) );
  1008.   }
  1009.   if (experimentalX==5) {
  1010.     better =   (over_noise <  best_over_noise)
  1011.       || ((over_noise == best_over_noise)&&(tot_noise < best_tot_noise));
  1012.   }
  1013.   if (experimentalX==6) {
  1014.     better = (over_noise < best_over_noise)
  1015.            ||( (over_noise == best_over_noise)
  1016.              &&( (max_noise < best_max_noise)
  1017.                ||( (max_noise == best_max_noise)
  1018.                  &&(tot_noise <= best_tot_noise)
  1019.                  )
  1020.                ) 
  1021.      );
  1022.   }
  1023.   return better;
  1024. }
  1025. int VBR_compare(
  1026. int best_over,FLOAT8 best_tot_noise,FLOAT8 best_over_noise,FLOAT8 best_max_noise,
  1027. int over,FLOAT8 tot_noise, FLOAT8 over_noise, FLOAT8 max_noise)
  1028. {
  1029.   /*
  1030.     noise is given in decibals (db) relative to masking thesholds.
  1031.     over_noise:  sum of quantization noise > masking
  1032.     tot_noise:   sum of all quantization noise
  1033.     max_noise:   max quantization noise 
  1034.    */
  1035.   int better=0;
  1036.   better = ((over <= best_over) &&
  1037.     (over_noise<=best_over_noise) &&
  1038.     (tot_noise<=best_tot_noise) &&
  1039.     (max_noise<=best_max_noise));
  1040.   return better;
  1041. }
  1042.