/********************************************************************
 *                                                                  *
 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
 *                                                                  *
 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
 * by the XIPHOPHORUS Company http://www.xiph.org/                  *
 *                                                                  *
 ********************************************************************

 function: normalized modified discrete cosine transform
           power of two length transform only [64 <= n ]
 last mod: $Id: mdct.c,v 1.28 2001/12/21 14:52:35 segher Exp $

 Original algorithm adapted long ago from _The use of multirate filter
 banks for coding of high quality digital audio_, by T. Sporer,
 K. Brandenburg and B. Edler, collection of the European Signal
 Processing Conference (EUSIPCO), Amsterdam, June 1992, Vol.1, pp
 211-214

 The below code implements an algorithm that no longer looks much like
 that presented in the paper, but the basic structure remains if you
 dig deep enough to see it.

 This module DOES NOT INCLUDE code to generate/apply the window
 function.  Everybody has their own weird favorite including me... I
 happen to like the properties of y=sin(2PI*sin^2(x)), but others may
 vehemently disagree.

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

/* this can also be run as an integer transform by uncommenting a
   define in mdct.h; the integerization is a first pass and although
   it's likely stable for Vorbis, the dynamic range is constrained and
   roundoff isn't done (so it's noisy).  Consider it functional, but
   only a starting point.  There's no point on a machine with an FPU */

#include <stdio.h>
#include <string.h>
#include "mdct.h"
#include "mdct_lookup.h"
#include "os.h"
#include "Log.h"

/* build lookups for trig functions; also pre-figure scaling and
   some window function algebra. */

void mdct_init(mdct_lookup *lookup,int n){
  lookup->n=n;
  switch(n){
  case 64:
    lookup->log2n=6;
    lookup->trig =(int32_t*)triglook_64;
    break;
  case 128:
    lookup->log2n=7;
    lookup->trig=(int32_t*)triglook_128;
    break;
  case 256:
    lookup->log2n=8;
    lookup->trig=(int32_t*)triglook_256;
    break;
  case 512:
    lookup->log2n=9;
    lookup->trig=(int32_t*)triglook_512;
    break;
  case 1024:
    lookup->log2n=10;
    lookup->trig=(int32_t*)triglook_1024;
    break;
  case 2048:
    lookup->log2n=11;
    lookup->trig=(int32_t*)triglook_2048;
    break;
  case 4096:
    lookup->log2n=12;
    lookup->trig=(int32_t*)triglook_4096;
    break;
  case 8192:
    lookup->log2n=13;
    lookup->trig=(int32_t*)triglook_8192;
    break;
  default:
    /* die horribly */
    memset(lookup,0,sizeof(*lookup));
  }
  lookup->trig2 = lookup->trig + (n>>3) + 2;
  lookup->trig3 = lookup->trig2 + (n>>1);
}

/* 8 point butterfly (in place, 4 register) */
extern void mdct_butterfly_32(DATA_TYPE *x);
// N/stage point generic N stage butterfly (in place, 2 register)
extern void mdct_butterfly_generic(DATA_TYPE *T, DATA_TYPE *x,
					  int points, int trigint);
static void mdct_butterflies(mdct_lookup *init,
			     DATA_TYPE *x,
			     int points){

  DATA_TYPE *T=init->trig;
  int stages=init->log2n-5;
  int i,j;

  if(--stages>0){
    mdct_butterfly_generic(T,x,points,4);
  }

  for(i=1;--stages>0;i++){
    for(j=0;j<(1<<i);j++)
      mdct_butterfly_generic(T,x+(points>>i)*j,points>>i,4<<i);
  }

  for(j=0;j<points;j+=32)
    mdct_butterfly_32(x+j);
}

void mdct_clear(mdct_lookup *l){
 IGNORE(l);
}

extern void mdct_bit(DATA_TYPE *x, DATA_TYPE *T3, int log, int n);
/*
static void mdct_bitreverse(mdct_lookup *init,
			    DATA_TYPE *x){
  int        n       = init->n;
  int       *bit     = init->bitrev;
  DATA_TYPE *w0      = x;
  DATA_TYPE *w1      = x = w0+(n>>1);
  DATA_TYPE *T       = init->trig3;

  do{
    DATA_TYPE *x0    = x+*bit++;
    DATA_TYPE *x1    = x+*bit++;

    REG_TYPE  r0     = x0[1]  - x1[1];
    REG_TYPE  r1     = x0[0]  + x1[0];
    REG_TYPE  r2     = MULT_NORM(r1, T[0]) + MULT_NORM(r0, T[1]);
    REG_TYPE  r3     = MULT_NORM(r1, T[1]) - MULT_NORM(r0, T[0]);

	      w1    -= 2;

              r0     = HALVE(x0[1] + x1[1]);
              r1     = HALVE(x0[0] - x1[0]);

	      w0[0]  = r0     + r2;
	      w1[0]  = r0     - r2;
	      w0[1]  = r1     + r3;
	      w1[1]  = r3     - r1;

	      T     += 2;
	      w0    += 2;

  }while(w0<w1);
}
*/
extern void mdct_b0(DATA_TYPE *out, DATA_TYPE *in, int n, DATA_TYPE *trig);
extern void mdct_b1(DATA_TYPE* out, int n2, int n4, DATA_TYPE *T2);
/*
static void mdct_b1(DATA_TYPE* out, int n2, int n4, DATA_TYPE *T)
{
    DATA_TYPE *oX1=out+n2+n4;
    DATA_TYPE *oX2=out+n2+n4;
    DATA_TYPE *iX =out;

    do{
      oX1-=4;

      oX1[3]  =  MULT_NORM (iX[0], T[1]) - MULT_NORM(iX[1], T[0]);
      oX2[0]  = -MULT_NORM (iX[0], T[0]) - MULT_NORM(iX[1], T[1]);

      oX1[2]  =  MULT_NORM (iX[2], T[3]) - MULT_NORM(iX[3], T[2]);
      oX2[1]  = -MULT_NORM (iX[2], T[2]) - MULT_NORM(iX[3], T[3]);

      oX1[1]  =  MULT_NORM (iX[4], T[5]) - MULT_NORM(iX[5], T[4]);
      oX2[2]  = -MULT_NORM (iX[4], T[4]) - MULT_NORM(iX[5], T[5]);

      oX1[0]  =  MULT_NORM (iX[6], T[7]) - MULT_NORM(iX[7], T[6]);
      oX2[3]  = -MULT_NORM (iX[6], T[6]) - MULT_NORM(iX[7], T[7]);

      oX2+=4;
      iX    +=   8;
      T     +=   8;
    }while(iX<oX1);

    iX=out+n2+n4;
    oX1=out+n4;
    oX2=oX1;

    do{
      oX1-=4;
      iX-=4;

      oX2[0] = -(oX1[3] = iX[3]);
      oX2[1] = -(oX1[2] = iX[2]);
      oX2[2] = -(oX1[1] = iX[1]);
      oX2[3] = -(oX1[0] = iX[0]);

      oX2+=4;
    }while(oX2<iX);

    iX=out+n2+n4;
    oX1=out+n2+n4;
    oX2=out+n2;
    do{
      oX1-=4;
      oX1[0]= iX[3];
      oX1[1]= iX[2];
      oX1[2]= iX[1];
      oX1[3]= iX[0];
      iX+=4;
    }while(oX1>oX2);
}
*/

void mdct_backward(mdct_lookup *init, DATA_TYPE *in, DATA_TYPE *out){
  int n2=init->n>>1;
  int n4=init->n>>2;

  // rotate
  mdct_b0(out, in, init->n, init->trig);

  mdct_butterflies(init,out+n2,n2);
//  mdct_bitreverse(init,out);
  mdct_bit(out, init->trig3, init->log2n, init->n);

  // rotate + window
  mdct_b1(out, n2, n4, init->trig2);
}
