mpegimdc.c

mp3解码源代码

/****************************************************************************/

/*
 *	mpegimdc.c -- MPEG Audio IMDCT hybrid filter
 *
 *	Author  :   St閜hane TAVENARD
 *
 *	(C) Copyright 1997-1999 St閜hane TAVENARD
 *	    All Rights Reserved
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; either version 2 of the License, or
 *	(at your option) any later version.
 *	
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *	
 *	You should have received a copy of the GNU General Public License
 *	along with this program; if not, write to the Free Software
 *	Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/****************************************************************************/

#include "defs.h"
#include "mpegaud.h"
#include "mpegimdc.h"
#include <math.h>

#ifdef MPEGAUD_INT
#define IMDCT_INT
#ifdef USE_IMDCT_TABLE
#include "imdct.h"
#endif
#else
// This one is faster for float (but not accurate for int math)
#define NEW_IMDCT	
#endif

#ifdef IMDCT_TABLE_GEN
#include <stdlib.h>
#endif

/****************************************************************************/

#ifdef NEW_IMDCT

#define SBLIMIT	MPA_SBLIMIT
#define SSLIMIT	MPA_SSLIMIT

#ifndef M_PI
#define M_PI	3.14159265358979
#endif

#ifdef IMDCT_INT

#define WIN_BITS 14
#define TFC_BITS 14
#define COS_BITS 14

#define WIN_MULT(val, windex)	(((val) * w[windex]) >> WIN_BITS)
#define T36_MULT(val, tindex)	(((val) * tfcos36[tindex]) >> TFC_BITS)
#define T12_MULT(val, tindex)	(((val) * tfcos12[tindex]) >> TFC_BITS)
#define COS_MULT(val, cosval)	(((val) * (cosval)) >> COS_BITS)
#define CM(val, cindex)		((val) * c[cindex]))
#define CM1(val)		((val) << COS_BITS)
#define CF(val)			((val) >> (COS_BITS+SH_BITS))

#define IMDCT_TYPE	INT32
#define WIN_TYPE	INT32
#define COS_TYPE	INT16

#else

#define WIN_MULT(val, windex)	((val) * w[windex])
#define T36_MULT(val, tindex)	((val) * tfcos36[tindex ])
#define T12_MULT(val, tindex)	((val) * tfcos12[tindex ])
#define COS_MULT(val, cosval)	((val) * (cosval))
#define CM(val, cindex)		((val) * c[cindex])
#define CM1(val)		(val)
#define CF(val)			(val)

#define IMDCT_TYPE REAL
#define WIN_TYPE REAL
#define COS_TYPE REAL

#endif

/****************************************************************************/

static COS_TYPE COS1[12][6];
static WIN_TYPE win[4][36];
static WIN_TYPE win1[4][36];
static COS_TYPE COS9[9];
static COS_TYPE COS6_1, COS6_2;
static IMDCT_TYPE tfcos36[9];
static IMDCT_TYPE tfcos12[3];

/****************************************************************************/

static void dct36(MPEGAUD_FRACT_TYPE *inbuf, MPEGIMDCT_BLOCK_TYPE *prev,
	WIN_TYPE *wintab, MPEGAUD_FRACT_TYPE *tsbuf)
{
	MPEGAUD_FRACT_TYPE *in = inbuf;

#define I0(i)	in[2*(i)+0]
#define I1(i)	in[2*(i)+1]

	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];

	{

#define MACRO0(v) { \
	IMDCT_TYPE tmp = sum0 + sum1; \
	sum0 -= sum1; \
	ts[SBLIMIT*(8-(v))] = prev[8-(v)] + WIN_MULT(sum0, 8-(v)); \
	prev[8-(v)] = WIN_MULT(tmp, 26-(v)); \
	ts[SBLIMIT*(9+(v))] = prev[9+(v)] + WIN_MULT(sum0, 9+(v)); \
	prev[9+(v)] = WIN_MULT(tmp, 27+(v));
}

#define MACRO1(v) { \
	IMDCT_TYPE sum0, sum1; \
	sum0 = tmp1a + tmp2a; \
	sum1 = T36_MULT(tmp1b + tmp2b, v); \
	MACRO0(v);
}

#define MACRO2(v) { \
	IMDCT_TYPE sum0, sum1; \
	sum0 = tmp2a - tmp1a; \
	sum1 = T36_MULT(tmp2b - tmp1b, v); \
	MACRO0(v);
}

		const COS_TYPE *c = COS9;
		WIN_TYPE *w = wintab;
		MPEGAUD_FRACT_TYPE *ts = tsbuf;
		IMDCT_TYPE ta33, ta66, tb33, tb66;

		ta33 = CM(I0(3), 3);
		ta66 = CM(I0(6), 6);
		tb33 = CM(I1(3), 3);
		tb66 = CM(I1(6), 6);

		{
			IMDCT_TYPE tmp1a, tmp2a, tmp1b, tmp2b;

			tmp1a = CF(CM(I0(1), 1) + ta33 + CM(I0(5), 5) +
				   CM(I0(7), 7));
			tmp1b = CF(CM(I1(1), 1) + tb33 + CM(I1(5), 5) +
				   CM(I1(7), 7));
			tmp2a = CF(CM1(I0(0)) + CM(I0(2), 2) +
				   CM(I0(4), 4) + ta66 + CM(I0(8), 8));
			tmp2b = CF(CM1(I1(0)) + CM(I1(2), 2) +
				   CM(I1(4), 4) + tb66 + CM(I1(8), 8));
			MACRO1(0);
			MACRO2(8);
		}

		{
			IMDCT_TYPE tmp1a, tmp2a, tmp1b, tmp2b;
			tmp1a = CF(CM(I0(1) - I0(5) - I0(7), 3));
			tmp1b = CF(CM(I1(1) - I1(5) - I1(7), 3));
			tmp2a = CF(CM(I0(2) - I0(4) - I0(8), 6) +
				   CM1(I0(0) - I0(6)));
			tmp2b = CF(CM(I1(2) - I1(4) - I1(8), 6) +
				   CM1(I1(0) - I1(6)));
			MACRO1(1);
			MACRO2(7);
		}

		{
			IMDCT_TYPE tmp1a, tmp2a, tmp1b, tmp2b;
			tmp1a = CF(CM(I0(1), 5) - ta33 - CM(I0(5), 7) +
				   CM(I0(7), 1));
			tmp1b = CF(CM(I1(1), 5) - tb33 - CM(I1(5), 7) +
				   CM(I1(7), 1));
			tmp2a = CF(CM1(I0(0)) - CM(I0(2), 8) -
				   CM(I0(4), 2) + ta66 + CM(I0(8), 4));
			tmp2b = CF(CM1(I1(0)) - CM(I1(2), 8) -
				   CM(I1(4), 2) + tb66 + CM(I1(8), 4));
			MACRO1(2);
			MACRO2(6);
		}

		{
			IMDCT_TYPE tmp1a, tmp2a, tmp1b, tmp2b;
			tmp1a = CF(CM(I0(1), 7) - ta33 + CM(I0(5), 1) -
				   CM(I0(7), 5));
			tmp1b = CF(CM(I1(1), 7) - tb33 + CM(I1(5), 1) -
				   CM(I1(7), 5));
			tmp2a = CF(CM1(I0(0)) - CM(I0(2), 4) +
				   CM(I0(4), 8) + ta66 - CM(I0(8), 2));
			tmp2b = CF(CM1(I1(0)) - CM(I1(2), 4) +
				   CM(I1(4), 8) + tb66 - CM(I1(8), 2));
			MACRO1(3);
			MACRO2(5);
		}

		{
			IMDCT_TYPE sum0, sum1;
			sum0 = CF(CM1 (I0(0) - I0(2) + I0(4) - I0(6) + I0(8)));
			sum1 = T36_MULT(CF(CM1(I1(0) - I1(2) + I1(4) - I1(6) +
				I1(8))), 4);
			MACRO0(4);
		}
	}
}

/****************************************************************************/

static void dct12(MPEGAUD_FRACT_TYPE *in, MPEGIMDCT_BLOCK_TYPE *prev,
	WIN_TYPE *w, MPEGAUD_FRACT_TYPE *ts)
{
#define DCT12_PART1 \
	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 = COS_MULT(in2, COS6_1); \
	in3 = COS_MULT(in3, COS6_1);

#define DCT12_PART2 \
	in0 += COS_MULT(in4, COS6_2); \
	\
	in4 = in0 + in2; \
	in0 -= in2; \
	\
	in1 += COS_MULT(in5, COS6_2); \
	\
	in5 = T12_MULT((in1 + in3), 0); \
	in1 = T12_MULT((in1 - in3), 2); \
	\
	in3 = in4 + in5; \
	in4 -= in5; \
	\
	in2 = in0 + in1; \
	in0 -= in1;

	{
		IMDCT_TYPE in0, in1, in2, in3, in4, in5;
		ts[SBLIMIT * 0] = prev[0];
		ts[SBLIMIT * 1] = prev[1];
		ts[SBLIMIT * 2] = prev[2];
		ts[SBLIMIT * 3] = prev[3];
		ts[SBLIMIT * 4] = prev[4];
		ts[SBLIMIT * 5] = prev[5];

		DCT12_PART1

		{
			IMDCT_TYPE tmp0, tmp1 = (in0 - in4);
			{
				IMDCT_TYPE tmp2 = T12_MULT((in1 - in5), 1);
				tmp0 = tmp1 + tmp2;
				tmp1 -= tmp2;
			}
			ts[(17 - 1) * SBLIMIT] =
				prev[17 - 1] + WIN_MULT(tmp0, 11 - 1);
			ts[(12 + 1) * SBLIMIT] =
				prev[12 + 1] + WIN_MULT(tmp0, 6 + 1);
			ts[(6 + 1) * SBLIMIT] =
				prev[6 + 1] + WIN_MULT(tmp1, 1);
			ts[(11 - 1) * SBLIMIT] =
				prev[11 - 1] + WIN_MULT(tmp1, 5 - 1);
		}

		DCT12_PART2

		ts[(17 - 0) * SBLIMIT] = prev[17 - 0] + WIN_MULT(in2, 11 - 0);
		ts[(12 + 0) * SBLIMIT] = prev[12 + 0] + WIN_MULT(in2, 6 + 0);
		ts[(12 + 2) * SBLIMIT] = prev[12 + 2] + WIN_MULT(in3, 6 + 2);
		ts[(17 - 2) * SBLIMIT] = prev[17 - 2] + WIN_MULT(in3, 11 - 2);

		ts[(6 + 0) * SBLIMIT] = prev[6 + 0] + WIN_MULT(in0, 0);
		ts[(11 - 0) * SBLIMIT] = prev[11 - 0] + WIN_MULT(in0, 5 - 0);
		ts[(6 + 2) * SBLIMIT] = prev[6 + 2] + WIN_MULT(in4, 2);
		ts[(11 - 2) * SBLIMIT] = prev[11 - 2] + WIN_MULT(in4, 5 - 2);
	}

	in++;

	{
		IMDCT_TYPE in0, in1, in2, in3, in4, in5;

		DCT12_PART1

		{
			IMDCT_TYPE tmp0, tmp1 = (in0 - in4);
			{
				IMDCT_TYPE tmp2 = T12_MULT((in1 - in5), 1);
				tmp0 = tmp1 + tmp2;
				tmp1 -= tmp2;
			}
			prev[5 - 1] = WIN_MULT(tmp0, 11 - 1);
			prev[0 + 1] = WIN_MULT(tmp0, 6 + 1);
			ts[(12 + 1) * SBLIMIT] += WIN_MULT(tmp1, 1);
			ts[(17 - 1) * SBLIMIT] += WIN_MULT(tmp1, 5 - 1);
		}

		DCT12_PART2

		prev[5 - 0] = WIN_MULT(in2, 11 - 0);
		prev[0 + 0] = WIN_MULT(in2, 6 + 0);
		prev[0 + 2] = WIN_MULT(in3, 6 + 2);
		prev[5 - 2] = WIN_MULT(in3, 11 - 2);

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

	in++;

	{
		IMDCT_TYPE in0, in1, in2, in3, in4, in5;
		prev[12] = prev[13] = prev[14] = prev[15] = prev[16] =
			prev[17] = (MPEGIMDCT_BLOCK_TYPE) 0;

		DCT12_PART1

		{
			IMDCT_TYPE tmp0, tmp1 = (in0 - in4);
			{
				IMDCT_TYPE tmp2 = T12_MULT((in1 - in5), 1);
				tmp0 = tmp1 + tmp2;
				tmp1 -= tmp2;
			}
			prev[11 - 1] = WIN_MULT(tmp0, 11 - 1);
			prev[6 + 1] = WIN_MULT(tmp0, 6 + 1);
			prev[0 + 1] += WIN_MULT(tmp1, 1);
			prev[5 - 1] += WIN_MULT(tmp1, 5 - 1);
		}

		DCT12_PART2

		prev[11 - 0] = WIN_MULT(in2, 11 - 0);
		prev[6 + 0] = WIN_MULT(in2, 6 + 0);
		prev[6 + 2] = WIN_MULT(in3, 6 + 2);
		prev[11 - 2] = WIN_MULT(in3, 11 - 2);

		prev[0 + 0] += WIN_MULT(in0, 0);
		prev[5 - 0] += WIN_MULT(in0, 5 - 0);
		prev[0 + 2] += WIN_MULT(in4, 2);
		prev[5 - 2] += WIN_MULT(in4, 5 - 2);
	}
}

/****************************************************************************/

/*
 *	Apply the hybrid imdct to a granule
 *	Return 0 if Ok
*/

int MPEGIMDCT_hybrid(MPEGIMDCT *mpegimdct, MPEGAUD_FRACT_TYPE *in,
	MPEGAUD_FRACT_TYPE *out, INT16 block_type,
	BOOL mixed, INT16 ch, INT16 sb_max)
{
	static BOOL init = FALSE;
	MPEGIMDCT_BLOCK_TYPE *prev;
	INT16 bt1, bt2;
	INT16 sb;

#define O(i)	out[i*32]

	if (!init) {
		int i, j;

#ifdef IMDCT_INT
#define KW	((double)(1<<(WIN_BITS-1))+0)
#define KC	(double)(1<<COS_BITS)
#define KT	((double)(1<<(TFC_BITS-1))+0)
#else
#define KW	0.5
#define KC	1
#define KT	0.5
#endif

		for (i = 0; i < 18; i++) {
			win[0][i] = win[1][i] = KW * sin(M_PI / 72.0 *
				(double) (2 * (i + 0) + 1)) / cos(M_PI *
				(double) (2 * (i + 0) + 19) / 72.0);
			win[0][i + 18] = win[3][i + 18] = KW * sin(M_PI / 72.0 *
				(double) (2 * (i + 18) + 1)) / cos(M_PI *
				(double) (2 * (i + 18) + 19) / 72.0);
		}
		for (i = 0; i < 6; i++) {
			win[1][i + 18] = KW / cos(M_PI *
				(double) (2 * (i + 18) + 19) / 72.0);
			win[3][i + 12] = KW / cos(M_PI *
				(double) (2 * (i + 12) + 19) / 72.0);
			win[1][i + 24] = KW * sin(M_PI / 24.0 *
				(double) (2 * i + 13)) / cos(M_PI *
				(double) (2 * (i + 24) + 19) / 72.0);
			win[1][i + 30] = win[3][i] = 0.0;
			win[3][i + 6] = KW * sin(M_PI / 24.0 *
				(double) (2 * i + 1)) / cos(M_PI *
				(double) (2 * (i + 6) + 19) / 72.0);
		}

		for (i = 0; i < 9; i++)
			COS9[i] = KC * cos(M_PI / 18.0 * (double) i);

		for (i = 0; i < 9; i++)
			tfcos36[i] =
				KT / cos(M_PI * (double) (i * 2 + 1) / 36.0);
		for (i = 0; i < 3; i++)
			tfcos12[i] =
				KT / cos(M_PI * (double) (i * 2 + 1) / 12.0);

		COS6_1 = KC * cos(M_PI / 6.0 * (double) 1);
		COS6_2 = KC * cos(M_PI / 6.0 * (double) 2);

		for (i = 0; i < 12; i++) {
			win[2][i] = KW * sin(M_PI / 24.0 *
				(double) (2 * i + 1)) / cos(M_PI *
				(double) (2 * i + 7) / 24.0);
			for (j = 0; j < 6; j++)
				COS1[i][j] = KC * cos(M_PI / 24.0 *
					(double) ((2 * i + 7) * (2 * j + 1)));
		}

		for (j = 0; j < 4; j++) {
			static const int len[4] = { 36, 36, 12, 36 };
			for (i = 0; i < len[j]; i += 2)
				win1[j][i] = +win[j][i];
			for (i = 1; i < len[j]; i += 2)
				win1[j][i] = -win[j][i];
		}
		init = TRUE;
	}

	prev = mpegimdct->prevblk[ch];

	bt1 = (mixed) ? 0 : block_type;
	bt2 = block_type;

	if (bt2 == 2) {
		if (!bt1) {
			dct36(in, prev, win[0], out);
			in += SSLIMIT;
			dct36(in, prev + 18, win1[0], out + 1);
			in += SSLIMIT;
		} else {
			dct12(in, prev, win[2], out);
			in += SSLIMIT;
			dct12(in, prev + 18, win1[2], out + 1);
			in += SSLIMIT;
		}
		prev += 36;
		out += 2;
		for (sb = 2; sb < sb_max; sb += 2, out += 2, prev += 36) {
			dct12(in, prev, win[2], out);
			in += SSLIMIT;
			dct12(in, prev + 18, win1[2], out + 1);
			in += SSLIMIT;
		}
	} else {
		dct36(in, prev, win[bt1], out);
		in += SSLIMIT;
		dct36(in, prev + 18, win1[bt1], out + 1);
		in += SSLIMIT;
		prev += 36;
		out += 2;
		for (sb = 2; sb < sb_max; sb += 2, out += 2, prev += 36) {
			dct36(in, prev, win[bt2], out);
			in += SSLIMIT;
			dct36(in, prev + 18, win1[bt2], out + 1);
			in += SSLIMIT;
		}
	}

	for (; sb < SBLIMIT; sb++) {
		O(0) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(1) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(2) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(3) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(4) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(5) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(6) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(7) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(8) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(9) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(10) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(11) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(12) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(13) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(14) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(15) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(16) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		O(17) = *prev;
		*prev++ = (MPEGIMDCT_BLOCK_TYPE) 0;
		out++;
	}

	return(0);
}

/****************************************************************************/
#else /* !NEW_IMDCT */
/****************************************************************************/

#ifndef PI
#define PI	3.14159265358979323846
#endif

#ifdef IMDCT_INT
#define IMDCT_IO_TYPE	MPEGAUD_FRACT_TYPE
#define IMDCT_TYPE	INT32
#define IMDCT_BITS	14
#define WIN_TYPE	INT16
#define WIN_BITS	14
#define WIN_MULT(t, w)	(((t) * (w))>>WIN_BITS)
#else
#define IMDCT_IO_TYPE	REAL
#define IMDCT_TYPE	REAL
#define IMDCT_BITS	0
#define WIN_TYPE	REAL
#define WIN_BITS	0
#define WIN_MULT(t, w)	((t) * (w))
#endif

/****************************************************************************/

static int imdct_l(IMDCT_IO_TYPE *x, IMDCT_IO_TYPE *out,
	MPEGIMDCT_BLOCK_TYPE *prev, WIN_TYPE *win)