psych.c

MPEG-4编解码的实现(包括MPEG4视音频编解码)

	double dataL[2048], dataR[2048];

	for (i = 0; i < psyInfoL->size*2; i++) {
		a = psyInfoL->transBuff[i];
		b = psyInfoR->transBuff[i];
		dataL[i] = (a+b)*SQRT2*0.5;
		dataR[i] = (a-b)*SQRT2*0.5;
	}

	/* Calculate magnitude of new data */
	for (i = 0; i < psyInfoL->size; i++) {
		a = dataL[i];
		b = dataL[i+psyInfoL->size];
		psyInfoL->energyMS[i] = 0.5 * (a*a + b*b);

		a = dataR[i];
		b = dataR[i+psyInfoL->size];
		psyInfoR->energyMS[i] = 0.5 * (a*a + b*b);
	}

	for (j = 0; j < 8; j++) {

		for (i = 0; i < psyInfoL->sizeS*2; i++) {
			a = psyInfoL->transBuffS[j][i];
			b = psyInfoR->transBuffS[j][i];
			dataL[i] = (a+b)*SQRT2*0.5;
			dataR[i] = (a-b)*SQRT2*0.5;
		}

		/* Calculate magnitude of new data */
		for (i = 0; i < psyInfoL->sizeS; i++) {
			a = dataL[i];
			b = dataL[i+psyInfoL->sizeS];
			psyInfoL->energySMS[j][i] = 0.5 * (a*a + b*b);

			a = dataR[i];
			b = dataR[i+psyInfoL->sizeS];
			psyInfoR->energySMS[j][i] = 0.5 * (a*a + b*b);
		}
	}
}

/* addition of simultaneous masking */
__inline double mask_add(double m1, double m2, int k, int b, double *ath)
{
	static const double 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 double 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
	};

	static const double 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;
	double m;

	if (m1 == 0) return m2;

	if (b < 0) b = -b;

	i = (int)(10*log10(m2 / m1)/10*16);
	m = 10*log10((m1+m2)/ath[k]);

	if (i < 0) i = -i;

	if (b <= 3) { /* approximately, 1 bark = 3 partitions */
		if (i > 8) return m1+m2;
		return (m1+m2)*table2[i];
	}

	if (m<15) {
		if (m > 0) {
			double f=1.0,r;
			if (i > 24) return m1+m2;
			if (i > 13) f = 1; else f = table3[i];
			r = (m-0)/15;
			return (m1+m2)*(table1[i]*r+f*(1-r));
		}
		if (i > 13) return m1+m2;
		return (m1+m2)*table3[i];
	}

	if (i > 24) return m1+m2;
	return (m1+m2)*table1[i];
}

static void PsyThreshold(GlobalPsyInfo *gpsyInfo, PsyInfo *psyInfo, int *cb_width_long,
						 int num_cb_long, int *cb_width_short, int num_cb_short)
{
	int b, bb, w, low, high, j;
	double tmp, ecb;

	double e[MAX_NPART];
	double c[MAX_NPART];
	double maxi[MAX_NPART];
	double avg[MAX_NPART];
	double eb;

	double nb_tmp[1024], epart, npart;

	double tot, mx, estot[8];
	double pe = 0.0;

	/* Energy in each partition and weighted unpredictability */
	high = 0;
	for (b = 0; b < gpsyInfo->psyPart->len; b++)
	{
		double m, a;
		low = high;
		high += gpsyInfo->psyPart->width[b];

		eb = psyInfo->energy[low];
		m = a = eb;

		for (w = low+1; w < high; w++)
		{
			double el = psyInfo->energy[w];
			eb += el;
			a += el;
			m = m < el ? el : m;
		}
		e[b] = eb;
		maxi[b] = m;
		avg[b] = a / gpsyInfo->psyPart->width[b];
	}

	for (b = 0; b < gpsyInfo->psyPart->len; b++)
	{
		static double tab[20] = {
			1,0.79433,0.63096,0.63096,0.63096,0.63096,0.63096,0.25119,0.11749,0.11749,
			0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749
		};
		int c1,c2,t;
		double m, a, tonality;

		c1 = c2 = 0;
		m = a = 0;
		for(w = b-1; w <= b+1; w++)
		{
			if (w >= 0 && w < gpsyInfo->psyPart->len) {
				c1++;
				c2 += gpsyInfo->psyPart->width[w];
				a += avg[w];
				m = m < maxi[w] ? maxi[w] : m;
			}
		}

		a /= c1;
		tonality = (a == 0) ? 0 : (m / a - 1)/(c2-1);

		t = (int)(20*tonality);
		if (t > 19) t = 19;
		psyInfo->tonality[b] = tab[t];
		c[b] = e[b] * tab[t];
	}

	/* Convolve the partitioned energy and unpredictability
	   with the spreading function */
	for (b = 0; b < gpsyInfo->psyPart->len; b++)
	{
		ecb = 0;
		for (bb = gpsyInfo->sprInd[b][0]; bb < gpsyInfo->sprInd[b][1]; bb++)
		{
			ecb = mask_add(ecb, gpsyInfo->spreading[b][bb] * c[bb], bb, bb-b, gpsyInfo->ath);
		}
		ecb *= 0.158489319246111;

		/* Actual energy threshold */
		psyInfo->nb[b] = NS_INTERP(min(ecb, 2*psyInfo->lastNb[b]), ecb, 1/*pcfact*/);
/*
		psyInfo->nb[b] = max(psyInfo->nb[b], gpsyInfo->ath[b]);
*/
		psyInfo->lastNb[b] = ecb;

		/* Perceptual entropy */
		tmp = gpsyInfo->psyPart->width[b]
			* log((psyInfo->nb[b] + 0.0000000001)
			/ (e[b] + 0.0000000001));
		tmp = min(0,tmp);

		pe -= tmp;
	}

	high = 0;
	for (b = 0; b < gpsyInfo->psyPart->len; b++)
	{
		low = high;
		high += gpsyInfo->psyPart->width[b];

		for (w = low; w < high; w++)
		{
			nb_tmp[w] = psyInfo->nb[b] / gpsyInfo->psyPart->width[b];
		}
	}

	high = 0;
	for (b = 0; b < num_cb_long; b++)
	{
		low = high;
		high += cb_width_long[b];

		epart = psyInfo->energy[low];
		npart = nb_tmp[low];
		for (w = low+1; w < high; w++)
		{
			epart += psyInfo->energy[w];

			if (nb_tmp[w] < npart)
				npart = nb_tmp[w];
		}
		npart *= cb_width_long[b];

		psyInfo->maskThr[b] = psyInfo->maskThrNext[b];
		psyInfo->maskEn[b] = psyInfo->maskEnNext[b];
		tmp = npart / epart;
		psyInfo->maskThrNext[b] = npart;
		psyInfo->maskEnNext[b] = epart;
	}

	/* Short windows */
	for (j = 0; j < 8; j++)
	{
		/* Energy in each partition and weighted unpredictability */
		high = 0;
		for (b = 0; b < gpsyInfo->psyPartS->len; b++)
		{
			low = high;
			high += gpsyInfo->psyPartS->width[b];

			eb = psyInfo->energyS[j][low];

			for (w = low+1; w < high; w++)
			{
				double el = psyInfo->energyS[j][w];
				eb += el;
			}
			e[b] = eb;
		}

		estot[j] = 0.0;

		/* Convolve the partitioned energy and unpredictability
		with the spreading function */
		for (b = 0; b < gpsyInfo->psyPartS->len; b++)
		{
			ecb = 0;
			for (bb = gpsyInfo->sprIndS[b][0]; bb <= gpsyInfo->sprIndS[b][1]; bb++)
			{
				ecb += gpsyInfo->spreadingS[b][bb] * e[bb];
			}

			/* Actual energy threshold */
			psyInfo->nbS[j][b] = max(1e-6, ecb);
/*
			psyInfo->nbS[j][b] = max(psyInfo->nbS[j][b], gpsyInfo->athS[b]);
*/

			estot[j] += e[b];
		}

		if (estot[j] != 0.0)
			estot[j] /= gpsyInfo->psyPartS->len;

		high = 0;
		for (b = 0; b < gpsyInfo->psyPartS->len; b++)
		{
			low = high;
			high += gpsyInfo->psyPartS->width[b];

			for (w = low; w < high; w++)
			{
				nb_tmp[w] = psyInfo->nbS[j][b] / gpsyInfo->psyPartS->width[b];
			}
		}

		high = 0;
		for (b = 0; b < num_cb_short; b++)
		{
			low = high;
			high += cb_width_short[b];

			epart = psyInfo->energyS[j][low];
			npart = nb_tmp[low];
			for (w = low+1; w < high; w++)
			{
				epart += psyInfo->energyS[j][w];

				if (nb_tmp[w] < npart)
					npart = nb_tmp[w];
			}
			npart *= cb_width_short[b];

			psyInfo->maskThrS[j][b] = psyInfo->maskThrNextS[j][b];
			psyInfo->maskEnS[j][b] = psyInfo->maskEnNextS[j][b];
			psyInfo->maskThrNextS[j][b] = npart;
			psyInfo->maskEnNextS[j][b] = epart;
		}
	}

	tot = mx = estot[0];
	for (j = 1; j < 8; j++) {
		tot += estot[j];
		mx = max(mx, estot[j]);
	}

#ifdef _DEBUG
	printf("%4f %2.2f ", pe, mx/tot);
#endif

	tot = max(tot, 1.e-12);
	if (((mx/tot) > 0.35) && (pe > 1800.0) || ((mx/tot) > 0.5) || (pe > 3000.0)) {
		psyInfo->block_type = ONLY_SHORT_WINDOW;
		psyInfo->threeInARow++;
	} else if ((psyInfo->lastEnr > 0.5) || (psyInfo->lastPe > 3000.0)) {
		psyInfo->block_type = ONLY_SHORT_WINDOW;
		psyInfo->threeInARow++;
	} else if (psyInfo->threeInARow >= 3) {
		psyInfo->block_type = ONLY_SHORT_WINDOW;
		psyInfo->threeInARow = 0;
	} else {
		psyInfo->block_type = ONLY_LONG_WINDOW;
	}

 	psyInfo->lastEnr = mx/tot;
	psyInfo->pe = psyInfo->lastPe;
	psyInfo->lastPe = pe;
}

static void PsyThresholdMS(ChannelInfo *channelInfoL, GlobalPsyInfo *gpsyInfo,
						   PsyInfo *psyInfoL, PsyInfo *psyInfoR,
						   int *cb_width_long, int num_cb_long, int *cb_width_short,
						   int num_cb_short)
{
	int b, bb, w, low, high, j;
	double ecb, tmp1, tmp2;

	double nb_tmpM[1024];
	double nb_tmpS[1024];
	double epartM, epartS, npartM, npartS;

	double nbM[MAX_NPART];
	double nbS[MAX_NPART];
	double eM[MAX_NPART];
	double eS[MAX_NPART];
	double cM[MAX_NPART];
	double cS[MAX_NPART];

	double mld;

#ifdef _DEBUG
	int ms_used = 0;
	int ms_usedS = 0;
#endif

	/* Energy in each partition and weighted unpredictability */
	high = 0;
	for (b = 0; b < gpsyInfo->psyPart->len; b++)
	{
		double mid, side, ebM, ebS;
		low = high;
		high += gpsyInfo->psyPart->width[b];

		mid  = psyInfoL->energyMS[low];
		side = psyInfoR->energyMS[low];

		ebM = mid;
		ebS = side;

		for (w = low+1; w < high; w++)
		{
			mid  = psyInfoL->energyMS[w];
			side = psyInfoR->energyMS[w];

			ebM += mid;
			ebS += side;
		}
		eM[b] = ebM;
		eS[b] = ebS;
		cM[b] = ebM * min(psyInfoL->tonality[b], psyInfoR->tonality[b]);
		cS[b] = ebS * min(psyInfoL->tonality[b], psyInfoR->tonality[b]);
	}

	/* Convolve the partitioned energy and unpredictability
	   with the spreading function */
	for (b = 0; b < gpsyInfo->psyPart->len; b++)
	{
		/* Mid channel */

		ecb = 0;
		for (bb = gpsyInfo->sprInd[b][0]; bb <= gpsyInfo->sprInd[b][1]; bb++)
		{
			ecb = mask_add(ecb, gpsyInfo->spreading[bb][b] * cM[bb], bb, bb-b, gpsyInfo->ath);
		}
		ecb *= 0.158489319246111;

		/* Actual energy threshold */
		nbM[b] = NS_INTERP(min(ecb, 2*psyInfoL->lastNbMS[b]), ecb, 1/*pcfact*/);
/*
		nbM[b] = max(nbM[b], gpsyInfo->ath[b]);
*/
		psyInfoL->lastNbMS[b] = ecb;


		/* Side channel */

		ecb = 0;
		for (bb = gpsyInfo->sprInd[b][0]; bb <= gpsyInfo->sprInd[b][1]; bb++)
		{
			ecb = mask_add(ecb, gpsyInfo->spreading[bb][b] * cS[bb], bb, bb-b, gpsyInfo->ath);
		}
		ecb *= 0.158489319246111;

		/* Actual energy threshold */
		nbS[b] = NS_INTERP(min(ecb, 2*psyInfoR->lastNbMS[b]), ecb, 1/*pcfact*/);
/*
		nbS[b] = max(nbS[b], gpsyInfo->ath[b]);
*/
		psyInfoR->lastNbMS[b] = ecb;

		if (psyInfoL->nb[b] <= 1.58*psyInfoR->nb[b]
			&& psyInfoR->nb[b] <= 1.58*psyInfoL->nb[b]) {

			mld = gpsyInfo->mld[b]*eM[b];
			tmp1 = max(nbM[b], min(nbS[b],mld));

			mld = gpsyInfo->mld[b]*eS[b];
			tmp2 = max(nbS[b], min(nbM[b],mld));

			nbM[b] = tmp1;
			nbS[b] = tmp2;
		}
	}

	high = 0;
	for (b = 0; b < gpsyInfo->psyPart->len; b++)
	{
		low = high;
		high += gpsyInfo->psyPart->width[b];

		for (w = low; w < high; w++)
		{
			nb_tmpM[w] = nbM[b] / gpsyInfo->psyPart->width[b];
			nb_tmpS[w] = nbS[b] / gpsyInfo->psyPart->width[b];
		}
	}

	high = 0;
	for (b = 0; b < num_cb_long; b++)
	{
		low = high;
		high += cb_width_long[b];