newmdct.c

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

/*
 *	MP3 window subband -> subband filtering -> mdct routine
 *
 *	Copyright (c) 1999 Takehiro TOMINAGA
 *
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/*
 *         Special Thanks to Patrick De Smet for your advices.
 */


#include "util.h"
#include "l3side.h"
#include "newmdct.h"

#define SCALE 32768

static FLOAT8 enwindow[] = 
{
  3.5780907e-02,1.7876148e-02,3.134727e-03,2.457142e-03,
    9.71317e-04,  2.18868e-04, 1.01566e-04,  1.3828e-05,

  3.5758972e-02, 3.401756e-03,  9.83715e-04,   9.9182e-05,
      -4.77e-07,  1.03951e-04,  9.53674e-04, 2.841473e-03,
     1.2398e-05,  1.91212e-04, 2.283096e-03,1.6994476e-02,
  1.8756866e-02, 2.630711e-03,  2.47478e-04,   1.4782e-05,

  3.5694122e-02, 3.643036e-03,  9.91821e-04,   9.6321e-05,
      -4.77e-07,  1.05858e-04,  9.30786e-04, 2.521515e-03,
     1.1444e-05,  1.65462e-04, 2.110004e-03,1.6112804e-02,
  1.9634247e-02, 2.803326e-03,  2.77042e-04,   1.6689e-05,

  3.5586357e-02, 3.858566e-03,  9.95159e-04,   9.3460e-05,
      -4.77e-07,  1.07288e-04,  9.02653e-04, 2.174854e-03,
     1.0014e-05,  1.40190e-04, 1.937389e-03,1.5233517e-02,
  2.0506859e-02, 2.974033e-03,  3.07560e-04,   1.8120e-05,

  3.5435200e-02, 4.049301e-03,  9.94205e-04,   9.0599e-05,
      -4.77e-07,  1.08242e-04,  8.68797e-04, 1.800537e-03,
      9.060e-06,  1.16348e-04, 1.766682e-03,1.4358521e-02,
  2.1372318e-02,  3.14188e-03,  3.39031e-04,   1.9550e-05,

  3.5242081e-02, 4.215240e-03,  9.89437e-04,   8.7261e-05,
      -4.77e-07,  1.08719e-04,  8.29220e-04, 1.399517e-03,
      8.106e-06,   9.3937e-05, 1.597881e-03,1.3489246e-02,
  2.2228718e-02, 3.306866e-03,  3.71456e-04,   2.1458e-05,

  3.5007000e-02, 4.357815e-03,  9.80854e-04,   8.3923e-05,
      -4.77e-07,  1.08719e-04,   7.8392e-04,  9.71317e-04,
      7.629e-06,   7.2956e-05, 1.432419e-03,1.2627602e-02,
  2.3074150e-02, 3.467083e-03,  4.04358e-04,   2.3365e-05,

  3.4730434e-02, 4.477024e-03,  9.68933e-04,   8.0585e-05,
      -9.54e-07,  1.08242e-04,  7.31945e-04,  5.15938e-04,
      6.676e-06,   5.2929e-05, 1.269817e-03,1.1775017e-02,
  2.3907185e-02, 3.622532e-03,  4.38213e-04,   2.5272e-05,

  3.4412861e-02, 4.573822e-03,  9.54151e-04,   7.6771e-05,
      -9.54e-07,  1.06812e-04,  6.74248e-04,   3.3379e-05,
      6.199e-06,   3.4332e-05, 1.111031e-03,1.0933399e-02,
  2.4725437e-02, 3.771782e-03,  4.72546e-04,   2.7657e-05,

  3.4055710e-02, 4.649162e-03,  9.35555e-04,   7.3433e-05,
      -9.54e-07,  1.05381e-04,  6.10352e-04, -4.75883e-04,
      5.245e-06,   1.7166e-05,  9.56535e-04,1.0103703e-02,
  2.5527000e-02, 3.914356e-03,  5.07355e-04,   3.0041e-05,

  3.3659935e-02, 4.703045e-03,  9.15051e-04,   7.0095e-05,
      -9.54e-07,  1.02520e-04,  5.39303e-04,-1.011848e-03,
      4.768e-06,     9.54e-07,  8.06808e-04, 9.287834e-03,
  2.6310921e-02, 4.048824e-03,  5.42164e-04,   3.2425e-05,

  3.3225536e-02, 4.737377e-03,  8.91685e-04,   6.6280e-05,
     -1.431e-06,   9.9182e-05,  4.62532e-04,-1.573563e-03,
      4.292e-06,  -1.3828e-05,  6.61850e-04, 8.487225e-03,
  2.7073860e-02, 4.174709e-03,  5.76973e-04,   3.4809e-05,

  3.2754898e-02, 4.752159e-03,  8.66413e-04,   6.2943e-05,
     -1.431e-06,   9.5367e-05,  3.78609e-04,-2.161503e-03,
      3.815e-06,   -2.718e-05,  5.22137e-04, 7.703304e-03,
  2.7815342e-02, 4.290581e-03,  6.11782e-04,   3.7670e-05,

  3.2248020e-02, 4.748821e-03,  8.38757e-04,   5.9605e-05,
     -1.907e-06,   9.0122e-05,  2.88486e-04,-2.774239e-03,
      3.338e-06,  -3.9577e-05,  3.88145e-04, 6.937027e-03,
  2.8532982e-02, 4.395962e-03,  6.46591e-04,   4.0531e-05,

  3.1706810e-02, 4.728317e-03,  8.09669e-04,    5.579e-05,
     -1.907e-06,   8.4400e-05,  1.91689e-04,-3.411293e-03,
      3.338e-06,  -5.0545e-05,  2.59876e-04, 6.189346e-03,
  2.9224873e-02, 4.489899e-03,  6.80923e-04,   4.3392e-05,

  3.1132698e-02, 4.691124e-03,  7.79152e-04,   5.2929e-05,
     -2.384e-06,   7.7724e-05,   8.8215e-05,-4.072189e-03,
      2.861e-06,  -6.0558e-05,  1.37329e-04, 5.462170e-03,
  2.9890060e-02, 4.570484e-03,  7.14302e-04,   4.6253e-05,

  3.0526638e-02, 4.638195e-03,  7.47204e-04,   4.9591e-05,
   4.756451e-03,   2.1458e-05,  -6.9618e-05,    2.384e-06
};

static FLOAT8 sb_sample[2][2][18][SBLIMIT];
static FLOAT8 mm[16][SBLIMIT - 1];

#define NS 12
#define NL 36

static const int all[] = {0,2,3,5,6,8,9,11,12,14,15,17};
static FLOAT8 ca[8], cs[8];
static FLOAT8 cos_s[NS / 2][NS / 2];
static FLOAT8 cos_l[(NL / 2) * 12 + (NL / 6) * 4 + (NL / 18) * 2];
static FLOAT8 win[4][36];

#define work (&win[2][4])

/************************************************************************
*
* window_subband()
*
* PURPOSE:  Overlapping window on PCM samples
*
* SEMANTICS:
* 32 16-bit pcm samples are scaled to fractional 2's complement and
* concatenated to the end of the window buffer #x#. The updated window
* buffer #x# is then windowed by the analysis window #c# to produce the
* windowed sample #z#
*
************************************************************************/

static void window_subband(short *xk, FLOAT8 d[SBLIMIT], FLOAT8 *in)
{
    int i;
    FLOAT8 s, t, *wp;
    wp = enwindow;
    {
	t  =  xk[255];
	t += (xk[223] - xk[287]) * *wp++;
	t += (xk[191] + xk[319]) * *wp++;
	t += (xk[159] - xk[351]) * *wp++;
	t += (xk[127] + xk[383]) * *wp++;
	t += (xk[ 95] - xk[415]) * *wp++;
	t += (xk[ 63] + xk[447]) * *wp++;
	t += (xk[ 31] - xk[479]) * *wp++;
	in[15] = t;
    }

    for (i = 14; i >= 0; --i) {
	short *x1 = &xk[i];
	short *x2 = &xk[-i];
	FLOAT8 w;

	s = x2[270]; t = x1[240];
	w = *wp++; s += x2[334] * w; t += x1[176] * w;
	w = *wp++; s += x2[398] * w; t += x1[112] * w;
	w = *wp++; s += x2[462] * w; t += x1[ 48] * w;
	w = *wp++; s += x2[ 14] * w; t += x1[496] * w;
	w = *wp++; s += x2[ 78] * w; t += x1[432] * w;
	w = *wp++; s += x2[142] * w; t += x1[368] * w;
	w = *wp++; s += x2[206] * w; t += x1[304] * w;

	w = *wp++; s += x1[ 16] * w; t -= x2[494] * w;
	w = *wp++; s += x1[ 80] * w; t -= x2[430] * w;
	w = *wp++; s += x1[144] * w; t -= x2[366] * w;
	w = *wp++; s += x1[208] * w; t -= x2[302] * w;
	w = *wp++; s -= x1[272] * w; t += x2[238] * w;
	w = *wp++; s -= x1[336] * w; t += x2[174] * w;
	w = *wp++; s -= x1[400] * w; t += x2[110] * w;
	w = *wp++; s -= x1[464] * w; t += x2[ 46] * w;

	in[30 - i] = s;
	in[i] = t;
    }

    {
	s  = xk[239];
	s += xk[175] * *wp++;
	s += xk[111] * *wp++;
	s += xk[ 47] * *wp++;
	s -= xk[303] * *wp++;
	s -= xk[367] * *wp++;
	s -= xk[431] * *wp++;
	s -= xk[495] * *wp++;
	/* in[-1] = s;  */
    }

    in++;
    wp = &mm[0][0];
    for (i = 15; i >= 0; --i) {
	int j;
	FLOAT8 s0 = s; /* mm[i][0] is always 1 */
	FLOAT8 s1 = t * *wp++;
	for (j = 14; j >= 0; j--) {
	    s0 += *wp++ * *in++;
	    s1 += *wp++ * *in++;
	}
	in -= 30;
	d[i     ] = s0 + s1;
	d[31 - i] = s0 - s1;
    }
}


/*-------------------------------------------------------------------*/
/*                                                                   */
/*   Function: Calculation of the MDCT                               */
/*   In the case of long blocks (type 0,1,3) there are               */
/*   36 coefficents in the time domain and 18 in the frequency       */
/*   domain.                                                         */
/*   In the case of short blocks (type 2) there are 3                */
/*   transformations with short length. This leads to 12 coefficents */
/*   in the time and 6 in the frequency domain. In this case the     */
/*   results are stored side by side in the vector out[].            */
/*                                                                   */
/*   New layer3                                                      */
/*                                                                   */
/*-------------------------------------------------------------------*/

static void mdct_short(FLOAT8 *out, FLOAT8 *in)
{
    int m;
    for (m = NS / 2 - 1; m >= 0; --m) {
	int l;
	FLOAT8 a0, a1, a2, a3, a4, a5;
	a0 = cos_s[m][0];
	a1 = cos_s[m][1];
	a2 = cos_s[m][2];
	a3 = cos_s[m][3];
	a4 = cos_s[m][4];
	a5 = cos_s[m][5];
	for (l = 2; l >= 0; l--) {
	    out[3 * m + l] =
		a0 * in[6 * l    ] +
		a1 * in[6 * l + 1] +
		a2 * in[6 * l + 2] +
		a3 * in[6 * l + 3] +
		a4 * in[6 * l + 4] +
		a5 * in[6 * l + 5];
	}
    }
}

static void mdct_long(FLOAT8 *out, FLOAT8 *in)
{
    FLOAT8 s0, s1, s2, s3, s4, s5;
    int j = sizeof(all) / sizeof(int) - 1;
    FLOAT8 *cos_l0 = cos_l;
    do {
	out[all[j]] =
	    in[ 0] * cos_l0[ 0] +
	    in[ 1] * cos_l0[ 1] +
	    in[ 2] * cos_l0[ 2] +
	    in[ 3] * cos_l0[ 3] +
	    in[ 4] * cos_l0[ 4] +
	    in[ 5] * cos_l0[ 5] +
	    in[ 6] * cos_l0[ 6] +
	    in[ 7] * cos_l0[ 7] +
	    in[ 8] * cos_l0[ 8] +
	    in[ 9] * cos_l0[ 9] +
	    in[10] * cos_l0[10] +
	    in[11] * cos_l0[11] +
	    in[12] * cos_l0[12] +
	    in[13] * cos_l0[13] +
	    in[14] * cos_l0[14] +
	    in[15] * cos_l0[15] +
	    in[16] * cos_l0[16] +
	    in[17] * cos_l0[17];
	cos_l0 += 18;
    } while (--j >= 0);

    s0 = in[0] + in[ 5] + in[15];
    s1 = in[1] + in[ 4] + in[16];
    s2 = in[2] + in[ 3] + in[17];
    s3 = in[6] - in[ 9] + in[14];
    s4 = in[7] - in[10] + in[13];
    s5 = in[8] - in[11] + in[12];

    /* 16 */
    out[16] =
	s0 * cos_l0[0] + s1 * cos_l0[1] + s2 * cos_l0[2] +
	s3 * cos_l0[3] + s4 * cos_l0[4] + s5 * cos_l0[5];
    cos_l0 += 6;

    /* 10 */
    out[10] =
	s0 * cos_l0[0] + s1 * cos_l0[1] + s2 * cos_l0[2] +
	s3 * cos_l0[3] + s4 * cos_l0[4] + s5 * cos_l0[5];
    cos_l0 += 6;

    /* 7 */
    out[7] =
	s0 * cos_l0[0] + s1 * cos_l0[1] + s2 * cos_l0[2] +
	s3 * cos_l0[3] + s4 * cos_l0[4] + s5 * cos_l0[5];
    cos_l0 += 6;

    /* 1 */
    out[1] =
	s0 * cos_l0[0] + s1 * cos_l0[1] + s2 * cos_l0[2] +
	s3 * cos_l0[3] + s4 * cos_l0[4] + s5 * cos_l0[5];
    cos_l0 += 6;