psymodel.c

MP3编码程序和资料

      if (sfreq == freq_tp/freq_scale )
	{
	  cbmax = cbmax_tp;
	  for(i=0,k2=0;i<cbmax_tp;i++)
	    {
	      j = (int) *p++;
	      numlines_l[i] = (int) *p++;
	      minval[i] = exp(-((*p++) ) * LN_TO_LOG10);
	      /* qthr_l[i] = *p++ */ p++;
	      /* norm_l[i] = *p++*/ p++;
	      /* bval_l[i] = *p++; */ p++;
	      if (j!=i)
		{
		  ERRORF("1. please check \"psy_data\"");
		  return -1;
		}
	    }
	}
      else
	p += cbmax_tp * 6;
    }

  *npart_l_orig = cbmax;

  /* Read short block data */
  for(loop=0;loop<6;loop++)
    {
      freq_tp = *p++;
      cbmax_tp = (int) *p++;
      cbmax_tp++;
      if (sfreq == freq_tp/freq_scale )
	{
	  cbmax = cbmax_tp;
	  for(i=0,k2=0;i<cbmax_tp;i++)
	    {
	      j = (int) *p++;
	      numlines_s[i] = (int) *p++;
	      /* qthr_s[i] = *p++*/  p++;         
	      /* norm_s[i] =*p++ */ p++;         
	      SNR[i] = *p++;            
	      /* bval_s[i] = *p++ */ p++;
	      if (j!=i)
		{
		  ERRORF("3. please check \"psy_data\"");
		  return -1;
		}
	    }
	}
      else
	p += cbmax_tp * 6;
    }
  *npart_s_orig = cbmax;

  /* MPEG1 SNR_s data is given in db, convert to energy */
  if (gfp->version == 1) {
    for ( i = 0;i < *npart_s_orig; i++ ) {
      SNR[i]=exp( (FLOAT8) SNR[i] * LN_TO_LOG10 );
    }
  }




  /* Read long block data for converting threshold calculation 
     partitions to scale factor bands */

  for(loop=0;loop<6;loop++)
    {
      freq_tp = *p++;
      sbmax =  (int) *p++;
      sbmax++;

      if (sfreq == freq_tp/freq_scale)
	{
	  for(i=0;i<sbmax;i++)
	    {
	      j = (int) *p++;
	      p++;             
	      bu_l[i] = (int) *p++;
	      bo_l[i] = (int) *p++;
	      w1_l[i] = (FLOAT8) *p++;
	      w2_l[i] = (FLOAT8) *p++;
	      if (j!=i)
		{ ERRORF("30:please check \"psy_data\"\n");
		return -1;
		}

	      if (i!=0)
		if ( (fabs(1.0-w1_l[i]-w2_l[i-1]) > 0.01 ) )
		  {
		    ERRORF("31l: please check \"psy_data.\"\n"
                           "w1,w2: %f %f \n",w1_l[i],w2_l[i-1]);
		    return -1;
		  }
	    }
	}
      else
	p += sbmax * 6;
    }

  /* Read short block data for converting threshold calculation 
     partitions to scale factor bands */

  for(loop=0;loop<6;loop++)
    {
      freq_tp = *p++;
      sbmax = (int) *p++;
      sbmax++;

      if (sfreq == freq_tp/freq_scale)
	{
	  for(i=0;i<sbmax;i++)
	    {
	      j = (int) *p++;
	      p++;
	      bu_s[i] = (int) *p++;
	      bo_s[i] = (int) *p++;
	      w1_s[i] = *p++;
	      w2_s[i] = *p++;
	      if (j!=i)
		{ ERRORF("30:please check \"psy_data\"\n");
		return -1;
		}

	      if (i!=0)
		if ( (fabs(1.0-w1_s[i]-w2_s[i-1]) > 0.01 ) )
		  { 
                  ERRORF("31s: please check \"psy_data.\"\n"
                         "w1,w2: %f %f \n",w1_s[i],w2_s[i-1]);
		  return -1;
		  }
	    }
	}
      else
	p += sbmax * 6;
    }

  /******************************************************************/
  /* done reading table data */
  /******************************************************************/




#undef NOTABLES
#ifdef NOTABLES
  /* compute numlines */
  j=0;
  for(i=0;i<CBANDS;i++)
    {
      FLOAT8 ji, bark1,bark2,delbark=.34;
      int k,j2;

      j2 = j;
      j2 = Min(j2,BLKSIZE/2);
      
      do {
	/* find smallest j2 >= j so that  (bark - bark_l[i-1]) < delbark */
	ji = j;
	bark1 = freq2bark(sfreq*ji/BLKSIZE);
	
	++j2;
	ji = j2;
	bark2  = freq2bark(sfreq*ji/BLKSIZE);
      } while ((bark2 - bark1) < delbark  && j2<=BLKSIZE/2);

      /*
      DEBUGF("%i old n=%i  %f old numlines:  %i   new=%i (%i,%i) (%f,%f) \n",
i,*npart_l_orig,freq,numlines_l[i],j2-j,j,j2-1,bark1,bark2);
      */
      for (k=j; k<j2; ++k)
	partition[k]=i;
      numlines_l[i]=(j2-j);
      j = j2;
      if (j > BLKSIZE/2) break;
    }
  *npart_l_orig = i;

  /* compute which partition bands are in which scalefactor bands */
  { int i1,i2,sfb,start,end;
    FLOAT8 freq1,freq2;
    for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
      start = gfc->scalefac_band.l[ sfb ];
      end   = gfc->scalefac_band.l[ sfb+1 ];
      freq1 = sfreq*(start-.5)/(2*576);
      freq2 = sfreq*(end-1+.5)/(2*576);
		     
      i1 = floor(.5 + BLKSIZE*freq1/sfreq);
      if (i1<0) i1=0;
      i2 = floor(.5 + BLKSIZE*freq2/sfreq);
      if (i2>BLKSIZE/2) i2=BLKSIZE/2;

      DEBUGF("old: (%i,%i)  new: (%i,%i) %i %i \n",bu_l[sfb],bo_l[sfb],
	     partition[i1],partition[i2],i1,i2);

      w1_l[sfb]=.5;
      w2_l[sfb]=.5;
      bu_l[sfb]=partition[i1];
      bo_l[sfb]=partition[i2];

    }
  }

#endif


  /* compute bark value and ATH of each critical band */
  j=0;
  for(i=0;i<*npart_l_orig;i++)
    {
      FLOAT8 ji, freq, mval, bark;
      int k;

      ji = j;
      bark = freq2bark(sfreq*ji/BLKSIZE);

      ji = j + (numlines_l[i]-1);
      bark = .5*(bark + freq2bark(sfreq*ji/BLKSIZE));
      bval_l[i]=bark;
      //      j += numlines_l[i];


      gfc->ATH_partitionbands[i]=1e99;
      for (k=0; k < numlines_l[i]; ++k) {
	FLOAT8 freq = sfreq*j/(1000.0*BLKSIZE);
	assert( freq < 25 );
	//	freq = Min(.1,freq);    /* ignore ATH below 100hz */
	freq= ATHformula(freq);  
	freq += -20; /* scale to FFT units */
	freq = pow( 10.0, freq/10.0 );  /* convert from db -> energy */
	freq *= numlines_l[i];
	gfc->ATH_partitionbands[i]=Min(gfc->ATH_partitionbands[i],freq);
	++j;
      }


      /*  minval formula needs work
      ji = j;
      freq = sfreq*ji/BLKSIZE;
      mval = Max(27.0 - freq/50.0, 0.0);
      //      mval = exp(-mval * LN_TO_LOG10);

      DEBUGF("%2i old minval=%f  new = %f  ",
	     i,-log(minval[i])/LN_TO_LOG10,mval);
      if (mval < minval[i]) 
      DEBUGF("less masking than ISO tables \n");
      else
      DEBUGF("more masking than ISO tables \n");
      */
    }



  /************************************************************************
   * Now compute the spreading function, s[j][i], the value of the spread-*
   * ing function, centered at band j, for band i, store for later use    *
   ************************************************************************/
  /* i.e.: sum over j to spread into signal barkval=i  
     NOTE: i and j are used opposite as in the ISO docs */
  for(i=0;i<*npart_l_orig;i++)
    {
      FLOAT8 tempx,x,tempy,temp;
      for(j=0;j<*npart_l_orig;j++)
	{
          if (i>=j) tempx = (bval_l[i] - bval_l[j])*3.0;
	  else    tempx = (bval_l[i] - bval_l[j])*1.5; 

	  if(tempx>=0.5 && tempx<=2.5)
	    {
	      temp = tempx - 0.5;
	      x = 8.0 * (temp*temp - 2.0 * temp);
	    }
	  else x = 0.0;
	  tempx += 0.474;
	  tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);

#ifdef NEWS3
	  if (j>=i) tempy = (bval_l[j] - bval_l[i])*(-15);
	  else    tempy = (bval_l[j] - bval_l[i])*25;
	  x=0; 
#endif
	  /*
	  if ((i==cbmax/2)  && (fabs(bval_l[j] - bval_l[i])) < 3) {
	    DEBUGF("bark=%f   x+tempy = %f  \n",bval_l[j] - bval_l[i],x+tempy);
	  }
	  */

	  if (tempy <= -60.0) s3_l[i][j] = 0.0;
	  else                s3_l[i][j] = exp( (x + tempy)*LN_TO_LOG10 ); 
	}
    }


  /************************************************************************/
  /* SHORT BLOCKS */
  /************************************************************************/

#ifdef NOTABLES
  /* compute numlines */
  j=0;
  for(i=0;i<CBANDS;i++)
    {
      FLOAT8 ji, bark1,bark2,delbark=.34;
      int k,j2;

      j2 = j;
      j2 = Min(j2,BLKSIZE_s/2);
      
      do {
	/* find smallest j2 >= j so that  (bark - bark_s[i-1]) < delbark */
	ji = j;
	bark1  = freq2bark(sfreq*ji/BLKSIZE_s);
	
	++j2;
	ji = j2;
	bark2  = freq2bark(sfreq*ji/BLKSIZE_s);

      } while ((bark2 - bark1) < delbark  && j2<=BLKSIZE_s/2);

      /*
      DEBUGF("%i old n=%i  %f old numlines:  %i   new=%i (%i,%i) (%f,%f) \n",
i,*npart_s_orig,freq,numlines_s[i],j2-j,j,j2-1,bark1,bark2);
      */
      for (k=j; k<j2; ++k)
	partition[k]=i;
      numlines_s[i]=(j2-j);
      j = j2;
      if (j > BLKSIZE_s/2) break;
    }
  *npart_s_orig = i;

  /* compute which partition bands are in which scalefactor bands */
  { int i1,i2,sfb,start,end;
    FLOAT8 freq1,freq2;
    for ( sfb = 0; sfb < SBMAX_s; sfb++ ) {
      start = gfc->scalefac_band.s[ sfb ];
      end   = gfc->scalefac_band.s[ sfb+1 ];
      freq1 = sfreq*(start-.5)/(2*192);
      freq2 = sfreq*(end-1+.5)/(2*192);
		     
      i1 = floor(.5 + BLKSIZE_s*freq1/sfreq);
      if (i1<0) i1=0;
      i2 = floor(.5 + BLKSIZE_s*freq2/sfreq);
      if (i2>BLKSIZE_s/2) i2=BLKSIZE_s/2;

      DEBUGF("old: (%i,%i)  new: (%i,%i) %i %i \n",bu_s[sfb],bo_s[sfb],
	     partition[i1],partition[i2],i1,i2);

      w1_s[sfb]=.5;
      w2_s[sfb]=.5;
      bu_s[sfb]=partition[i1];
      bo_s[sfb]=partition[i2];

    }
  }

#endif



  /* compute bark values of each critical band */
  j = 0;
  for(i=0;i<*npart_s_orig;i++)
    {
      FLOAT8 ji, bark, freq, snr;
      ji = j;
      bark = freq2bark(sfreq*ji/BLKSIZE_s);

      ji = j + numlines_s[i] -1;
      bark = .5*(bark+freq2bark(sfreq*ji/BLKSIZE_s));
      /*
      DEBUGF("%i %i bval_s = %f  %f  numlines=%i  formula=%f \n",i,j,bval_s[i],freq,numlines_s[i],bark);
      */
      bval_s[i]=bark;
      j += numlines_s[i];

      /* SNR formula needs work 
      ji = j;
      freq = sfreq*ji/BLKSIZE;
      snr  = .14 + freq/80000;
      DEBUGF("%2i old SNR=%f  new = %f \n ",
	     i,SNR[i],snr);
      */
    }






  /************************************************************************
   * Now compute the spreading function, s[j][i], the value of the spread-*
   * ing function, centered at band j, for band i, store for later use    *
   ************************************************************************/
  for(i=0;i<*npart_s_orig;i++)
    {
      FLOAT8 tempx,x,tempy,temp;
      for(j=0;j<*npart_s_orig;j++)
	{
          if (i>=j) tempx = (bval_s[i] - bval_s[j])*3.0;
	  else    tempx = (bval_s[i] - bval_s[j])*1.5; 

	  if(tempx>=0.5 && tempx<=2.5)
	    {
	      temp = tempx - 0.5;
	      x = 8.0 * (temp*temp - 2.0 * temp);
	    }
	  else x = 0.0;
	  tempx += 0.474;
	  tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
#ifdef NEWS3
	  if (j>=i) tempy = (bval_s[j] - bval_s[i])*(-15);
	  else    tempy = (bval_s[j] - bval_s[i])*25;
	  x=0; 
#endif
	  if (tempy <= -60.0) s3_s[i][j] = 0.0;
	  else                s3_s[i][j] = exp( (x + tempy)*LN_TO_LOG10 );
	}
    }




  /* compute: */
  /* npart_l_orig   = number of partition bands before convolution */
  /* npart_l  = number of partition bands after convolution */

  assert(*npart_l_orig <=CBANDS);
  assert(*npart_s_orig<=CBANDS);

  
  *npart_l=bo_l[SBPSY_l-1]+1;
  *npart_s=bo_s[SBPSY_s-1]+1;
  
  /* if npart_l = npart_l_orig + 1, we can fix that below.  else: */
  assert(*npart_l <= *npart_l_orig+1);
  assert(*npart_s <= *npart_s_orig+1);


  /* MPEG2 tables are screwed up 
   * the mapping from paritition bands to scalefactor bands will use
   * more paritition bands than we have.  
   * So we will not compute these fictitious partition bands by reducing
   * npart_l below.  */
  if (*npart_l > *npart_l_orig) {
    *npart_l=*npart_l_orig;
    bo_l[SBPSY_l-1]=(*npart_l)-1;
    w2_l[SBPSY_l-1]=1.0;
  }

  if (*npart_s > *npart_s_orig) {
    *npart_s=*npart_s_orig;
    bo_s[SBPSY_s-1]=(*npart_s)-1;
    w2_s[SBPSY_s-1]=1.0;
  }


    /* setup stereo demasking thresholds */
    /* formula reverse enginerred from plot in paper */
    for ( i = 0; i < SBPSY_s; i++ ) {
      FLOAT8 arg,mld;
      arg = freq2bark(sfreq*gfc->scalefac_band.s[i]/(2*192));
      arg = (Min(arg, 15.5)/15.5);

      mld = 1.25*(1-cos(PI*arg))-2.5;
      gfc->mld_s[i] = pow(10.0,mld);
    }
    for ( i = 0; i < SBPSY_l; i++ ) {
      FLOAT8 arg,mld;
      arg = freq2bark(sfreq*gfc->scalefac_band.l[i]/(2*576));
      arg = (Min(arg, 15.5)/15.5);

      mld = 1.25*(1-cos(PI*arg))-2.5;
      gfc->mld_l[i] = pow(10.0,mld);
    }

  
  return 0;
}