psymodel.c

音频编码

	while (--kk > 0)
	    ecb += fftenergy_s[sblock][j++];
	eb[b] = ecb;
    }
    for (j = b = 0; b < gfc->npart_s; b++) {
	int kk = gfc->s3ind_s[b][0];
	FLOAT ecb = gfc->s3_ss[j++] * eb[kk++];
	while (kk <= gfc->s3ind_s[b][1])
	    ecb += gfc->s3_ss[j++] * eb[kk++];

	thr[b] = Min( ecb, rpelev_s  * gfc->nb_s1[chn][b] );
	if (gfc->blocktype_old[chn & 1] == SHORT_TYPE ) {
	    thr[b] = Min(thr[b], rpelev2_s * gfc->nb_s2[chn][b]);
	}
	thr[b] = Max( thr[b],
                  Min(gfc->ATH->cb[gfc->bm_s[b]] * athlower,
                      thr[b] * 2) );
	gfc->nb_s2[chn][b] = gfc->nb_s1[chn][b];
	gfc->nb_s1[chn][b] = ecb;
    assert( thr[b] >= 0 );
    }
}

static void
block_type_set(
    lame_global_flags * gfp,
    int *uselongblock,
    int *blocktype_d,
    int *blocktype
    )
{
    lame_internal_flags *gfc=gfp->internal_flags;
    int chn;

    if (gfp->short_blocks == short_block_coupled
	/* force both channels to use the same block type */
	/* this is necessary if the frame is to be encoded in ms_stereo.  */
	/* But even without ms_stereo, FhG  does this */
	&& !(uselongblock[0] && uselongblock[1]))
        uselongblock[0] = uselongblock[1] = 0;

    /* update the blocktype of the previous granule, since it depends on what
     * happend in this granule */
    for (chn=0; chn<gfc->channels_out; chn++) {
	blocktype[chn] = NORM_TYPE;
	/* disable short blocks */
	if (gfp->short_blocks == short_block_dispensed)
	    uselongblock[chn]=1;
	if (gfp->short_blocks == short_block_forced)
	    uselongblock[chn]=0;

	if (uselongblock[chn]) {
	    /* no attack : use long blocks */
	    assert( gfc->blocktype_old[chn] != START_TYPE );
	    if (gfc->blocktype_old[chn] == SHORT_TYPE)
		blocktype[chn] = STOP_TYPE;
	} else {
	    /* attack : use short blocks */
	    blocktype[chn] = SHORT_TYPE;
	    if (gfc->blocktype_old[chn] == NORM_TYPE) {
		gfc->blocktype_old[chn] = START_TYPE;
	    }
	    if (gfc->blocktype_old[chn] == STOP_TYPE)
		gfc->blocktype_old[chn] = SHORT_TYPE;
	}

	blocktype_d[chn] = gfc->blocktype_old[chn];  /* value returned to calling program */
	gfc->blocktype_old[chn] = blocktype[chn];    /* save for next call to l3psy_anal */
    }
}

static void 
determine_block_type( lame_global_flags * gfp, FLOAT fftenergy_s[3][HBLKSIZE_s], int uselongblock[],
    int chn, int gr_out, FLOAT* pe )
{
    lame_internal_flags* gfc = gfp->internal_flags;
    int j;
	/*************************************************************** 
	 * determine the block type (window type) based on L & R channels
	 ***************************************************************/
	/* compute PE for all 4 channels */
	    FLOAT mn,mx,ma=0,mb=0,mc=0;
            for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++) {
		ma += fftenergy_s[0][j];
		mb += fftenergy_s[1][j];
		mc += fftenergy_s[2][j];
	    }
	    mn = Min(ma,mb);
	    mn = Min(mn,mc);
	    mx = Max(ma,mb);
	    mx = Max(mx,mc);
#if defined(HAVE_GTK)
	    if (gfp->analysis) {
		gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn];
		gfc->ers_save[chn]=(mx/(1e-12+mn));
	    }
#endif
	    /* bit allocation is based on pe.  */
	    if (mx>mn) {
		FLOAT tmp = FAST_LOG_X(mx/(1e-12+mn), 400.0);
		if (tmp > *pe) *pe = tmp;
	    }

	    /* block type is based just on L or R channel */      
	    if (chn<2) {
		uselongblock[chn] = 1;

		/* tuned for t1.wav.  doesnt effect most other samples */
		if (*pe > 3000) 
		    uselongblock[chn]=0;
	
		if ( mx > 30*mn ) 
		{/* big surge of energy - always use short blocks */
		    uselongblock[chn] = 0;
		}
		else if ((mx > 10*mn) && (*pe > 1000))
		{/* medium surge, medium pe - use short blocks */
		    uselongblock[chn] = 0;
		}
	    }
}

int L3psycho_anal( 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])
{
    lame_internal_flags *gfc=gfp->internal_flags;

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

    /* convolution   */
    FLOAT eb[CBANDS+1];
    FLOAT cb[CBANDS];
    FLOAT thr[CBANDS+1];

    /* ratios    */
    FLOAT ms_ratio_l=0, ms_ratio_s=0;

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

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

    /*  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;

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

    for (chn=0; chn<numchn; chn++) {
	FLOAT (*wsamp_l)[BLKSIZE];
	FLOAT (*wsamp_s)[3][BLKSIZE_s];
	energy[chn] = gfc->tot_ener[chn];

	/* there is a one granule delay.  Copy maskings computed last call
	 * into masking_ratio to return to calling program.
	 */
	if (chn < 2) {
	    /* LR maskings  */
	    percep_entropy            [chn]       = gfc -> pe  [chn];
	    masking_ratio    [gr_out] [chn]  .en  = gfc -> en  [chn];
	    masking_ratio    [gr_out] [chn]  .thm = gfc -> thm [chn];
	} else {
	    /* MS maskings  */
	    percep_MS_entropy         [chn-2]     = gfc -> pe  [chn]; 
	    masking_MS_ratio [gr_out] [chn-2].en  = gfc -> en  [chn];
	    masking_MS_ratio [gr_out] [chn-2].thm = gfc -> thm [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 unpredicatability of first six spectral lines
	 *********************************************************************/
	for ( j = 0; j < CW_LOWER_INDEX; j++ ) {
	    /* calculate unpredictability measure cw */
	    FLOAT a2, b2, r1, r2;
	    FLOAT numre, numim, den;

	    a2 = gfc-> ax_sav[chn][1][j];
	    b2 = gfc-> bx_sav[chn][1][j];

	    r2 = gfc-> rx_sav[chn][1][j];
	    r1 = gfc-> rx_sav[chn][1][j] = gfc-> rx_sav[chn][0][j];

	    /* square (x1,y1) */
	    if (r1 != 0.0) {
		FLOAT a1 = gfc-> ax_sav[chn][1][j] = gfc-> ax_sav[chn][0][j];
		FLOAT b1 = gfc-> bx_sav[chn][1][j] = gfc-> bx_sav[chn][0][j];
		den = r1*r1;
		numre = a1*b1;
		numim = den-b1*b1;
	    } else {
		/* no aging is needed for ax_sav[chn][0][j] and that of bx
		   because if r1=0, r2 should be 0 for next time. */
		den = numre = 1.0;
		numim = 0.0;
	    }

	    /* multiply by (x2,-y2) */
	    if (r2 != 0.0) {
		FLOAT tmp2 = (numim+numre)*(a2+b2)*0.5f;
		FLOAT tmp1 = -a2*numre+tmp2;
		numre =      -b2*numim+tmp2;
		numim = tmp1;
		den *= r2;
	    }

	    r1 = 2.0f*r1-r2;
	    r2 = gfc-> rx_sav[chn][0][j] = sqrt(fftenergy[j]);
	    r2 = r2+fabs(r1);
	    if (r2 != 0) {
		FLOAT an = gfc-> ax_sav[chn][0][j] = wsamp_l[0][j];
		FLOAT bn = gfc-> bx_sav[chn][0][j] = j==0 ? wsamp_l[0][0] : wsamp_l[0][BLKSIZE-j];  
		den = r1/den*2.0;
		numre = (an+bn)-numre*den;
		numim = (an-bn)-numim*den;
		r2 = sqrt(numre*numre+numim*numim)/(r2*2.0);
	    }
	    gfc->cw[j] = r2;
	}

	/**********************************************************************
	 *     compute unpredicatibility of next 200 spectral lines
	 *********************************************************************/
	for (; j < gfc->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 = fftenergy_s[0][k];
	    if (r1 != 0.0) {
		FLOAT a1 = (*wsamp_s)[0][k]; 
		FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */
		numre = a1*b1;
		numim = r1-b1*b1;
		den = r1;
		r1 = sqrt(r1);
	    } else {
		den = numre = 1.0;
		numim = 0.0;
	    }

	    /* multiply by (x2,-y2) */
	    r2 = fftenergy_s[2][k];
	    if (r2 != 0.0) {
		FLOAT a2 = (*wsamp_s)[2][k]; 
		FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k];

		FLOAT tmp2 = (numim+numre)*(a2+b2)*0.5f;
		FLOAT tmp1 = tmp2-a2*numre;
		numre =      tmp2-b2*numim;
		numim = tmp1;

		r2 = sqrt(r2);
		den *= r2;
	    }

	    /* r-prime factor */
	    rn = sqrt(fftenergy_s[1][k])+fabs(2*r1-r2);
	    if (rn != 0) {
		FLOAT an = (*wsamp_s)[1][k]; 
		FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k];
		den = (2*r1-r2)/den*2.0f;
		numre = (an+bn)-numre*den;
		numim = (an-bn)-numim*den;
		rn = sqrt(numre*numre+numim*numim)/(rn*2.0f);
	    }
	    gfc->cw[j+1] = gfc->cw[j+2] = gfc->cw[j+3] = gfc->cw[j] = rn;
	}

	/**********************************************************************
	 *    Calculate the energy and the unpredictability in the threshold
	 *    calculation partitions
	 *********************************************************************/
	b = 0;
	for (j = 0; j < gfc->cw_upper_index
		 && gfc->numlines_l[b] && b < gfc->npart_l; ) {
	    FLOAT ebb, cbb;

	    ebb = NON_LINEAR_SCALE_ITEM(fftenergy[j]);
	    cbb = NON_LINEAR_SCALE_ITEM(fftenergy[j] * gfc->cw[j]);
	    j++;

	    for (i = gfc->numlines_l[b] - 1; i > 0; i--) {
		ebb += NON_LINEAR_SCALE_ITEM(fftenergy[j]);
		/* XXX: should "* gfc->cw[j])" be outside of the scaling? */
		cbb += NON_LINEAR_SCALE_ITEM(fftenergy[j] * gfc->cw[j]);
		j++;
	    }
	    eb[b] = NON_LINEAR_SCALE_SUM(ebb);
	    cb[b] = NON_LINEAR_SCALE_SUM(cbb);
	    b++;
	}

	for (; b < gfc->npart_l; b++ ) {
	    FLOAT ebb = NON_LINEAR_SCALE_ITEM(fftenergy[j++]);
	    assert(gfc->numlines_l[b]);
	    for (i = gfc->numlines_l[b] - 1; i > 0; i--) {
		ebb += NON_LINEAR_SCALE_ITEM(fftenergy[j++]);
	    }
	    eb[b] = NON_LINEAR_SCALE_SUM(ebb);
	    /* XXX: should the "* .4" be outside of the scaling? */
	    cb[b] = NON_LINEAR_SCALE_SUM(ebb * 0.4);
	}

	/**********************************************************************
	 *      convolve the partitioned energy and unpredictability
	 *      with the spreading function, s3_l[b][k](packed into s3_ll)
	 *********************************************************************/
	/*  calculate percetual entropy */
	gfc->pe[chn] = 0;
	k = 0;
	for ( b = 0;b < gfc->npart_l; b++ ) {
	    FLOAT tbb,ecb,ctb;
	    int kk;
	    ecb = ctb = 0.;
	    for (kk = gfc->s3ind[b][0]; kk <= gfc->s3ind[b][1]; kk++ ) {
		/* sprdngf for Layer III */
		ecb += gfc->s3_ll[k] * eb[kk];
		ctb += gfc->s3_ll[k] * cb[kk];
		k++;
	    }

/* calculate the tonality of each threshold calculation partition 
 * calculate the SNR in each threshold calculation partition 
 * tonality = -0.299 - .43*log(ctb/ecb);
 * tonality = 0:           use NMT   (lots of masking)
 * tonality = 1:           use TMN   (little masking)
 */
	    tbb = ecb;
	    if (tbb != 0.0) {
		tbb = ctb / tbb;
		/* convert to tonality index */
		/* tonality small:   tbb=1 */
		/* tonality large:   tbb=-.299 */
		tbb = CONV1 + FAST_LOG_X(tbb, CONV2);
		if (tbb < 0.0) tbb = exp(-LN_TO_LOG10*NMT);
		else if (tbb > 1.0) tbb = exp(-LN_TO_LOG10*TMN);
		else tbb = exp(-LN_TO_LOG10 * ( (TMN-NMT)*tbb + NMT ));
	    }

/* at this point, tbb represents the amount the spreading function
 * will be reduced.  The smaller the value, the less masking.
 * minval[] = 1 (0db)     says just use tbb.
 * minval[]= .01 (-20db)  says reduce spreading function by at least 20db.
 */
	    tbb = Min(gfc->minval[b], tbb);
	    /* stabilize tonality estimation */
	    if (gfc->PSY->tonalityPatch && b > 5) {
		FLOAT const x = 1.8699422;
		FLOAT w = gfc->PSY->prvTonRed[b/2] * x;
		if (tbb > w) 
		    tbb = w;
		gfc->PSY->prvTonRed[b] = tbb;
	    }
	    ecb *= tbb;

	    /* long block pre-echo control.   */
	    /* rpelev=2.0, rpelev2=16.0 */
	    /* note: all surges in PE are because of this pre-echo formula
	     * for thr[b].  If it this is not used, PE is always around 600
	     */
	    /* dont use long block pre-echo control if previous granule was
	     * a short block.  This is to avoid the situation:   
	     * frame0:  quiet (very low masking)  
	     * frame1:  surge  (triggers short blocks)
	     * frame2:  regular frame. looks like pre-echo when compared to
	     *          frame0, but all pre-echo was in frame1.
	     */
	    /* chn=0,1   L and R channels
	       chn=2,3   S and M channels.  
	    */
	    thr[b] = Min(ecb, rpelev*gfc->nb_1[chn][b]);
	    if (gfc->blocktype_old[chn & 1] != SHORT_TYPE
		&& thr[b] > rpelev2*gfc->nb_2[chn][b])
		thr[b] = rpelev2*gfc->nb_2[chn][b];

	    gfc->nb_2[chn][b] = gfc->nb_1[chn][b];
	    gfc->nb_1[chn][b] = ecb;

	    ecb = Max(thr[b], gfc->ATH->cb[b]*gfc->ATH->adjust);
	    if (ecb < eb[b])
		gfc->pe[chn] -= gfc->numlines_l[b] * FAST_LOG(ecb / eb[b]);
	}

    determine_block_type( gfp, fftenergy_s, uselongblock, chn, gr_out, &gfc->pe[chn] );

	/* compute masking thresholds for long blocks */
	convert_partition2scalefac_l(gfc, eb, thr, chn);

	/* 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;
	    enn = thmm = 0.0;
	    for (sb = 0; sb < SBMAX_s; sb++) {
		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 */
        assert( enn >= 0 );
        assert( thmm >= 0 );
		gfc->en [chn].s[sb][sblock] = enn;
		gfc->thm[chn].s[sb][sblock] = thmm;
		enn  = 0.5 * eb[b];
		thmm = 0.5 * thr[b];
        assert( enn >= 0 );
        assert( thmm >= 0 );
	    }
	    gfc->en [chn].s[sb-1][sblock] += enn;
	    gfc->thm[chn].s[sb-1][sblock] += thmm;
	}
    } /* end loop over chn */

    if (gfp->interChRatio != 0.0)
	calc_interchannel_masking(gfp, gfp->interChRatio);

    if (gfp->mode == JOINT_STEREO) {
	FLOAT db,x1,x2,sidetot=0,tot=0;
	msfix1(gfc);
	if (gfp->msfix != 0.0)
	    ns_msfix(gfc, gfp->msfix, gfp->ATHlower*gfc->ATH->adjust);

	/* determin ms_ratio from masking thresholds*/
	/* use ms_stereo (ms_ratio < .35) if average thresh. diff < 5 db */
	for (sb= SBMAX_l/4 ; sb< SBMAX_l; sb ++ ) {
	    x1 = Min(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
	    x2 = Max(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
	    /* thresholds difference in db */
	    if (x2 >= 1000*x1)  db=3;
	    else db = FAST_LOG10(x2/x1);  
	    /*  DEBUGF(gfc,"db = %f %e %e  \n",db,gfc->thm[0].l[sb],gfc->thm[1].l[sb]);*/
	    sidetot += db;
	    tot++;
	}
	ms_ratio_l= (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
	ms_ratio_l = Min(ms_ratio_l,0.5);

	sidetot=0; tot=0;
	for ( sblock = 0; sblock < 3; sblock++ )
	    for ( sb = SBMAX_s/4; sb < SBMAX_s; sb++ ) {
		x1 = Min(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
		x2 = Max(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
		/* thresholds difference in db */
		if (x2 >= 1000*x1)  db=3;
		else db = FAST_LOG10(x2/x1);
		sidetot += db;
		tot++;
	    }
	ms_ratio_s = (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
	ms_ratio_s = Min(ms_ratio_s,.5);
    }

    /*************************************************************** 
     * determine final block type
     ***************************************************************/
    block_type_set(gfp, uselongblock, blocktype_d, blocktype);

    if (blocktype_d[0]==SHORT_TYPE && blocktype_d[1]==SHORT_TYPE)
	*ms_ratio = gfc->ms_ratio_s_old;
    else
	*ms_ratio = gfc->ms_ratio_l_old;

    gfc->ms_ratio_s_old = ms_ratio_s;
    gfc->ms_ratio_l_old = ms_ratio_l;

    /* we dont know the block type of this frame yet - assume long */
    *ms_ratio_next = ms_ratio_l;

    return 0;
}



/* mask_add optimization */
/* init the limit values used to avoid computing log in mask_add when it is not necessary */

/* For example, with i = 10*log10(m2/m1)/10*16         (= log10(m2/m1)*16)
 *
 * abs(i)>8 is equivalent (as i is an integer) to
 * abs(i)>=9
 * i>=9 || i<=-9
 * equivalent to (as i is the biggest integer smaller than log10(m2/m1)*16 
 * or the smallest integer bigger than log10(m2/m1)*16 depending on the sign of log10(m2/m1)*16)
 * log10(m2/m1)>=9/16 || log10(m2/m1)<=-9/16
 * exp10 is strictly increasing thus this is equivalent to
 * m2/m1 >= 10^(9/16) || m2/m1<=10^(-9/16) which are comparisons to constants
 */