layer3.c

xmms-1.2.10.tar.gz学习使用的就下吧

	/* 31 alias-reduction operations between each pair of sub-bands */
	/* with 8 butterflies between each pair                         */

	{
		int sb;
		real *xr1 = (real *) xr[1];

		if (sblim < 1 || sblim > SBLIMIT)
			return;
		
		for (sb = sblim; sb; sb--, xr1 += 10)
		{
			int ss;
			real *cs = aa_cs, *ca = aa_ca;
			real *xr2 = xr1;

			for (ss = 7; ss >= 0; ss--)
			{	/* upper and lower butterfly inputs */
				register real bu = *--xr2, bd = *xr1;

				*xr2 = (bu * (*cs)) - (bd * (*ca));
				*xr1++ = (bd * (*cs++)) + (bu * (*ca++));
			}
		}
	}
}

/* 
 This is an optimized DCT from Jeff Tsay's maplay 1.2+ package.
 Saved one multiplication by doing the 'twiddle factor' stuff
 together with the window mul. (MH)

 This uses Byeong Gi Lee's Fast Cosine Transform algorithm, but the
 9 point IDCT needs to be reduced further. Unfortunately, I don't
 know how to do that, because 9 is not an even number. - Jeff.

 ****************************************************************

 9 Point Inverse Discrete Cosine Transform

 This piece of code is Copyright 1997 Mikko Tommila and is freely usable
 by anybody. The algorithm itself is of course in the public domain.

 Again derived heuristically from the 9-point WFTA.

 The algorithm is optimized (?) for speed, not for small rounding errors or
 good readability.

 36 additions, 11 multiplications

 Again this is very likely sub-optimal.

 The code is optimized to use a minimum number of temporary variables,
 so it should compile quite well even on 8-register Intel x86 processors.
 This makes the code quite obfuscated and very difficult to understand.

 References:
 [1] S. Winograd: "On Computing the Discrete Fourier Transform",
     Mathematics of Computation, Volume 32, Number 141, January 1978,
     Pages 175-199
*/

/*------------------------------------------------------------------*/
/*                                                                  */
/*    Function: Calculation of the inverse MDCT                     */
/*                                                                  */
/*------------------------------------------------------------------*/

#define MACRO(v)						\
do {								\
	real tmpval;						\
								\
	tmpval = tmp[(v)] + tmp[17-(v)];			\
	out2[9+(v)] = tmpval * w[27+(v)];			\
	out2[8-(v)] = tmpval * w[26-(v)];			\
	tmpval = tmp[(v)] - tmp[17-(v)];			\
	ts[SBLIMIT*(8-(v))] = out1[8-(v)] + tmpval * w[8-(v)];	\
	ts[SBLIMIT*(9+(v))] = out1[9+(v)] + tmpval * w[9+(v)];	\
} while (0)

#ifndef USE_SIMD
static
#endif
void mpg123_dct36(real * inbuf, real * o1, real * o2, real * wintab, real * tsbuf)
{
	real tmp[18];
	register real *in = inbuf;

	in[17] += in[16];
	in[16] += in[15];
	in[15] += in[14];
	in[14] += in[13];
	in[13] += in[12];
	in[12] += in[11];
	in[11] += in[10];
	in[10] += in[9];
	in[9] += in[8];
	in[8] += in[7];
	in[7] += in[6];
	in[6] += in[5];
	in[5] += in[4];
	in[4] += in[3];
	in[3] += in[2];
	in[2] += in[1];
	in[1] += in[0];

	in[17] += in[15];
	in[15] += in[13];
	in[13] += in[11];
	in[11] += in[9];
	in[9] += in[7];
	in[7] += in[5];
	in[5] += in[3];
	in[3] += in[1];


	{
		real t3;
		{
			real t0, t1, t2;

			t0 = COS6_2 * (in[8] + in[16] - in[4]);
			t1 = COS6_2 * in[12];

			t3 = in[0];
			t2 = t3 - t1 - t1;
			tmp[1] = tmp[7] = t2 - t0;
			tmp[4] = t2 + t0 + t0;
			t3 += t1;

			t2 = COS6_1 * (in[10] + in[14] - in[2]);
			tmp[1] -= t2;
			tmp[7] += t2;
		}
		{
			real t0, t1, t2;

			t0 = cos9[0] * (in[4] + in[8]);
			t1 = cos9[1] * (in[8] - in[16]);
			t2 = cos9[2] * (in[4] + in[16]);

			tmp[2] = tmp[6] = t3 - t0 - t2;
			tmp[0] = tmp[8] = t3 + t0 + t1;
			tmp[3] = tmp[5] = t3 - t1 + t2;
		}
	}
	{
		real t1, t2, t3;

		t1 = cos18[0] * (in[2] + in[10]);
		t2 = cos18[1] * (in[10] - in[14]);
		t3 = COS6_1 * in[6];

		{
			real t0 = t1 + t2 + t3;
			tmp[0] += t0;
			tmp[8] -= t0;
		}

		t2 -= t3;
		t1 -= t3;

		t3 = cos18[2] * (in[2] + in[14]);

		t1 += t3;
		tmp[3] += t1;
		tmp[5] -= t1;

		t2 -= t3;
		tmp[2] += t2;
		tmp[6] -= t2;
	}

	{
		real t0, t1, t2, t3, t4, t5, t6, t7;

		t1 = COS6_2 * in[13];
		t2 = COS6_2 * (in[9] + in[17] - in[5]);

		t3 = in[1] + t1;
		t4 = in[1] - t1 - t1;
		t5 = t4 - t2;

		t0 = cos9[0] * (in[5] + in[9]);
		t1 = cos9[1] * (in[9] - in[17]);

		tmp[13] = (t4 + t2 + t2) * tfcos36[17 - 13];
		t2 = cos9[2] * (in[5] + in[17]);

		t6 = t3 - t0 - t2;
		t0 += t3 + t1;
		t3 += t2 - t1;

		t2 = cos18[0] * (in[3] + in[11]);
		t4 = cos18[1] * (in[11] - in[15]);
		t7 = COS6_1 * in[7];

		t1 = t2 + t4 + t7;
		tmp[17] = (t0 + t1) * tfcos36[17 - 17];
		tmp[9] = (t0 - t1) * tfcos36[17 - 9];
		t1 = cos18[2] * (in[3] + in[15]);
		t2 += t1 - t7;

		tmp[14] = (t3 + t2) * tfcos36[17 - 14];
		t0 = COS6_1 * (in[11] + in[15] - in[3]);
		tmp[12] = (t3 - t2) * tfcos36[17 - 12];

		t4 -= t1 + t7;

		tmp[16] = (t5 - t0) * tfcos36[17 - 16];
		tmp[10] = (t5 + t0) * tfcos36[17 - 10];
		tmp[15] = (t6 + t4) * tfcos36[17 - 15];
		tmp[11] = (t6 - t4) * tfcos36[17 - 11];
	}

	{
		register real *out2 = o2;
		register real *w = wintab;
		register real *out1 = o1;
		register real *ts = tsbuf;

		MACRO(0);
		MACRO(1);
		MACRO(2);
		MACRO(3);
		MACRO(4);
		MACRO(5);
		MACRO(6);
		MACRO(7);
		MACRO(8);
	}
}

/*
 * new DCT12
 */
static void dct12(real * in, real * rawout1, real * rawout2, register real * wi, register real * ts)
{

#define DCT12_PART1()				\
do {						\
	in5 = in[5*3];				\
	in5 += (in4 = in[4*3]);			\
	in4 += (in3 = in[3*3]);			\
	in3 += (in2 = in[2*3]);			\
	in2 += (in1 = in[1*3]);			\
	in1 += (in0 = in[0*3]);			\
						\
	in5 += in3; in3 += in1;			\
						\
	in2 *= COS6_1;				\
	in3 *= COS6_1;				\
} while (0) 

#define DCT12_PART2()				\
do {						\
	in0 += in4 * COS6_2;			\
						\
	in4 = in0 + in2;			\
	in0 -= in2;				\
						\
	in1 += in5 * COS6_2;			\
						\
	in5 = (in1 + in3) * tfcos12[0];		\
	in1 = (in1 - in3) * tfcos12[2];		\
						\
	in3 = in4 + in5;			\
	in4 -= in5;				\
						\
	in2 = in0 + in1;			\
	in0 -= in1;				\
} while (0)


	{
		real in0, in1, in2, in3, in4, in5;
		register real *out1 = rawout1;

		ts[SBLIMIT * 0] = out1[0];
		ts[SBLIMIT * 1] = out1[1];
		ts[SBLIMIT * 2] = out1[2];
		ts[SBLIMIT * 3] = out1[3];
		ts[SBLIMIT * 4] = out1[4];
		ts[SBLIMIT * 5] = out1[5];

		DCT12_PART1();

		{
			real tmp0, tmp1 = (in0 - in4);

			{
				real tmp2 = (in1 - in5) * tfcos12[1];

				tmp0 = tmp1 + tmp2;
				tmp1 -= tmp2;
			}
			ts[(17 - 1) * SBLIMIT] = out1[17 - 1] + tmp0 * wi[11 - 1];
			ts[(12 + 1) * SBLIMIT] = out1[12 + 1] + tmp0 * wi[6 + 1];
			ts[(6 + 1) * SBLIMIT] = out1[6 + 1] + tmp1 * wi[1];
			ts[(11 - 1) * SBLIMIT] = out1[11 - 1] + tmp1 * wi[5 - 1];
		}

		DCT12_PART2();

		ts[(17 - 0) * SBLIMIT] = out1[17 - 0] + in2 * wi[11 - 0];
		ts[(12 + 0) * SBLIMIT] = out1[12 + 0] + in2 * wi[6 + 0];
		ts[(12 + 2) * SBLIMIT] = out1[12 + 2] + in3 * wi[6 + 2];
		ts[(17 - 2) * SBLIMIT] = out1[17 - 2] + in3 * wi[11 - 2];

		ts[(6 + 0) * SBLIMIT] = out1[6 + 0] + in0 * wi[0];
		ts[(11 - 0) * SBLIMIT] = out1[11 - 0] + in0 * wi[5 - 0];
		ts[(6 + 2) * SBLIMIT] = out1[6 + 2] + in4 * wi[2];
		ts[(11 - 2) * SBLIMIT] = out1[11 - 2] + in4 * wi[5 - 2];
	}

	in++;

	{
		real in0, in1, in2, in3, in4, in5;
		register real *out2 = rawout2;

		DCT12_PART1();

		{
			real tmp0, tmp1 = (in0 - in4);

			{
				real tmp2 = (in1 - in5) * tfcos12[1];

				tmp0 = tmp1 + tmp2;
				tmp1 -= tmp2;
			}
			out2[5 - 1] = tmp0 * wi[11 - 1];
			out2[0 + 1] = tmp0 * wi[6 + 1];
			ts[(12 + 1) * SBLIMIT] += tmp1 * wi[1];
			ts[(17 - 1) * SBLIMIT] += tmp1 * wi[5 - 1];
		}

		DCT12_PART2();

		out2[5 - 0] = in2 * wi[11 - 0];
		out2[0 + 0] = in2 * wi[6 + 0];
		out2[0 + 2] = in3 * wi[6 + 2];
		out2[5 - 2] = in3 * wi[11 - 2];

		ts[(12 + 0) * SBLIMIT] += in0 * wi[0];
		ts[(17 - 0) * SBLIMIT] += in0 * wi[5 - 0];
		ts[(12 + 2) * SBLIMIT] += in4 * wi[2];
		ts[(17 - 2) * SBLIMIT] += in4 * wi[5 - 2];
	}

	in++;

	{
		real in0, in1, in2, in3, in4, in5;
		register real *out2 = rawout2;

		out2[12] = out2[13] = out2[14] = out2[15] = out2[16] = out2[17] = 0.0;

		DCT12_PART1();

		{
			real tmp0, tmp1 = (in0 - in4);

			{
				real tmp2 = (in1 - in5) * tfcos12[1];

				tmp0 = tmp1 + tmp2;
				tmp1 -= tmp2;
			}
			out2[11 - 1] = tmp0 * wi[11 - 1];
			out2[6 + 1] = tmp0 * wi[6 + 1];
			out2[0 + 1] += tmp1 * wi[1];
			out2[5 - 1] += tmp1 * wi[5 - 1];
		}

		DCT12_PART2();

		out2[11 - 0] = in2 * wi[11 - 0];
		out2[6 + 0] = in2 * wi[6 + 0];
		out2[6 + 2] = in3 * wi[6 + 2];
		out2[11 - 2] = in3 * wi[11 - 2];

		out2[0 + 0] += in0 * wi[0];
		out2[5 - 0] += in0 * wi[5 - 0];
		out2[0 + 2] += in4 * wi[2];
		out2[5 - 2] += in4 * wi[5 - 2];
	}
}

/*
 * III_hybrid
 */
static void III_hybrid(real fsIn[SBLIMIT][SSLIMIT],
		       real tsOut[SSLIMIT][SBLIMIT], int ch,
		       struct gr_info_s *gr_info, struct frame *fr)
{
	static real block[2][2][SBLIMIT * SSLIMIT] = { {{0,}} };
	static int blc[2] = { 0, 0 };

	real *tspnt = (real *) tsOut;
	real *rawout1, *rawout2;
	int bt, sb = 0;

	{
		int b = blc[ch];
		rawout1 = block[b][ch];
		b = -b + 1;
		rawout2 = block[b][ch];
		blc[ch] = b;
	}

	if (gr_info->mixed_block_flag)
	{
		sb = 2;
		DCT36(fsIn[0],rawout1,rawout2,win[0],tspnt);
		DCT36(fsIn[1],rawout1+18,rawout2+18,win1[0],tspnt+1);
		rawout1 += 36;
		rawout2 += 36;
		tspnt += 2;
	}

	bt = gr_info->block_type;
	if (bt == 2)
	{
		for (; sb < gr_info->maxb; sb += 2, tspnt += 2, rawout1 += 36, rawout2 += 36)
		{
			dct12(fsIn[sb], rawout1, rawout2, win[2], tspnt);
			dct12(fsIn[sb + 1], rawout1 + 18, rawout2 + 18, win1[2], tspnt + 1);
		}
	}
	else
	{
		for (; sb < gr_info->maxb; sb += 2, tspnt += 2, rawout1 += 36, rawout2 += 36)
		{
			DCT36(fsIn[sb], rawout1, rawout2, win[bt], tspnt);
			DCT36(fsIn[sb+1], rawout1+18, rawout2+18, win1[bt], tspnt+1);
		}
	}

	for (; sb < SBLIMIT; sb++, tspnt++)
	{
		int i;
		for (i = 0; i < SSLIMIT; i++)
		{
			tspnt[i * SBLIMIT] = *rawout1++;
			*rawout2++ = 0.0;
		}
	}
}

/*
 * main layer3 handler
 */
int mpg123_do_layer3(struct frame *fr)
{
	int gr, ch, ss;
	int scalefacs[2][39];	/* max 39 for short[13][3] mode, mixed: 38, long: 22 */
	struct III_sideinfo sideinfo;
	int stereo = fr->stereo;
	int single = fr->single;
	int ms_stereo, i_stereo;
	int sfreq = fr->sampling_frequency;
	int stereo1, granules;

	if (stereo == 1)
	{			/* stream is mono */
		stereo1 = 1;
		single = 0;
	}
	else if (single >= 0)	/* stream is stereo, but force to mono */
		stereo1 = 1;
	else
		stereo1 = 2;

	if (fr->mode == MPG_MD_JOINT_STEREO)
	{
		ms_stereo = (fr->mode_ext & 0x2) >> 1;
		i_stereo = fr->mode_ext & 0x1;
	}
	else
		ms_stereo = i_stereo = 0;

	granules = fr->lsf ? 1 : 2;
	if (!III_get_side_info(&sideinfo, stereo, ms_stereo, sfreq, single, fr->lsf))
		return 0;

	mpg123_set_pointer(sideinfo.main_data_begin);

	for (gr = 0; gr < granules; gr++)
	{
		real hybridIn[2][SBLIMIT][SSLIMIT];
		real hybridOut[2][SSLIMIT][SBLIMIT];

		{
			struct gr_info_s *gr_info = &(sideinfo.ch[0].gr[gr]);
			long part2bits;

			if (fr->lsf)
				part2bits = III_get_scale_factors_2(scalefacs[0], gr_info, 0);
			else
				part2bits = III_get_scale_factors_1(scalefacs[0], gr_info);

			if (III_dequantize_sample(hybridIn[0], scalefacs[0], gr_info, sfreq, part2bits))
				return 0;
		}

		if (stereo == 2)
		{
			struct gr_info_s *gr_info = &(sideinfo.ch[1].gr[gr]);
			long part2bits;

			if (fr->lsf)
				part2bits = III_get_scale_factors_2(scalefacs[1], gr_info, i_stereo);
			else
				part2bits = III_get_scale_factors_1(scalefacs[1], gr_info);

			if (III_dequantize_sample(hybridIn[1], scalefacs[1], gr_info, sfreq, part2bits))
				return 0;

			if (ms_stereo)
			{
				int i;
				int maxb = sideinfo.ch[0].gr[gr].maxb;

				if (sideinfo.ch[1].gr[gr].maxb > maxb)
					maxb = sideinfo.ch[1].gr[gr].maxb;
				for (i = 0; i < SSLIMIT * maxb; i++)
				{
					real tmp0 = ((real *) hybridIn[0])[i];
					real tmp1 = ((real *) hybridIn[1])[i];
					((real *) hybridIn[0])[i] = tmp0 + tmp1;
					((real *) hybridIn[1])[i] = tmp0 - tmp1;
				}
			}

			if (i_stereo)
				III_i_stereo(hybridIn, scalefacs[1], gr_info, sfreq, ms_stereo, fr->lsf);

			if (ms_stereo || i_stereo || (single == 3))
			{
				if (gr_info->maxb > sideinfo.ch[0].gr[gr].maxb)
					sideinfo.ch[0].gr[gr].maxb = gr_info->maxb;
				else
					gr_info->maxb = sideinfo.ch[0].gr[gr].maxb;
			}

			switch (single)
			{
				case 3:
					{
						register int i;
						register real *in0 = (real *) hybridIn[0],
							*in1 = (real *) hybridIn[1];
						for (i = 0; i < SSLIMIT * gr_info->maxb; 	     i++, in0++)
							*in0 = (*in0 + *in1++);	/* *0.5 done by pow-scale */
					}
					break;
				case 1:
					{
						register int i;
						register real *in0 = (real *) hybridIn[0],
							*in1 = (real *) hybridIn[1];
						for (i = 0; i < SSLIMIT * gr_info->maxb; i++)
							*in0++ = *in1++;
					}
					break;
			}
		}
		if (mpg123_info->eq_active)
		{
			int i, sb;

			if (single < 0)
			{
				for (sb = 0, i = 0; sb < SBLIMIT; sb++)
				{
					for (ss = 0; ss < SSLIMIT; ss++)
					{
						hybridIn[0][sb][ss] *= mpg123_info->eq_mul[i];
						hybridIn[1][sb][ss] *= mpg123_info->eq_mul[i++];
					}
				}
			}
			else
			{
				for (sb = 0, i = 0; sb < SBLIMIT; sb++)
				{
					for (ss = 0; ss < SSLIMIT; ss++)
						hybridIn[0][sb][ss] *= mpg123_info->eq_mul[i++];
				}
			}
		}

#ifdef USE_SIMD
		if (fr->synth_type == SYNTH_MMX && single < 0)
		{
			int i, sb;

			for (sb = 0, i = 0; sb < SBLIMIT; sb++)
			{
				for (ss = 0; ss < SSLIMIT; ss++)
				{
					hybridIn[0][sb][ss] *= 16384.0;
					hybridIn[1][sb][ss] *= 16384.0;
				}
			}
		}
#endif
		for (ch = 0; ch < stereo1; ch++)
		{
			struct gr_info_s *gr_info = &(sideinfo.ch[ch].gr[gr]);

			III_antialias(hybridIn[ch], gr_info);
			if (gr_info->maxb < 1 || gr_info->maxb > SBLIMIT)
				return 0;
			III_hybrid(hybridIn[ch], hybridOut[ch], ch, gr_info, fr);
		}

		for (ss = 0; ss < SSLIMIT; ss++)
		{
			if (single >= 0)
			{
				(fr->synth_mono) (hybridOut[0][ss], mpg123_pcm_sample, &mpg123_pcm_point);
			}
			else
			{
				int p1 = mpg123_pcm_point;

				(fr->synth) (hybridOut[0][ss], 0, mpg123_pcm_sample, &p1);
				(fr->synth) (hybridOut[1][ss], 1, mpg123_pcm_sample, &mpg123_pcm_point);
			}
		}

		if (mpg123_info->output_audio)
		{
			mpg123_ip.add_vis_pcm(mpg123_ip.output->written_time(),
					      mpg123_cfg.resolution == 16 ? FMT_S16_NE : FMT_U8,
					      mpg123_cfg.channels == 2 ? fr->stereo : 1,
					      mpg123_pcm_point, mpg123_pcm_sample);
			while (mpg123_ip.output->buffer_free() < mpg123_pcm_point &&
			       mpg123_info->going && mpg123_info->jump_to_time == -1)
				xmms_usleep(10000);
			if (mpg123_info->going && mpg123_info->jump_to_time == -1)
				mpg123_ip.output->write_audio(mpg123_pcm_sample, mpg123_pcm_point);
		}
		mpg123_pcm_point = 0;
	}
	return 1;
}