psymodel.c

音频编码

#define I1LIMIT 8   /* as in if(i>8)  */ 
#define I2LIMIT 23  /* as in if(i>24) -> changed 23 */ 
#define MLIMIT  15  /* as in if(m<15) */ 

static FLOAT ma_max_i1;
static FLOAT ma_max_i2;
static FLOAT ma_max_m;



static void init_mask_add_max_values(void)
{
    ma_max_i1 = pow(10,(I1LIMIT+1)/16.0);
    ma_max_i2 = pow(10,(I2LIMIT+1)/16.0);
    ma_max_m  = pow(10,(MLIMIT)/10.0);
}






/* addition of simultaneous masking   Naoki Shibata 2000/7 */
inline static FLOAT mask_add(FLOAT m1,FLOAT m2,int k,int b, lame_internal_flags * const gfc)
{
  static const FLOAT table1[] = {
    3.3246 *3.3246 ,3.23837*3.23837,3.15437*3.15437,3.00412*3.00412,2.86103*2.86103,2.65407*2.65407,2.46209*2.46209,2.284  *2.284  ,
    2.11879*2.11879,1.96552*1.96552,1.82335*1.82335,1.69146*1.69146,1.56911*1.56911,1.46658*1.46658,1.37074*1.37074,1.31036*1.31036,
    1.25264*1.25264,1.20648*1.20648,1.16203*1.16203,1.12765*1.12765,1.09428*1.09428,1.0659 *1.0659 ,1.03826*1.03826,1.01895*1.01895,
    1
  };

  static const FLOAT table2[] = {
    1.33352*1.33352,1.35879*1.35879,1.38454*1.38454,1.39497*1.39497,1.40548*1.40548,1.3537 *1.3537 ,1.30382*1.30382,1.22321*1.22321,
    1.14758*1.14758,
    1
  };

  static const FLOAT table3[] = {
    2.35364*2.35364,2.29259*2.29259,2.23313*2.23313,2.12675*2.12675,2.02545*2.02545,1.87894*1.87894,1.74303*1.74303,1.61695*1.61695,
    1.49999*1.49999,1.39148*1.39148,1.29083*1.29083,1.19746*1.19746,1.11084*1.11084,1.03826*1.03826
  };


  int i;
  FLOAT ratio;


  if (m2 > m1) {
      if (m2 < (m1*ma_max_i2))
        ratio = m2/m1;
      else
        return (m1+m2);
  } else {
      if (m1 >= (m2*ma_max_i2))
          return (m1+m2);
      ratio = m1/m2;
  }
  /*i = abs(10*log10(m2 / m1)/10*16);
  m = 10*log10((m1+m2)/gfc->ATH->cb[k]);*/


  /* Should always be true, just checking */
  assert(m1>=0);
  assert(m2>=0);
  assert(gfc->ATH->cb[k]>=0);


  m1 += m2;

  if ((unsigned int)(b+3) <= 3+3) {  /* approximately, 1 bark = 3 partitions */
      /* 65% of the cases */
      /* originally 'if(i > 8)' */
      if (ratio >= ma_max_i1) {
	  /* 43% of the total */
	  return m1;
      }

      /* 22% of the total */
      i = FAST_LOG10_X(ratio,16.0);
      return m1*table2[i];
  }

  /* m<15 equ log10((m1+m2)/gfc->ATH->cb[k])<1.5
   * equ (m1+m2)/gfc->ATH->cb[k]<10^1.5
   * equ (m1+m2)<10^1.5 * gfc->ATH->cb[k]
   */

  i = FAST_LOG10_X(ratio, 16.0);
  m2 = gfc->ATH->cb[k]*gfc->ATH->adjust;
  if (m1 < ma_max_m*m2)  {
      /* 3% of the total */
      /* Originally if (m > 0) { */
      if (m1 > m2) {
	  FLOAT f, r;

	  f = 1.0;
	  if (i <= 13) f = table3[i];

	  r = FAST_LOG10_X(m1 / m2, 10.0/15.0);
	  return m1 * ((table1[i]-f)*r+f);
      }

      if (i > 13) return m1;

      return m1*table3[i];
  }


  /* 10% of total */
  return m1*table1[i];
}



static inline FLOAT NS_INTERP(FLOAT x, FLOAT y, FLOAT r)
{
    /* was pow((x),(r))*pow((y),1-(r))*/
    if(r==1.0)
        return x;              /* 99.7% of the time */
    if(r==0.0)
	return y;
    if(y>0.0)
        return pow(x/y,r)*y;   /* rest of the time */ 
    return 0.0;                /* never happens */ 
}



static void nsPsy2dataRead(
    FILE *fp,
    FLOAT *eb2,
    FLOAT *eb,
    int chn,
    int npart_l
    )
{
    int b;
    for(;;) {
	static const char chname[] = {'L','R','M','S'};
	char c;

	fscanf(fp, "%c",&c);
	for (b=0; b < npart_l; b++) {
	    double e;
	    fscanf(fp, "%lf",&e);
	    eb2[b] = e;
	}

	if (feof(fp)) abort();
	if (c == chname[chn]) break;
	abort();
    }

    eb2[62] = eb2[61];
    for (b=0; b < npart_l; b++ )
	eb2[b] = eb2[b] * eb[b];
}

static FLOAT
pecalc_s(
    III_psy_ratio *mr,
    FLOAT masking_lower
    )
{
    FLOAT pe_s;
    const static FLOAT regcoef_s[] = {
	11.8, /* these values are tuned only for 44.1kHz... */
	13.6,
	17.2,
	32,
	46.5,
	51.3,
	57.5,
	67.1,
	71.5,
	84.6,
	97.6,
	130,
/*	255.8 */
    };
    int sb, sblock;

    pe_s = 1236.28/4;
    for (sb = 0; sb < SBMAX_s - 1; sb++) {
	for (sblock=0;sblock<3;sblock++) {
	    FLOAT x;
	    if (regcoef_s[sb] == 0.0
		|| mr->thm.s[sb][sblock] <= 0.0
		|| mr->en.s[sb][sblock]
		<= (x = mr->thm.s[sb][sblock] * masking_lower))
		continue;

	    if (mr->en.s[sb][sblock] > x*1e10)
		pe_s += regcoef_s[sb] * (10.0 * LOG10);
	    else
		pe_s += regcoef_s[sb] * FAST_LOG10(mr->en.s[sb][sblock] / x);
	}
    }

    return pe_s;
}

static FLOAT
pecalc_l(
    III_psy_ratio *mr,
    FLOAT masking_lower
    )
{
    FLOAT pe_l;
    const static FLOAT regcoef_l[] = {
	6.8, /* these values are tuned only for 44.1kHz... */
	5.8,
	5.8,
	6.4,
	6.5,
	9.9,
	12.1,
	14.4,
	15,
	18.9,
	21.6,
	26.9,
	34.2,
	40.2,
	46.8,
	56.5,
	60.7,
	73.9,
	85.7,
	93.4,
	126.1,
/*	241.3 */
    };
    int sb;

    pe_l = 1124.23/4;
    for (sb = 0; sb < SBMAX_l - 1; sb++) {
	FLOAT x;
	if (mr->thm.l[sb] <= 0.0
	    || mr->en.l[sb] <= (x = mr->thm.l[sb]*masking_lower))
	    continue;

	if (mr->en.l[sb] > x*1e10)
	    pe_l += regcoef_l[sb] * (10.0 * LOG10);
	else
	    pe_l += regcoef_l[sb] * FAST_LOG10(mr->en.l[sb] / x);
    }

    return pe_l;
}





int L3psycho_anal_ns( lame_global_flags * gfp,
                    const sample_t *buffer[2], int gr_out, 
                    FLOAT *ms_ratio,
                    FLOAT *ms_ratio_next,
		    III_psy_ratio masking_ratio[2][2],
		    III_psy_ratio masking_MS_ratio[2][2],
		    FLOAT percep_entropy[2],FLOAT percep_MS_entropy[2], 
		    FLOAT energy[4], 
                    int blocktype_d[2])
{
/* to get a good cache performance, one has to think about
 * the sequence, in which the variables are used.  
 * (Note: these static variables have been moved to the gfc-> struct,
 * and their order in memory is layed out in util.h)
 */
    lame_internal_flags *gfc=gfp->internal_flags;

    /* fft and energy calculation   */
    FLOAT wsamp_L[2][BLKSIZE];
    FLOAT wsamp_S[2][3][BLKSIZE_s];

    /* block type  */
    int blocktype[2],uselongblock[2];

    /* usual variables like loop indices, etc..    */
    int numchn, chn;
    int b, i, j, k;
    int	sb,sblock;

    /* variables used for --nspsytune */
    FLOAT ns_hpfsmpl[2][576];
    FLOAT pcfact;

    numchn = gfc->channels_out;
    /* chn=2 and 3 = Mid and Side channels */
    if (gfp->mode == JOINT_STEREO) numchn=4;

    if (gfp->VBR==vbr_off) pcfact = gfc->ResvMax == 0 ? 0 : ((FLOAT)gfc->ResvSize)/gfc->ResvMax*0.5;
    else if (gfp->VBR == vbr_rh  ||  gfp->VBR == vbr_mtrh  ||  gfp->VBR == vbr_mt) {
	    /*static const FLOAT pcQns[10]={1.0,1.0,1.0,0.8,0.6,0.5,0.4,0.3,0.2,0.1};
	    pcfact = pcQns[gfp->VBR_q];*/
	    pcfact = 0.6;
    } else pcfact = 1.0;

    /**********************************************************************
     *  Apply HPF of fs/4 to the input signal.
     *  This is used for attack detection / handling.
     **********************************************************************/
    /* Don't copy the input buffer into a temporary buffer */
    /* unroll the loop 2 times */
    for(chn=0;chn<gfc->channels_out;chn++) {
	static const FLOAT fircoef[] = {
	    -8.65163e-18*2, -0.00851586*2, -6.74764e-18*2, 0.0209036*2,
	    -3.36639e-17*2, -0.0438162 *2, -1.54175e-17*2, 0.0931738*2,
	    -5.52212e-17*2, -0.313819  *2
	};

	/* apply high pass filter of fs/4 */
	const sample_t * const firbuf = &buffer[chn][576-350-NSFIRLEN+192];
	for (i=0;i<576;i++) {
	    FLOAT sum1, sum2;
	    sum1 = firbuf[i + 10];
	    sum2 = 0.0;
	    for (j=0;j<(NSFIRLEN-1)/2;j+=2) {
		sum1 += fircoef[j  ] * (firbuf[i+j  ]+firbuf[i+NSFIRLEN-j  ]);
		sum2 += fircoef[j+1] * (firbuf[i+j+1]+firbuf[i+NSFIRLEN-j-1]);
	    }
	    ns_hpfsmpl[chn][i] = sum1 + sum2;
	}
	masking_ratio    [gr_out] [chn]  .en  = gfc -> en  [chn];
	masking_ratio    [gr_out] [chn]  .thm = gfc -> thm [chn];
	if (numchn > 2) {
	    /* MS maskings  */
	    /*percep_MS_entropy         [chn-2]     = gfc -> pe  [chn];  */
	    masking_MS_ratio [gr_out] [chn].en  = gfc -> en  [chn+2];
	    masking_MS_ratio [gr_out] [chn].thm = gfc -> thm [chn+2];
	}
    }

    for (chn=0; chn<numchn; chn++) {
	FLOAT (*wsamp_l)[BLKSIZE];
	FLOAT (*wsamp_s)[3][BLKSIZE_s];
	FLOAT en_subshort[12];
	FLOAT attack_intensity[12];
	int ns_uselongblock = 1;
	FLOAT attackThreshold;
	FLOAT max[CBANDS],avg[CBANDS];
	int ns_attacks[4] = {0};
	FLOAT fftenergy[HBLKSIZE];
	FLOAT fftenergy_s[3][HBLKSIZE_s];
	/* convolution   */
	FLOAT eb[CBANDS+1],eb2[CBANDS];
	FLOAT thr[CBANDS+1];


    /*This is the masking table:
      According to tonality, values are going from 0dB (TMN)
      to 9.3dB (NMT).
      After additive masking computation, 8dB are added, so
      final values are going from 8dB to 17.3dB
    */
    static const FLOAT tab[] = {
        1.0/*pow(10, -0)*/,
        0.79433/*pow(10, -0.1)*/,
        0.63096/*pow(10, -0.2)*/,
        0.63096/*pow(10, -0.2)*/,
        0.63096/*pow(10, -0.2)*/,
        0.63096/*pow(10, -0.2)*/,
        0.63096/*pow(10, -0.2)*/,
        0.25119/*pow(10, -0.6)*/,
	    0.11749/*pow(10, -0.93)*/
	};

    /*  rh 20040301: the following loops do access one off the limits
     *  so I increase  the array dimensions by one and initialize the
     *  accessed values to zero
     */
    assert( gfc->npart_s <= CBANDS );
    assert( gfc->npart_l <= CBANDS );
    eb [gfc->npart_s] = 0;
    thr[gfc->npart_s] = 0;
    eb [gfc->npart_l] = 0;
    thr[gfc->npart_l] = 0;
    
	/*************************************************************** 
	 * determine the block type (window type)
	 ***************************************************************/
	/* calculate energies of each sub-shortblocks */
	for (i=0; i<3; i++) {
	    en_subshort[i] = gfc->nsPsy.last_en_subshort[chn][i+6];
	    attack_intensity[i]
		= en_subshort[i] / gfc->nsPsy.last_en_subshort[chn][i+4];
	}

	if (chn == 2) {
	    for(i=0;i<576;i++) {
		FLOAT l, r;
		l = ns_hpfsmpl[0][i];
		r = ns_hpfsmpl[1][i];
		ns_hpfsmpl[0][i] = l+r;
		ns_hpfsmpl[1][i] = l-r;
	    }
	}
	{
	    FLOAT *pf = ns_hpfsmpl[chn & 1];
	    for (i=0;i<9;i++) {
		FLOAT *pfe = pf + 576/9, p = 1.;
		for (; pf < pfe; pf++)
		    if (p < fabs(*pf))
			p = fabs(*pf);

		gfc->nsPsy.last_en_subshort[chn][i] = en_subshort[i+3] = p;
		if (p > en_subshort[i+3-2])
		    p = p / en_subshort[i+3-2];
		else if (en_subshort[i+3-2] > p*10.0)
		    p = en_subshort[i+3-2] / (p*10.0);
		else
		    p = 0.0;
		attack_intensity[i+3] = p;
	    }
	}
#if defined(HAVE_GTK)
	if (gfp->analysis) {
	    FLOAT x = attack_intensity[0];
	    for (i=1;i<12;i++) 
		if (x < attack_intensity[i])
		    x = attack_intensity[i];
	    gfc->pinfo->ers[gr_out][chn] = gfc->ers_save[chn];
	    gfc->ers_save[chn] = x;
	}
#endif
	/* compare energies between sub-shortblocks */
	attackThreshold = (chn == 3)
	    ? gfc->nsPsy.attackthre_s : gfc->nsPsy.attackthre;
	for (i=0;i<12;i++) 
	    if (!ns_attacks[i/3] && attack_intensity[i] > attackThreshold)
		ns_attacks[i/3] = (i % 3)+1;

	if (ns_attacks[0] && gfc->nsPsy.last_attacks[chn])
	    ns_attacks[0] = 0;

	if (gfc->nsPsy.last_attacks[chn] == 3 ||
	    ns_attacks[0] + ns_attacks[1] + ns_attacks[2] + ns_attacks[3]) {
	    ns_uselongblock = 0;

	    if (ns_attacks[1] && ns_attacks[0]) ns_attacks[1] = 0;
	    if (ns_attacks[2] && ns_attacks[1]) ns_attacks[2] = 0;
	    if (ns_attacks[3] && ns_attacks[2]) ns_attacks[3] = 0;
	}

    if (chn < 2) {
        uselongblock[chn] = ns_uselongblock;
    }
    else {
        if (ns_uselongblock == 0) {
            uselongblock[0] = uselongblock[1] = 0;
        }
    }

	/* there is a one granule delay.  Copy maskings computed last call
	 * into masking_ratio to return to calling program.
	 */
	energy[chn]=gfc->tot_ener[chn];

	/*********************************************************************
	 *  compute FFTs
	 *********************************************************************/
	wsamp_s = wsamp_S+(chn & 1);
	wsamp_l = wsamp_L+(chn & 1);
	compute_ffts(gfp, fftenergy, fftenergy_s,
		     wsamp_l, wsamp_s, gr_out, chn, buffer);

	/* compute masking thresholds for short blocks */
	for (sblock = 0; sblock < 3; sblock++) {
	    FLOAT enn, thmm;
	    compute_masking_s(gfc, fftenergy_s, eb, thr, chn, sblock,
			      gfp->ATHlower*gfc->ATH->adjust);
        b = -1;
	    for (sb = 0; sb < SBMAX_s; sb++) {
		enn = thmm = 0.0;
		while (++b < gfc->bo_s[sb]) {
		    enn  += eb[b];
		    thmm += thr[b];
		}
		enn  += 0.5 * eb[b];    /* for the last sfb b is larger than npart_s!! */
		thmm += 0.5 * thr[b];   /* rh 20040301 */
		gfc->en [chn].s[sb][sblock] = enn;
        
        assert( enn >= 0 );
        assert( thmm >= 0 );
        
		/****   short block pre-echo control   ****/
		thmm *= NS_PREECHO_ATT0;
		if (ns_attacks[sblock] >= 2 || ns_attacks[sblock+1] == 1) {
		    int idx = (sblock != 0) ? sblock-1 : 2;
		    double p = NS_INTERP(gfc->thm[chn].s[sb][idx],
					 thmm, NS_PREECHO_ATT1*pcfact);
		    thmm = Min(thmm,p);
		}

		if (ns_attacks[sblock] == 1) {
		    int idx = (sblock != 0) ? sblock-1 : 2;
		    double p = NS_INTERP(gfc->thm[chn].s[sb][idx],
					 thmm,NS_PREECHO_ATT2*pcfact);
		    thmm = Min(thmm,p);
		} else if ((sblock != 0 && ns_attacks[sblock-1] == 3)
			|| (sblock == 0 && gfc->nsPsy.last_attacks[chn] == 3)) {
		    int idx = (sblock != 2) ? sblock+1 : 0;
		    double p = NS_INTERP(gfc->thm[chn].s[sb][idx],
					 thmm,NS_PREECHO_ATT2*pcfact);
		    thmm = Min(thmm,p);
		}

		/* pulse like signal detection for fatboy.wav and so on */
		enn = en_subshort[sblock*3+3] + en_subshort[sblock*3+4]
		    + en_subshort[sblock*3+5];
		if (en_subshort[sblock*3+5]*6 < enn) {
		    thmm *= 0.5;
		    if (en_subshort[sblock*3+4]*6 < enn)
			thmm *= 0.5;
		}

		gfc->thm[chn].s[sb][sblock] = thmm;
	    }