idctf.c

MPEG-2 has 7 distinct parts as well. The first part is the Systems section which defines the contain

/* idctf.c, inverse fast discrete cosine transform, floating point version  */

/* Copyright (C) 1994, MPEG Software Simulation Group. All Rights Reserved. */

/*
 * Disclaimer of Warranty
 *
 * These software programs are available to the user without any license fee or
 * royalty on an "as is" basis.  The MPEG Software Simulation Group disclaims
 * any and all warranties, whether express, implied, or statuary, including any
 * implied warranties or merchantability or of fitness for a particular
 * purpose.  In no event shall the copyright-holder be liable for any
 * incidental, punitive, or consequential damages of any kind whatsoever
 * arising from the use of these programs.
 *
 * This disclaimer of warranty extends to the user of these programs and user's
 * customers, employees, agents, transferees, successors, and assigns.
 *
 * The MPEG Software Simulation Group does not represent or warrant that the
 * programs furnished hereunder are free of infringement of any third-party
 * patents.
 *
 * Commercial implementations of MPEG-1 and MPEG-2 video, including shareware,
 * are subject to royalty fees to patent holders.  Many of these patents are
 * general enough such that they are unavoidable regardless of implementation
 * design.
 *
 */

#include <math.h>
#include "config.h"

#define W1 1.38703984532  /* sqrt(2)*cos(1*pi/16) */
#define W2 1.30656296487  /* sqrt(2)*cos(2*pi/16) */
#define W3 1.17587560242  /* sqrt(2)*cos(3*pi/16) */

#define W4 0.707106781187 /* cos(4*pi/16) */

#define W5 0.785694958387 /* sqrt(2)*cos(5*pi/16) */
#define W6 0.541196100153 /* sqrt(2)*cos(6*pi/16) */
#define W7 0.275899379291 /* sqrt(2)*cos(7*pi/16) */


/* global declarations */
void init_idct _ANSI_ARGS_((void));
void idct _ANSI_ARGS_((short *block));

/* private data */
static short iclip[1024]; /* clipping table */
static short *iclp;
static double temp[64];

/* private prototypes */
static void idctrow _ANSI_ARGS_((short *src, double *dst));
static void idctcol _ANSI_ARGS_((double *src, short *dst));

/* row (horizontal) IDCT
 *
 *           7                       pi         1
 * dst[k] = sum c[l] * src[l] * cos( -- * ( k + - ) * l )
 *          l=0                      8          2
 *
 * where: c[0]    = 128
 *        c[1..7] = 128*sqrt(2)
 */

static void idctrow(src,dst)
short *src;
double *dst;
{
  double x0, x1, x2, x3, x4, x5, x6, x7, x8;

  x0 = src[0];
  x1 = src[4];
  x2 = src[6];
  x3 = src[2];
  x4 = src[1];
  x5 = src[7];
  x6 = src[5];
  x7 = src[3];

  /* first stage */
  x8 = W7*(x4+x5);
  x4 = x8 + (W1-W7)*x4;
  x5 = x8 - (W1+W7)*x5;
  x8 = W3*(x6+x7);
  x6 = x8 - (W3-W5)*x6;
  x7 = x8 - (W3+W5)*x7;
  
  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = W6*(x3+x2);
  x2 = x1 - (W2+W6)*x2;
  x3 = x1 + (W2-W6)*x3;
  x1 = x4 + x6;
  x4 -= x6;
  x6 = x5 + x7;
  x5 -= x7;

  /* third stage */
  x7 = x8 + x3;
  x8 -= x3;
  x3 = x0 + x2;
  x0 -= x2;
  x2 = W4*(x4+x5);
  x4 = W4*(x4-x5);

  /* fourth stage */
  dst[0] = x7 + x1;
  dst[1] = x3 + x2;
  dst[2] = x0 + x4;
  dst[3] = x8 + x6;
  dst[4] = x8 - x6;
  dst[5] = x0 - x4;
  dst[6] = x3 - x2;
  dst[7] = x7 - x1;
}

/* column (vertical) IDCT
 *
 *             7                         pi         1
 * dst[8*k] = sum c[l] * src[8*l] * cos( -- * ( k + - ) * l )
 *            l=0                        8          2
 *
 * where: c[0]    = 1/1024
 *        c[1..7] = (1/1024)*sqrt(2)
 */
static void idctcol(src,dst)
double *src;
short *dst;
{
  double x0, x1, x2, x3, x4, x5, x6, x7, x8;

  x0 = src[8*0];
  x1 = src[8*4];
  x2 = src[8*6];
  x3 = src[8*2];
  x4 = src[8*1];
  x5 = src[8*7];
  x6 = src[8*5];
  x7 = src[8*3];

  /* first stage */
  x8 = W7*(x4+x5);
  x4 = x8 + (W1-W7)*x4;
  x5 = x8 - (W1+W7)*x5;
  x8 = W3*(x6+x7);
  x6 = x8 - (W3-W5)*x6;
  x7 = x8 - (W3+W5)*x7;

  /* second stage */
  x8 = x0 + x1;
  x0 -= x1;
  x1 = W6*(x3+x2);
  x2 = x1 - (W2+W6)*x2;
  x3 = x1 + (W2-W6)*x3;
  x1 = x4 + x6;
  x4 -= x6;
  x6 = x5 + x7;
  x5 -= x7;

  /* third stage */
  x7 = x8 + x3;
  x8 -= x3;
  x3 = x0 + x2;
  x0 -= x2;
  x2 = W4*(x4+x5);
  x4 = W4*(x4-x5);

  /* fourth stage */
  dst[8*0] = iclp[(int)floor(0.125*(x7+x1)+0.5)];
  dst[8*1] = iclp[(int)floor(0.125*(x3+x2)+0.5)];
  dst[8*2] = iclp[(int)floor(0.125*(x0+x4)+0.5)];
  dst[8*3] = iclp[(int)floor(0.125*(x8+x6)+0.5)];
  dst[8*4] = iclp[(int)floor(0.125*(x8-x6)+0.5)];
  dst[8*5] = iclp[(int)floor(0.125*(x0-x4)+0.5)];
  dst[8*6] = iclp[(int)floor(0.125*(x3-x2)+0.5)];
  dst[8*7] = iclp[(int)floor(0.125*(x7-x1)+0.5)];
}

/* two dimensional inverse discrete cosine transform */
void idct(block)
short *block;
{
  int i;

  for (i=0; i<8; i++)
    idctrow(block+8*i,temp+8*i);

  for (i=0; i<8; i++)
    idctcol(temp+i,block+i);
}

void init_idct()
{
  int i;

  iclp = iclip+512;
  for (i= -512; i<512; i++)
    iclp[i] = (i<-256) ? -256 : ((i>255) ? 255 : i);
}