/* this file is a part of amp software, (C) tomislav uzelac 1996,1997
*/
 
/* transform.c  imdct and polyphase(DCT) transforms
 *
 * Created by: tomislav uzelac  May 1996
 * Karl Anders Oygard optimized this for speed, Mar 13 97
 */
#include <math.h>
#include <sys/types.h>
#include <unistd.h>
#include <string.h>

#include "audio.h"
#include "getdata.h"
#include "misc2.h"

#define TRANSFORM
#include "transform.h"

static inline float pass(float *x,float c,int is_short)
{
float q1,q2,tmp;

	q1=q2=*x++;
	if (!is_short) {
	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp; 
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp;
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp;   
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp;
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp;
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp;
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	}

	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp;
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	tmp=*x++; q1+=c*q2; q2=-q2; q1+=tmp; q2+=tmp;
	tmp=*x++; q2+=c*q1; q1=-q1; q2+=tmp; q1+=tmp;
	tmp=*x;   q1+=c*q2;         q1+=tmp;

	return q1;
}

/* imdct is done according to Chiang/Liu[96]
 */
void imdct(int win_type,int sb,int ch)
{
float *x,c,q;
int i,j;

	x=&xr[ch][18*sb];
	if (win_type!=2) {
		j=17;
		for (i=0;i<9;i++) {
			c=t_2cos[win_type][i];
			q=pass(x,c,0);  
			res[ch][sb][i]= q*t_sin[win_type][i] + s[ch][sb][i];
			res[ch][sb][j]= q*t_sin[win_type][j] + s[ch][sb][j];
			j--;
		}
		j=17;
		for (i=0;i<9;i++) {
			c=t_2cos[win_type][9+i]; /* t_2cos ima 18 clanova zbog interne simetrije */
			q=pass(x,c,0);
			s[ch][sb][i]= q * t_sin[win_type][18+i];
			s[ch][sb][j]= q * t_sin[win_type][18+j];
			j--;
		}
	} else {
		for (i=0;i<18;i++) res[ch][sb][i]=s[ch][sb][i];
		for (i=0;i<3;i++) {
			c=t_2cos[2][i]; q=pass(x,c,1);
			res[ch][sb][6+i]+=q*t_sin[2][i];
			res[ch][sb][11-i]+=q*t_sin[2][5-i];
		}
		for (i=0;i<3;i++) {
			c=t_2cos[2][3+i]; q=pass(x,c,1);
			res[ch][sb][12+i]+=q*t_sin[2][6+i];
			res[ch][sb][17-i]+=q*t_sin[2][11-i];
		}
		x+=6;
		for (i=0;i<3;i++) {
			c=t_2cos[2][i]; q=pass(x,c,1);
			res[ch][sb][12+i]+=q*t_sin[2][i];
			res[ch][sb][17-i]+=q*t_sin[2][5-i];
		}

		for (i=0;i<18;i++) s[ch][sb][i]=0;
		for (i=0;i<3;i++) {
			c=t_2cos[2][3+i]; q=pass(x,c,1);
			s[ch][sb][i]+=q*t_sin[2][6+i];
			s[ch][sb][5-i]+=q*t_sin[2][11-i];
		}
		x+=6;
		for (i=0;i<3;i++) {
			c=t_2cos[2][i]; q=pass(x,c,1);
			s[ch][sb][i]+=q*t_sin[2][i];
			s[ch][sb][5-i]+=q*t_sin[2][5-i];
		}
		for (i=0;i<3;i++) {
			c=t_2cos[2][3+i]; q=pass(x,c,1);
			s[ch][sb][6+i]+=q*t_sin[2][6+i];
			s[ch][sb][11-i]+=q*t_sin[2][11-i];
		}
	}

	if (sb&1) for (i=0;i<18;i++) if (i&1) res[ch][sb][i]=-res[ch][sb][i];

	/* if (sb&1) {
	 *      float *ptr=&res[ch][sb][1];
	 *      for (i=0;i<9;i++) {
	 *              *ptr=-*ptr;
	 *              ptr+=2;
	 *      }
	 * }
	 */
}

/* fast DCT according to Lee[84]
 * reordering according to Konstantinides[94]
 */ 
void poly(const int ch,int f)
{
static float u[2][2][16][32]; /* no v[][], it's redundant */
static int u_start=0; /* first element of u[][] */
static int u_div=0; /* which part of u[][] is currently used */
float *p,*q;
float c0,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15;
float c16,c17,c18,c19,c20,c21,c22,c23,c24,c25,c26,c27,c28,c29,c30,c31;
float d0,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14,d15;
float d16,d17,d18,d19,d20,d21,d22,d23,d24,d25,d26,d27,d28,d29,d30,d31;


	/* step 1: initial reordering and 1st (16 wide) butterflies
	*/

	d0 =res[ch][ 0][f]; d16=(d0  - res[ch][31][f]) *  b1; d0 += res[ch][31][f];
	d1 =res[ch][ 1][f]; d17=(d1  - res[ch][30][f]) *  b3; d1 += res[ch][30][f];
	d3 =res[ch][ 2][f]; d19=(d3  - res[ch][29][f]) *  b5; d3 += res[ch][29][f];
	d2 =res[ch][ 3][f]; d18=(d2  - res[ch][28][f]) *  b7; d2 += res[ch][28][f];
	d6 =res[ch][ 4][f]; d22=(d6  - res[ch][27][f]) *  b9; d6 += res[ch][27][f];
	d7 =res[ch][ 5][f]; d23=(d7  - res[ch][26][f]) * b11; d7 += res[ch][26][f];
	d5 =res[ch][ 6][f]; d21=(d5  - res[ch][25][f]) * b13; d5 += res[ch][25][f];
	d4 =res[ch][ 7][f]; d20=(d4  - res[ch][24][f]) * b15; d4 += res[ch][24][f];
	d12=res[ch][ 8][f]; d28=(d12 - res[ch][23][f]) * b17; d12+= res[ch][23][f];
	d13=res[ch][ 9][f]; d29=(d13 - res[ch][22][f]) * b19; d13+= res[ch][22][f];
	d15=res[ch][10][f]; d31=(d15 - res[ch][21][f]) * b21; d15+= res[ch][21][f];
	d14=res[ch][11][f]; d30=(d14 - res[ch][20][f]) * b23; d14+= res[ch][20][f];
	d10=res[ch][12][f]; d26=(d10 - res[ch][19][f]) * b25; d10+= res[ch][19][f];
	d11=res[ch][13][f]; d27=(d11 - res[ch][18][f]) * b27; d11+= res[ch][18][f];
	d9 =res[ch][14][f]; d25=(d9  - res[ch][17][f]) * b29; d9 += res[ch][17][f];
	d8 =res[ch][15][f]; d24=(d8  - res[ch][16][f]) * b31; d8 += res[ch][16][f];

	/* step 2: 8-wide butterflies
	*/
	c0 = d0 + d8 ; c8 = ( d0 - d8 ) *  b2;
	c1 = d1 + d9 ; c9 = ( d1 - d9 ) *  b6;
	c2 = d2 + d10; c10= ( d2 - d10) * b14;
	c3 = d3 + d11; c11= ( d3 - d11) * b10;
	c4 = d4 + d12; c12= ( d4 - d12) * b30;
	c5 = d5 + d13; c13= ( d5 - d13) * b26;
	c6 = d6 + d14; c14= ( d6 - d14) * b18;
	c7 = d7 + d15; c15= ( d7 - d15) * b22;
	
	c16=d16 + d24; c24= (d16 - d24) *  b2;
	c17=d17 + d25; c25= (d17 - d25) *  b6;
	c18=d18 + d26; c26= (d18 - d26) * b14;
	c19=d19 + d27; c27= (d19 - d27) * b10;
	c20=d20 + d28; c28= (d20 - d28) * b30;
	c21=d21 + d29; c29= (d21 - d29) * b26;
	c22=d22 + d30; c30= (d22 - d30) * b18;
	c23=d23 + d31; c31= (d23 - d31) * b22;
	
	/* step 3: 4-wide butterflies
	*/
	d0 = c0 + c4 ; d4 = ( c0 - c4 ) *  b4;
	d1 = c1 + c5 ; d5 = ( c1 - c5 ) * b12;
	d2 = c2 + c6 ; d6 = ( c2 - c6 ) * b28;
	d3 = c3 + c7 ; d7 = ( c3 - c7 ) * b20;
	
	d8 = c8 + c12; d12= ( c8 - c12) *  b4;
	d9 = c9 + c13; d13= ( c9 - c13) * b12;
	d10= c10+ c14; d14= (c10 - c14) * b28;
	d11= c11+ c15; d15= (c11 - c15) * b20;

	d16= c16+ c20; d20= (c16 - c20) *  b4;
	d17= c17+ c21; d21= (c17 - c21) * b12;
	d18= c18+ c22; d22= (c18 - c22) * b28;
	d19= c19+ c23; d23= (c19 - c23) * b20;
	
	d24= c24+ c28; d28= (c24 - c28) *  b4;
	d25= c25+ c29; d29= (c25 - c29) * b12;
	d26= c26+ c30; d30= (c26 - c30) * b28;
	d27= c27+ c31; d31= (c27 - c31) * b20;

	/* step 4: 2-wide butterflies
	*/
/**/    c0 = d0 + d2 ; c2 = ( d0 - d2 ) *  b8;
	c1 = d1 + d3 ; c3 = ( d1 - d3 ) * b24;
/**/    c4 = d4 + d6 ; c6 = ( d4 - d6 ) *  b8;
	c5 = d5 + d7 ; c7 = ( d5 - d7 ) * b24;
/**/    c8 = d8 + d10; c10= ( d8 - d10) *  b8;
	c9 = d9 + d11; c11= ( d9 - d11) * b24;
/**/    c12= d12+ d14; c14= (d12 - d14) *  b8;
	c13= d13+ d15; c15= (d13 - d15) * b24;
/**/    c16= d16+ d18; c18= (d16 - d18) *  b8;
	c17= d17+ d19; c19= (d17 - d19) * b24;  
/**/    c20= d20+ d22; c22= (d20 - d22) *  b8;
	c21= d21+ d23; c23= (d21 - d23) * b24;  
/**/    c24= d24+ d26; c26= (d24 - d26) *  b8;
	c25= d25+ d27; c27= (d25 - d27) * b24;  
/**/    c28= d28+ d30; c30= (d28 - d30) *  b8;
	c29= d29+ d31; c31= (d29 - d31) * b24;  
	
	/* step 5: 1-wide butterflies
	*/
	d0 = c0 + c1 ; d1 = ( c0 - c1 ) * b16;
	d2 = c2 + c3 ; d3 = ( c2 - c3 ) * b16;
	d4 = c4 + c5 ; d5 = ( c4 - c5 ) * b16;
	d6 = c6 + c7 ; d7 = ( c6 - c7 ) * b16;
	d8 = c8 + c9 ; d9 = ( c8 - c9 ) * b16;
	d10= c10+ c11; d11= (c10 - c11) * b16;
	d12= c12+ c13; d13= (c12 - c13) * b16;
	d14= c14+ c15; d15= (c14 - c15) * b16;
	d16= c16+ c17; d17= (c16 - c17) * b16;
	d18= c18+ c19; d19= (c18 - c19) * b16;
	d20= c20+ c21; d21= (c20 - c21) * b16;
	d22= c22+ c23; d23= (c22 - c23) * b16;
	d24= c24+ c25; d25= (c24 - c25) * b16;
	d26= c26+ c27; d27= (c26 - c27) * b16;
	d28= c28+ c29; d29= (c28 - c29) * b16;
	d30= c30+ c31; d31= (c30 - c31) * b16;
	
	/* step 6: final resolving & reordering 
	 * first 32 are used right away
	 */

	p=&u[ch][u_div][u_start][0]; 
	q=&u[ch][u_div][u_start][31];

/*16*/  *p++ =+d1 ;
	*q-- = -(*p++ =+d16 +d17 +d18 +d22 -d30);
	*q-- = -(*p++ =+d8 +d9 +d10 -d14);
	*q-- = -(*p++ =-d16 -d17 -d18 -d22 +d24 +d25 +d26);
/*20*/  *q-- = -(*p++ =+d4 +d5 -d6);
	*q-- = -(*p++ =+d16 +d17 +d18 +d20 +d21 -d24 -d25 -d26);
	*q-- = -(*p++ =-d8 -d9 -d10 +d12 +d13);
	*q-- = -(*p++ =-d16 -d17 -d18 -d20 -d21 +d28 +d29);
/*24*/  *q-- = -(*p++ =-d2 +d3);
	*q-- = -(*p++ =+d16 +d17 +d19 +d20 +d21 -d28 -d29);
	*q-- = -(*p++ =+d8 +d9 +d11 -d12 -d13);
	*q-- = -(*p++ =-d16 -d17 -d19 -d20 -d21 +d24 +d25 +d27);
/*28*/  *q-- = -(*p++ =-d4 -d5 +d7);
	*q-- = -(*p++ =+d16 +d17 +d19 +d23 -d24 -d25 -d27);
	*q-- = -(*p++ =-d8 -d9 -d11 +d15);
	*q-- = -(*p   =-d16 -d17 -d19 -d23 +d31);
	*q   = 0;

	/* the other 32 are stored for use with the next granule
	 */

	p=&u[ch][u_div ? 0 : 1][u_start][16];
	q=&u[ch][u_div ? 0 : 1][u_start][17];

/*0*/   *p-- = -2*d0;
	*p-- = *q++ = -(+d16 );
	*p-- = *q++ = -(+d8 );
	*p-- = *q++ = -(-d16 +d24 );
/*4*/   *p-- = *q++ = -(+d4 );
	*p-- = *q++ = -(+d16 +d20 -d24 );
	*p-- = *q++ = -(-d8 +d12 );
	*p-- = *q++ = -(-d16 -d20 +d28 );
/*8*/   *p-- = *q++ = -(+d2 );
	*p-- = *q++ = -(+d16 +d18 +d20 -d28 );
	*p-- = *q++ = -(+d8 +d10 -d12 );
	*p-- = *q++ = -(-d16 -d18 -d20 +d24 +d26 );
/*12*/  *p-- = *q++ = -(-d4 +d6 );
	*p-- = *q++ = -(+d16 +d18 +d22 -d24 -d26 );
	*p-- = *q++ = -(-d8 -d10 +d14 );
	*p-- = *q   = -(-d16 -d18 -d22 +d30 );  
	*p   = -d1;

	/* these are 2x bigger than they're supposed to, but we'll 
	 * compensate later for that.
	 */

	/* we're doing dewindowing and calculating final samples now
	 */

	{
	int i = u_start,j;
	short *putptr;
	float* u_ptr = (float*) u[ch][u_div];
	float out;

#ifdef MONO
		if (nch == 2 && ch == 1 && !AMIGA_MONO) putptr=&(pcm_sample[1]);
		else putptr=pcm_sample;
#else
		if (nch == 2 && ch == 1) putptr=&(pcm_sample[1]);
		else putptr=pcm_sample;
#endif

#ifdef FREQDIV
		for (j=0;j<32;j+=AMIGA_FDIV) {
#else
		for (j=0;j<32;j++) {
#endif
			out  = u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 0][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 1][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 2][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 3][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 4][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 5][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 6][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 7][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 8][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[ 9][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[10][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[11][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[12][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[13][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[14][j];
			out += u_ptr[((i++ & 0xf) << 5) + j] * t_dewindow[15][j];

			/* the ifs are, regrettably, necessary
			*/
			if (out > 2.0) *putptr=32767;
			else if (out < -2.0) *putptr=-32768;
			else *putptr=(short)(out*16383.5f);
#ifdef MONO
			if (nch == 2 && !AMIGA_MONO) putptr+=2;
#else
			if (nch == 2) putptr+=2;
#endif
#ifdef FREQDIV
			else putptr+=AMIGA_FDIV;
			i+=((AMIGA_FDIV-1)*16);
#else
			else putptr++;
#endif
		}
#ifdef MONO
		if (AMIGA_MONO) {
			u_start=(u_start-1)&0xf;
			u_div=u_div ? 0 : 1;
			}
		else
#endif
		if (ch || nch ==1) {
			u_start=(u_start-1)&0xf;
			u_div=u_div ? 0 : 1;
		}       
	}
}


/* little endian systems do not require any byte swapping whatsoever. 
 * big endian systems require byte swapping when writing .wav files, 
 * but not when playing directly
 */
void printout(void)
{
#ifdef FREQDIV
int i;
static short puttab[64]; /* a table to copy into */
static short *putptr;    /* a pointer to puttab  */
short *getptr;           /* pointer to pcm_sample */
#endif
if (A_WRITE_TO_FILE) {
#ifdef NO_BYTE_SWAPPING

#ifdef MONO
	if (nch!=1 && !AMIGA_MONO) fwrite(pcm_sample,2,64,out_file);
#else
	if (nch!=1) fwrite(pcm_sample,2,64,out_file);
#endif
	else
#ifndef FREQDIV
	    fwrite(pcm_sample,2,32,out_file);
#else
	if (AMIGA_FDIV==1)
	    fwrite(pcm_sample,2,32,out_file);
	else {
	getptr=pcm_sample;
	putptr=puttab;
	for (i=0;i<32;i++) {
	    *putptr=*getptr;
	    putptr++;
	    getptr+=AMIGA_FDIV;
	    }
	fwrite(puttab,2,AMIGA_FSIZE,out_file);
	}
#endif

#else

int i;
unsigned char tmp[128],*ptr;
	ptr=tmp;
	for (i=0;(nch == 1) ? i<32 : i<64;i++) {
		*ptr++=pcm_sample[i]&0xff;
		*ptr++=pcm_sample[i]>>8;
	}
	fwrite(tmp,1,(nch == 1) ? 64 : 128,out_file);

#endif
}

if (A_WRITE_TO_AUDIO) {
#ifdef A_SIF_HPUX
#define AUSIZ 4096
#else
#define AUSIZ 32768
#endif
static unsigned char au_buf[AUSIZ];
static int au_ptr=0;
	memcpy(au_buf+au_ptr,pcm_sample,(nch==1 ? 64:128));
	au_ptr+=(nch==1 ? 64:128);
	if (au_ptr==AUSIZ) {
		au_ptr=0;
#if defined(A_SIF_IRIX)
		ALwritesamps(audioport, au_buf, AUSIZ / 2);
#else /* A_SIF_IRIX */
		write(audio_fd,au_buf,AUSIZ);
#endif /* all others */
	}
}
}