psymodel.c

MP3编码的完整实现(源代码和使用例子都有)

    for ( j = 0; j < cw_lower_index; j++ )
      {	 /* calculate unpredictability measure cw */
	FLOAT an, a1, a2;
	FLOAT bn, b1, b2;
	FLOAT rn, r1, r2;
	FLOAT numre, numim, den;

	a2 = ax_sav[chn][1][j];
	b2 = bx_sav[chn][1][j];
	r2 = rx_sav[chn][1][j];
	a1 = ax_sav[chn][1][j] = ax_sav[chn][0][j];
	b1 = bx_sav[chn][1][j] = bx_sav[chn][0][j];
	r1 = rx_sav[chn][1][j] = rx_sav[chn][0][j];
	an = ax_sav[chn][0][j] = (*wsamp_l)[j];
	bn = bx_sav[chn][0][j] = j==0 ? (*wsamp_l)[0] : (*wsamp_l)[BLKSIZE-j];  
	rn = rx_sav[chn][0][j] = sqrt(energy[j]);

	{ /* square (x1,y1) */
	  if( r1 != 0 ) {
	    numre = (a1*b1);
	    numim = (a1*a1-b1*b1)*(FLOAT)0.5;
	    den = r1*r1;
	  } else {
	    numre = 1;
	    numim = 0;
	    den = 1;
	  }
	}
	
	{ /* multiply by (x2,-y2) */
	  if( r2 != 0 ) {
	    FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
	    FLOAT tmp1 = -a2*numre+tmp2;
	    numre =       -b2*numim+tmp2;
	    numim = tmp1;
	    den *= r2;
	  } else {
	    /* do nothing */
	  }
	}
	
	{ /* r-prime factor */
	  FLOAT tmp = (2*r1-r2)/den;
	  numre *= tmp;
	  numim *= tmp;
	}
	den=rn+fabs(2*r1-r2);
	if( den != 0 ) {
	  numre = (an+bn)*(FLOAT)0.5-numre;
	  numim = (an-bn)*(FLOAT)0.5-numim;
	  den = sqrt(numre*numre+numim*numim)/den;
	}
	cw[j] = den;
      }



    /**********************************************************************
     *     compute unpredicatibility of next 200 spectral lines            *
     **********************************************************************/ 
    for ( j = cw_lower_index; j < cw_upper_index; j += 4 )
      {/* calculate unpredictability measure cw */
	FLOAT rn, r1, r2;
	FLOAT numre, numim, den;
	
	k = (j+2) / 4; 
	
	{ /* square (x1,y1) */
	  r1 = energy_s[0][k];
	  if( r1 != 0 ) {
	    FLOAT a1 = (*wsamp_s)[0][k]; 
	    FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */
	    numre = (a1*b1);
	    numim = (a1*a1-b1*b1)*(FLOAT)0.5;
	    den = r1;
	    r1 = sqrt(r1);
	  } else {
	    numre = 1;
	    numim = 0;
	    den = 1;
	  }
	}
	
	
	{ /* multiply by (x2,-y2) */
	  r2 = energy_s[2][k];
	  if( r2 != 0 ) {
	    FLOAT a2 = (*wsamp_s)[2][k]; 
	    FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k];
	    
	    
	    FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
	    FLOAT tmp1 = -a2*numre+tmp2;
	    numre =       -b2*numim+tmp2;
	    numim = tmp1;
	    
	    r2 = sqrt(r2);
	    den *= r2;
	  } else {
	    /* do nothing */
	  }
	}
	
	{ /* r-prime factor */
	  FLOAT tmp = (2*r1-r2)/den;
	  numre *= tmp;
	  numim *= tmp;
	}
	
	rn = sqrt(energy_s[1][k]);
	den=rn+fabs(2*r1-r2);
	if( den != 0 ) {
	  FLOAT an = (*wsamp_s)[1][k]; 
	  FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k];
	  numre = (an+bn)*(FLOAT)0.5-numre;
	  numim = (an-bn)*(FLOAT)0.5-numim;
	  den = sqrt(numre*numre+numim*numim)/den;
	}
	
	cw[j+1] = cw[j+2] = cw[j+3] = cw[j] = den;
      }
    
#if 0
    for ( j = 14; j < HBLKSIZE-4; j += 4 )
      {/* calculate energy from short ffts */
	FLOAT8 tot,ave;
	k = (j+2) / 4; 
	for (tot=0, sblock=0; sblock < 3; sblock++)
	  tot+=energy_s[sblock][k];
	ave = energy[j+1]+ energy[j+2]+ energy[j+3]+ energy[j];
	ave /= 4.;
	/*
	  printf("energy / tot %i %5.2f   %e  %e\n",j,ave/(tot*16./3.),
	  ave,tot*16./3.);
	*/
	energy[j+1] = energy[j+2] = energy[j+3] =  energy[j]=tot;
      }
#endif
    
    
    
    
    
    
    
    
    /**********************************************************************
     *    Calculate the energy and the unpredictability in the threshold   *
     *    calculation partitions                                           *
     **********************************************************************/
#if 0
    for ( b = 0; b < CBANDS; b++ )
      {
	eb[b] = 0;
	cb[b] = 0;
      }
    for ( j = 0; j < HBLKSIZE; j++ )
      {
	int tp = partition_l[j];
	
	if ( tp >= 0 )
	  {
	    eb[tp] += energy[j];
	    cb[tp] += cw[j] * energy[j];
	  }
	assert(tp<npart_l_orig);
      }
#else
    b = 0;
    for (j = 0; j < cw_upper_index;)
      {
	FLOAT8 ebb, cbb;
	int i;

	ebb = energy[j];
	cbb = energy[j] * cw[j];
	j++;

	for (i = numlines_l[b] - 1; i > 0; i--)
	  {
	    ebb += energy[j];
	    cbb += energy[j] * cw[j];
	    j++;
	  }
	eb[b] = ebb;
	cb[b] = cbb;
	b++;
      }

    for (; b < npart_l_orig; b++ )
      {
	int i;
	FLOAT8 ebb = energy[j++];

	for (i = numlines_l[b] - 1; i > 0; i--)
	  {
	    ebb += energy[j++];
	  }
	eb[b] = ebb;
	cb[b] = ebb * 0.4;
      }
#endif

    /**********************************************************************
     *      convolve the partitioned energy and unpredictability           *
     *      with the spreading function, s3_l[b][k]                        *
     ******************************************************************** */
    pe[chn] = 0;		/*  calculate percetual entropy */
    for ( b = 0;b < npart_l; b++ )
      {
	FLOAT8 tbb,ecb,ctb;
	FLOAT8 temp_1; /* BUG of IS */

	ecb = 0;
	ctb = 0;
	for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ )
	  {
	    ecb += s3_l[b][k] * eb[k];	/* sprdngf for Layer III */
	    ctb += s3_l[b][k] * cb[k];
	  }

	/* calculate the tonality of each threshold calculation partition */
	/* calculate the SNR in each threshhold calculation partition */

	tbb = ecb;
	if (tbb != 0)
	  {
	    tbb = ctb / tbb;
	    if (tbb <= 0.04875584301)
	      {
		tbb = exp(-LN_TO_LOG10 * (TMN - NMT));
	      }
	    else if (tbb > 0.4989003827)
	      {
		tbb = 1;
	      }
	    else
	      {
		tbb = log(tbb);
		tbb = exp(((TMN - NMT)*(LN_TO_LOG10*0.299))
			+ ((TMN - NMT)*(LN_TO_LOG10*0.43 ))*tbb);  /* conv1=-0.299, conv2=-0.43 */
	      }
	  }

	tbb = Min(minval[b], tbb);
	ecb *= tbb;

	/* pre-echo control */
	/* rpelev=2.0, rpelev2=16.0 */
	temp_1 = Min(ecb, Min(rpelev*nb_1[chn][b],rpelev2*nb_2[chn][b]) );
	thr[b] = Max( qthr_l[b], temp_1 ); 
	nb_2[chn][b] = nb_1[chn][b];
	nb_1[chn][b] = ecb;

	/* note: all surges in PE are because of the above pre-echo formula
	 * for temp_1.  it this is not used, PE is always around 600
	 */

	if (thr[b] < eb[b])
	  {
	    /* there's no non sound portition, because thr[b] is
	     maximum of qthr_l and temp_1 */
	    pe[chn] -= numlines_l[b] * log(thr[b] / eb[b]);
	  }
      }


#ifdef HAVEGTK
    if (gfp->gtkflag) {
      FLOAT mn,mx,ma=0,mb=0,mc=0;

      for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
      {
        ma += energy_s[0][j];
        mb += energy_s[1][j];
        mc += energy_s[2][j];
      }
      mn = Min(ma,mb);
      mn = Min(mn,mc);
      mx = Max(ma,mb);
      mx = Max(mx,mc);

      pinfo->ers[gr_out][chn]=ers_save[chn];
      ers_save[chn]=mx/(1e-12+mn);
      pinfo->pe[gr_out][chn]=pe_save[chn];
      pe_save[chn]=pe[chn];
    }
#endif
    
    /*************************************************************** 
     * determine the block type (window type) based on L & R channels
     * 
     ***************************************************************/
    if (chn<2) {
      if (gfp->no_short_blocks){
	uselongblock[chn]=1;
      } else {
	/* tuned for t1.wav.  doesnt effect most other samples */
	if (pe[chn] > 3000) {
	  uselongblock[chn]=0;
	} else { 
	  FLOAT mn,mx,ma=0,mb=0,mc=0;
	
	  for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
	  {
	      ma += energy_s[0][j];
	      mb += energy_s[1][j];
	      mc += energy_s[2][j];
	  }
	  mn = Min(ma,mb);
	  mn = Min(mn,mc);
	  mx = Max(ma,mb);
	  mx = Max(mx,mc);

	  uselongblock[chn] = 1;
	  
	  if ( mx > 30*mn ) 
	  {/* big surge of energy - always use short blocks */
	    uselongblock[chn] = 0;
	  } 
	  else if ((mx > 10*mn) && (pe[chn] > 1000))
	  {/* medium surge, medium pe - use short blocks */
	    uselongblock[chn] = 0;
	  }
	} 
      }
    }



    /*************************************************************** 
     * compute masking thresholds for both short and long blocks
     ***************************************************************/
    /* longblock threshold calculation (part 2) */
    for ( sb = 0; sb < SBPSY_l; sb++ )
      {
	FLOAT8 enn = w1_l[sb] * eb[bu_l[sb]] + w2_l[sb] * eb[bo_l[sb]];
	FLOAT8 thmm = w1_l[sb] *thr[bu_l[sb]] + w2_l[sb] * thr[bo_l[sb]];
	for ( b = bu_l[sb]+1; b < bo_l[sb]; b++ )
	  {
	    enn  += eb[b];
	    thmm += thr[b];
	  }
	en[chn].l[sb] = enn;
	thm[chn].l[sb] = thmm;
      }
    
    
    /* threshold calculation for short blocks */
    for ( sblock = 0; sblock < 3; sblock++ )
      {
	j = 0;
	for ( b = 0; b < npart_s_orig; b++ )
	  {
	    int i;
	    FLOAT ecb = energy_s[sblock][j++];
	    for (i = numlines_s[b]; i > 0; i--)
	      {
		ecb += energy_s[sblock][j++];
	      }
	    eb[b] = ecb;
	  }

	for ( b = 0; b < npart_s; b++ )
	  {
	    FLOAT8 ecb = 0;
	    for ( k = s3ind_s[b][0]; k <= s3ind_s[b][1]; k++ )
	      {
		ecb += s3_s[b][k] * eb[k];
	      }
	    thr[b] = Max (qthr_s[b], ecb);
	  }

	for ( sb = 0; sb < SBPSY_s; sb++ )
	  {
	    FLOAT8 enn  = w1_s[sb] * eb[bu_s[sb]] + w2_s[sb] * eb[bo_s[sb]];
	    FLOAT8 thmm = w1_s[sb] *thr[bu_s[sb]] + w2_s[sb] * thr[bo_s[sb]];
	    for ( b = bu_s[sb]+1; b < bo_s[sb]; b++ )
	      {
		enn  += eb[b];
		thmm += thr[b];
	      }
	    en[chn].s[sb][sblock] = enn;
	    thm[chn].s[sb][sblock] = thmm;
	  }
      }
  } /* end loop over chn */


  /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */
  if ( numchn==4 /* mid/side and r/l */) {
    FLOAT8 rside,rmid,mld;
    int chmid=2,chside=3; 
    
    for ( sb = 0; sb < SBPSY_l; sb++ ) {
      /* use this fix if L & R masking differs by 2db or less */
      /* if db = 10*log10(x2/x1) < 2 */
      /* if (x2 < 1.58*x1) { */
      if (thm[0].l[sb] <= 1.58*thm[1].l[sb]
	  && thm[1].l[sb] <= 1.58*thm[0].l[sb]) {

	mld = mld_l[sb]*en[chside].l[sb];
	rmid = Max(thm[chmid].l[sb], Min(thm[chside].l[sb],mld));

	mld = mld_l[sb]*en[chmid].l[sb];
	rside = Max(thm[chside].l[sb],Min(thm[chmid].l[sb],mld));

	thm[chmid].l[sb]=rmid;
	thm[chside].l[sb]=rside;
      }
    }
    for ( sb = 0; sb < SBPSY_s; sb++ ) {
      for ( sblock = 0; sblock < 3; sblock++ ) {
	if (thm[0].s[sb][sblock] <= 1.58*thm[1].s[sb][sblock]
	    && thm[1].s[sb][sblock] <= 1.58*thm[0].s[sb][sblock]) {

	  mld = mld_s[sb]*en[chside].s[sb][sblock];
	  rmid = Max(thm[chmid].s[sb][sblock],Min(thm[chside].s[sb][sblock],mld));