/* *	MP3 quantization * *	Copyright (c) 1999 Mark Taylor * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2, or (at your option) * any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; see the file COPYING.  If not, write to * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */#include <assert.h>#include "util.h"#include "l3side.h"#include "quantize.h"#include "reservoir.h"#include "quantize-pvt.h"#include "gtkanal.h"#if (defined(__GNUC__) && defined(__i386__))#define USE_GNUC_ASM#endif#ifdef _MSC_VER#define USE_MSC_ASM#endif/********************************************************************* * XRPOW_FTOI is a macro to convert floats to ints.   * if XRPOW_FTOI(x) = nearest_int(x), then QUANTFAC(x)=adj43asm[x] *                                         ROUNDFAC= -0.0946 * * if XRPOW_FTOI(x) = floor(x), then QUANTFAC(x)=asj43[x]    *                                   ROUNDFAC=0.4054 *********************************************************************/#ifdef USE_GNUC_ASM#  define QUANTFAC(rx)  adj43asm[rx]#  define ROUNDFAC -0.0946#  define XRPOW_FTOI(src, dest) \     asm ("fistpl %0 " : "=m"(dest) : "t"(src) : "st")#elif defined (USE_MSC_ASM)#  define QUANTFAC(rx)  adj43asm[rx]#  define ROUNDFAC -0.0946#  define XRPOW_FTOI(src, dest) do { \     FLOAT8 src_ = (src); \     int dest_; \     { \       __asm fld src_ \       __asm fistp dest_ \     } \     (dest) = dest_; \   } while (0)#else#  define QUANTFAC(rx)  adj43[rx]#  define ROUNDFAC 0.4054#  define XRPOW_FTOI(src,dest) ((dest) = (int)(src))#endif#undef MAXQUANTERRORFLOAT8 calc_sfb_noise(FLOAT8 *xr, FLOAT8 *xr34, int bw, int sf){  int j;  FLOAT8 xfsf=0;  FLOAT8 sfpow,sfpow34;  sfpow = POW20(sf+210); /*pow(2.0,sf/4.0); */  sfpow34  = IPOW20(sf+210); /*pow(sfpow,-3.0/4.0);*/  for ( j=0; j < bw ; ++j) {    int ix;    FLOAT8 temp;#if 0    if (xr34[j]*sfpow34 > IXMAX_VAL) return -1;    ix=floor( xr34[j]*sfpow34);    temp = fabs(xr[j])- pow43[ix]*sfpow;    temp *= temp;    if (ix < IXMAX_VAL) {      temp2 = fabs(xr[j])- pow43[ix+1]*sfpow;      temp2 *=temp2;      if (temp2<temp) {	temp=temp2;	++ix;      }    }#else    if (xr34[j]*sfpow34 > IXMAX_VAL) return -1;    temp = xr34[j]*sfpow34;    XRPOW_FTOI(temp, ix);    XRPOW_FTOI(temp + QUANTFAC(ix), ix);    temp = fabs(xr[j])- pow43[ix]*sfpow;    temp *= temp;    #endif    #ifdef MAXQUANTERROR    xfsf = Max(xfsf,temp);#else    xfsf += temp;#endif  }#ifdef MAXQUANTERROR  return xfsf;#else  return xfsf/bw;#endif}FLOAT8 calc_sfb_noise_ave(FLOAT8 *xr, FLOAT8 *xr34, int bw,int sf){  int j;  FLOAT8 xfsf=0, xfsf_p1=0, xfsf_m1=0;  FLOAT8 sfpow34,sfpow34_p1,sfpow34_m1;  FLOAT8 sfpow,sfpow_p1,sfpow_m1;  sfpow = POW20(sf+210); /*pow(2.0,sf/4.0); */  sfpow34  = IPOW20(sf+210); /*pow(sfpow,-3.0/4.0);*/  sfpow_m1 = sfpow*.8408964153;  /* pow(2,(sf-1)/4.0) */  sfpow34_m1 = sfpow34*1.13878863476;       /* .84089 ^ -3/4 */  sfpow_p1 = sfpow*1.189207115;    sfpow34_p1 = sfpow34*0.878126080187;  for ( j=0; j < bw ; ++j) {    int ix;    FLOAT8 temp,temp_p1,temp_m1;    if (xr34[j]*sfpow34_m1 > IXMAX_VAL) return -1;    temp = xr34[j]*sfpow34;    XRPOW_FTOI(temp, ix);    XRPOW_FTOI(temp + QUANTFAC(ix), ix);    temp = fabs(xr[j])- pow43[ix]*sfpow;    temp *= temp;    temp_p1 = xr34[j]*sfpow34_p1;    XRPOW_FTOI(temp_p1, ix);    XRPOW_FTOI(temp_p1 + QUANTFAC(ix), ix);    temp_p1 = fabs(xr[j])- pow43[ix]*sfpow_p1;    temp_p1 *= temp_p1;        temp_m1 = xr34[j]*sfpow34_m1;    XRPOW_FTOI(temp_m1, ix);    XRPOW_FTOI(temp_m1 + QUANTFAC(ix), ix);    temp_m1 = fabs(xr[j])- pow43[ix]*sfpow_m1;    temp_m1 *= temp_m1;#ifdef MAXQUANTERROR    xfsf = Max(xfsf,temp);    xfsf_p1 = Max(xfsf_p1,temp_p1);    xfsf_m1 = Max(xfsf_m1,temp_m1);#else    xfsf += temp;    xfsf_p1 += temp_p1;    xfsf_m1 += temp_m1;#endif  }  if (xfsf_p1>xfsf) xfsf = xfsf_p1;  if (xfsf_m1>xfsf) xfsf = xfsf_m1;#ifdef MAXQUANTERROR  return xfsf;#else  return xfsf/bw;#endif}int find_scalefac(FLOAT8 *xr,FLOAT8 *xr34,int sfb,		     FLOAT8 l3_xmin,int bw){  FLOAT8 xfsf;  int i,sf,sf_ok,delsf;  /* search will range from sf:  -209 -> 45  */  sf = -82;  delsf = 128;  sf_ok=10000;  for (i=0; i<7; i++) {    delsf /= 2;    //      xfsf = calc_sfb_noise(xr,xr34,bw,sf);    xfsf = calc_sfb_noise_ave(xr,xr34,bw,sf);    if (xfsf < 0) {      /* scalefactors too small */      sf += delsf;    }else{      if (sf_ok==10000) sf_ok=sf;        if (xfsf > l3_xmin)  {	/* distortion.  try a smaller scalefactor */	sf -= delsf;      }else{	sf_ok = sf;	sf += delsf;      }    }  }   assert(sf_ok!=10000);  //  assert(delsf==1);  /* when for loop goes up to 7 */  return sf;}/*    sfb=0..5  scalefac < 16     sfb>5     scalefac < 8    ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;    ol_sf =  (cod_info->global_gain-210.0);    ol_sf -= 8*cod_info->subblock_gain[i];    ol_sf -= ifqstep*scalefac[gr][ch].s[sfb][i];*/int compute_scalefacs_short(int sf[SBPSY_s][3],gr_info *cod_info,int scalefac[SBPSY_s][3],unsigned int sbg[3]){  int maxrange,maxrange1,maxrange2,maxover;  int sfb,i;  int ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;  maxover=0;  maxrange1 = 15;  maxrange2 = 7;  for (i=0; i<3; ++i) {    int maxsf1=0,maxsf2=0,minsf=1000;    /* see if we should use subblock gain */    for ( sfb = 0; sfb < SBPSY_s; sfb++ ) {      if (sfb < 6) {	if (-sf[sfb][i]>maxsf1) maxsf1 = -sf[sfb][i];      } else {	if (-sf[sfb][i]>maxsf2) maxsf2 = -sf[sfb][i];      }      if (-sf[sfb][i]<minsf) minsf = -sf[sfb][i];    }    /* boost subblock gain as little as possible so we can     * reach maxsf1 with scalefactors      * 8*sbg >= maxsf1        */    maxsf1 = Max(maxsf1-maxrange1*ifqstep,maxsf2-maxrange2*ifqstep);    sbg[i]=0;    if (minsf >0 ) sbg[i] = floor(.125*minsf + .001);    if (maxsf1 > 0)  sbg[i]  = Max(sbg[i],maxsf1/8 + (maxsf1 % 8 != 0));    if (sbg[i] > 7) sbg[i]=7;    for ( sfb = 0; sfb < SBPSY_s; sfb++ ) {      sf[sfb][i] += 8*sbg[i];      if (sf[sfb][i] < 0) {	maxrange = sfb < 6 ? maxrange1 : maxrange2;	scalefac[sfb][i]=-sf[sfb][i]/ifqstep + (-sf[sfb][i]%ifqstep != 0);	if (scalefac[sfb][i]>maxrange) scalefac[sfb][i]=maxrange;		if (-(sf[sfb][i] + scalefac[sfb][i]*ifqstep) >maxover)  {	  maxover=-(sf[sfb][i] + scalefac[sfb][i]*ifqstep);	}      }    }  }  return maxover;}int max_range_short[SBPSY_s]={15, 15, 15, 15, 15, 15 ,  7,    7,    7,    7,   7,     7 };int max_range_long[SBPSY_l]={15,   15,  15,  15,  15,  15,  15,  15,  15,  15,  15,   7,   7,   7,   7,   7,   7,   7,    7,    7,    7};/*	  ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;	  ol_sf =  (cod_info->global_gain-210.0);	  ol_sf -= ifqstep*scalefac[gr][ch].l[sfb];	  if (cod_info->preflag && sfb>=11) 	  ol_sf -= ifqstep*pretab[sfb];*/int compute_scalefacs_long(int sf[SBPSY_l],gr_info *cod_info,int scalefac[SBPSY_l]){  int sfb;  int maxover;  int ifqstep = ( cod_info->scalefac_scale == 0 ) ? 2 : 4;    if (cod_info->preflag)    for ( sfb = 11; sfb < SBPSY_l; sfb++ )       sf[sfb] += pretab[sfb]*ifqstep;  maxover=0;  for ( sfb = 0; sfb < SBPSY_l; sfb++ ) {    if (sf[sfb]<0) {      /* ifqstep*scalefac >= -sf[sfb], so round UP */      scalefac[sfb]=-sf[sfb]/ifqstep  + (-sf[sfb] % ifqstep != 0);      if (scalefac[sfb] > max_range_long[sfb]) scalefac[sfb]=max_range_long[sfb];            /* sf[sfb] should now be positive: */      if (  -(sf[sfb] + scalefac[sfb]*ifqstep)  > maxover) {	maxover = -(sf[sfb] + scalefac[sfb]*ifqstep);      }    }  }  return maxover;}    /************************************************************************ * * quantize and encode with the given scalefacs and global gain * * compute scalefactors, l3_enc, and return number of bits needed to encode * * ************************************************************************/voidVBR_quantize_granule(lame_global_flags *gfp,                FLOAT8 xr34[576], int l3_enc[2][2][576],		     III_psy_ratio *ratio,      III_psy_xmin l3_xmin,                III_scalefac_t scalefac[2][2],int gr, int ch){  lame_internal_flags *gfc=gfp->internal_flags;  int status;  gr_info *cod_info;    III_side_info_t * l3_side;  l3_side = &gfc->l3_side;  cod_info = &l3_side->gr[gr].ch[ch].tt;  /* encode scalefacs */  if ( gfp->version == 1 )     status=scale_bitcount(&scalefac[gr][ch], cod_info);  else    status=scale_bitcount_lsf(&scalefac[gr][ch], cod_info);  if (status!=0) {    cod_info->part2_3_length = LARGE_BITS;    return;  }    /* quantize xr34 */  cod_info->part2_3_length = count_bits(gfp,l3_enc[gr][ch],xr34,cod_info);  if (cod_info->part2_3_length >= LARGE_BITS) return;  cod_info->part2_3_length += cod_info->part2_length;    if (gfc->use_best_huffman==1 && cod_info->block_type != SHORT_TYPE) {    best_huffman_divide(gfc, gr, ch, cod_info, l3_enc[gr][ch]);  }  return;}  /************************************************************************ * * VBR_noise_shaping() * * compute scalefactors, l3_enc, and return number of bits needed to encode * * return code:    0   scalefactors were found with all noise < masking * *               n>0   scalefactors required too many bits.  global gain *                     was decreased by n *                     If n is large, we should probably recompute scalefacs *                     with a lower quality. * *               n<0   scalefactors used less than minbits. *                     global gain was increased by n.   *                     If n is large, might want to recompute scalefacs *                     with a higher quality setting? * ************************************************************************/intVBR_noise_shaping( lame_global_flags *gfp, FLOAT8 xr[576], III_psy_ratio *ratio, int l3_enc[2][2][576], int *ath_over, int minbits, int maxbits, III_scalefac_t scalefac[2][2], int gr,int ch){  lame_internal_flags *gfc=gfp->internal_flags;  int       start,end,bw,sfb,l, i,j, vbrmax;  III_scalefac_t vbrsf;  III_scalefac_t save_sf;  int maxover0,maxover1,maxover0p,maxover1p,maxover,mover;  int ifqstep;  III_psy_xmin l3_xmin;  III_side_info_t * l3_side;  gr_info *cod_info;    FLOAT8 xr34[576];  int shortblock;  int global_gain_adjust=0;  l3_side = &gfc->l3_side;  cod_info = &l3_side->gr[gr].ch[ch].tt;  shortblock = (cod_info->block_type == SHORT_TYPE);  *ath_over = calc_xmin( gfp,xr, ratio, cod_info, &l3_xmin);    for(i=0;i<576;i++) {    FLOAT8 temp=fabs(xr[i]);    xr34[i]=sqrt(sqrt(temp)*temp);  }#define MAX_SF_DELTA 4    vbrmax=-10000;  if (shortblock) {    for ( j=0, sfb = 0; sfb < SBMAX_s; sfb++ )  {      for ( i = 0; i < 3; i++ ) {	start = gfc->scalefac_band.s[ sfb ];	end   = gfc->scalefac_band.s[ sfb+1 ];	bw = end - start;	vbrsf.s[sfb][i] = find_scalefac(&xr[j],&xr34[j],sfb,					l3_xmin.s[sfb][i],bw);	j += bw;      }    }    for ( sfb = 0; sfb < SBMAX_s; sfb++ )  {      for ( i = 0; i < 3; i++ ) {	if (sfb>0) 	  vbrsf.s[sfb][i] = Min(vbrsf.s[sfb-1][i]+MAX_SF_DELTA,vbrsf.s[sfb][i]);	if (sfb< SBMAX_s-1) 	  vbrsf.s[sfb][i] = Min(vbrsf.s[sfb+1][i]+MAX_SF_DELTA,vbrsf.s[sfb][i]);      }    }    for ( j=0, sfb = 0; sfb < SBMAX_s; sfb++ )  {      for ( i = 0; i < 3; i++ ) {