mpegimdc.c

Application (fix point) for playing MP3 data on embedded systems. mp3play is designed to be able to

/*-----------------------------------------------------------------------------

    File    :   MPEGIMDC.C

    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.

    #Rev|   Date   |                      Comment
    ----|----------|--------------------------------------------------------
    0   |27/05/1997| Initial revision                                     ST
    1   |05/06/1997| Added FPU Optimized version                          ST
    2   |02/05/1998| PPC Support                                          ST
    3   |24/05/1999| Added some const                                     ST
    4   |21/10/1999| Added imdct table generation/usage                   RS

    ------------------------------------------------------------------------

    MPEG Audio IMDCT hybrid filter

-----------------------------------------------------------------------------*/


#include "defs.h"
#include "mpegaud.h"
#include "mpegimdc.h"
#if defined(ASM_OPTIMIZE) || defined(COLDFIRE_ASM)
#ifdef MPEGAUD_INT
#include "mpegimda.h"
#else
#include "mpegimdf.h" // #1
#endif
#endif
#include <math.h>

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

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

#if defined(ASM_OPTIMIZE) || defined(COLDFIRE_ASM)

int MPEGIMDCT_hybrid( MPEGIMDCT *mpegimdct, MPEGAUD_FRACT_TYPE *in,
                      MPEGAUD_FRACT_TYPE *out, INT16 block_type, BOOL mixed,
                      INT16 ch, INT16 sb_max )
/*--------------------------------------------------------------------------
   Apply the hybrid imdct to a granule
   Return 0 if Ok
*/
{
#ifdef MPEGAUD_INT
   MPEGIMDA_hybrid( in, out, mpegimdct->prevblk[ ch ], block_type, mixed, sb_max );
#else
   MPEGIMDF_hybrid( in, out, mpegimdct->prevblk[ ch ], block_type, mixed, sb_max );
#endif

MPEGAUD_CHECK_DEMO;

   return 0;

} /* MPEGIMDCT_hybrid */


#else // -> NOT ASM_OPTIMIZE


#ifdef NEW_IMDCT

#define SBLIMIT MPA_SBLIMIT
#define SSLIMIT MPA_SSLIMIT

#ifndef M_PI // #2
#define M_PI 3.14159265358979
#endif

#ifdef IMDCT_INT

#define WIN_BITS 14 // 14
#define TFC_BITS 14 // 14
#define COS_BITS 14 // 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 )
{
  register 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); }

    register const COS_TYPE *c = COS9;
    register WIN_TYPE *w = wintab;
    register 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 );