/* ***** BEGIN LICENSE BLOCK *****  
 * Source last modified: $Id: sbrfft.c,v 1.1 2005/02/26 01:47:35 jrecker Exp $ 
 *   
 * Portions Copyright (c) 1995-2005 RealNetworks, Inc. All Rights Reserved.  
 *       
 * The contents of this file, and the files included with this file, 
 * are subject to the current version of the RealNetworks Public 
 * Source License (the "RPSL") available at 
 * http://www.helixcommunity.org/content/rpsl unless you have licensed 
 * the file under the current version of the RealNetworks Community 
 * Source License (the "RCSL") available at 
 * http://www.helixcommunity.org/content/rcsl, in which case the RCSL 
 * will apply. You may also obtain the license terms directly from 
 * RealNetworks.  You may not use this file except in compliance with 
 * the RPSL or, if you have a valid RCSL with RealNetworks applicable 
 * to this file, the RCSL.  Please see the applicable RPSL or RCSL for 
 * the rights, obligations and limitations governing use of the 
 * contents of the file. 
 *   
 * This file is part of the Helix DNA Technology. RealNetworks is the 
 * developer of the Original Code and owns the copyrights in the 
 * portions it created. 
 *   
 * This file, and the files included with this file, is distributed 
 * and made available on an 'AS IS' basis, WITHOUT WARRANTY OF ANY 
 * KIND, EITHER EXPRESS OR IMPLIED, AND REALNETWORKS HEREBY DISCLAIMS 
 * ALL SUCH WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES 
 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, QUIET 
 * ENJOYMENT OR NON-INFRINGEMENT. 
 *  
 * Technology Compatibility Kit Test Suite(s) Location:  
 *    http://www.helixcommunity.org/content/tck  
 *  
 * Contributor(s):  
 *   
 * ***** END LICENSE BLOCK ***** */  

/**************************************************************************************
 * Fixed-point HE-AAC decoder
 * Jon Recker (jrecker@real.com), Ken Cooke (kenc@real.com)
 * February 2005
 *
 * sbrfft.c - optimized FFT for SBR QMF filters
 **************************************************************************************/

#include "sbr.h"
#include "assembly.h"

#define SQRT1_2	0x5a82799a

/* swap RE{p0} with RE{p1} and IM{P0} with IM{P1} */
#define swapcplx(p0,p1) \
	t = p0; t1 = *(&(p0)+1); p0 = p1; *(&(p0)+1) = *(&(p1)+1); p1 = t; *(&(p1)+1) = t1

/* nfft = 32, hard coded since small, fixed size FFT
static const unsigned char bitrevtab32[9] = {
	0x01, 0x04, 0x03, 0x06, 0x00, 0x02, 0x05, 0x07, 0x00,
};
*/

/* twiddle table for radix 4 pass, format = Q31 */
static const int twidTabOdd32[8*6] = {
	0x40000000, 0x00000000, 0x40000000, 0x00000000, 0x40000000, 0x00000000, 0x539eba45, 0xe7821d59, 
	0x4b418bbe, 0xf383a3e2, 0x58c542c5, 0xdc71898d, 0x5a82799a, 0xd2bec333, 0x539eba45, 0xe7821d59, 
	0x539eba45, 0xc4df2862, 0x539eba45, 0xc4df2862, 0x58c542c5, 0xdc71898d, 0x3248d382, 0xc13ad060, 
	0x40000000, 0xc0000000, 0x5a82799a, 0xd2bec333, 0x00000000, 0xd2bec333, 0x22a2f4f8, 0xc4df2862, 
	0x58c542c5, 0xcac933ae, 0xcdb72c7e, 0xf383a3e2, 0x00000000, 0xd2bec333, 0x539eba45, 0xc4df2862, 
	0xac6145bb, 0x187de2a7, 0xdd5d0b08, 0xe7821d59, 0x4b418bbe, 0xc13ad060, 0xa73abd3b, 0x3536cc52, 
};

/**************************************************************************************
 * Function:    BitReverse32
 *
 * Description: Ken's fast in-place bit reverse
 *
 * Inputs:      buffer of 32 complex samples
 *
 * Outputs:     bit-reversed samples in same buffer
 *
 * Return:      none
**************************************************************************************/
static void BitReverse32(int *inout)
{
	int t, t1;

	swapcplx(inout[2],  inout[32]);
	swapcplx(inout[4],  inout[16]);
	swapcplx(inout[6],  inout[48]);
	swapcplx(inout[10], inout[40]);
	swapcplx(inout[12], inout[24]);
	swapcplx(inout[14], inout[56]);
	swapcplx(inout[18], inout[36]);
	swapcplx(inout[22], inout[52]);
	swapcplx(inout[26], inout[44]);
	swapcplx(inout[30], inout[60]);
	swapcplx(inout[38], inout[50]);
	swapcplx(inout[46], inout[58]);
}

/**************************************************************************************
 * Function:    R8FirstPass32
 *
 * Description: radix-8 trivial pass for decimation-in-time FFT (log2(N) = 5)
 *
 * Inputs:      buffer of (bit-reversed) samples
 *
 * Outputs:     processed samples in same buffer
 *
 * Return:      none
 *
 * Notes:       assumes 3 guard bits, gains 1 integer bit
 *              guard bits out = guard bits in - 3 (if inputs are full scale)
 *                or guard bits in - 2 (if inputs bounded to +/- sqrt(2)/2)
 *              see scaling comments in fft.c for base AAC
 *              should compile with no stack spills on ARM (verify compiled output)
 *              current instruction count (per pass): 16 LDR, 16 STR, 4 SMULL, 61 ALU
 **************************************************************************************/
static void R8FirstPass32(int *r0)
{
	int r1, r2, r3, r4, r5, r6, r7;
	int r8, r9, r10, r11, r12, r14;

	/* number of passes = fft size / 8 = 32 / 8 = 4 */
	r1 = (32 >> 3);
	do {

		r2 = r0[8];
		r3 = r0[9];
		r4 = r0[10];
		r5 = r0[11];
		r6 = r0[12];
		r7 = r0[13];
		r8 = r0[14];
		r9 = r0[15];

		r10 = r2 + r4;
		r11 = r3 + r5;
		r12 = r6 + r8;
		r14 = r7 + r9;

		r2 -= r4;
		r3 -= r5;
		r6 -= r8;
		r7 -= r9;

		r4 = r2 - r7;
		r5 = r2 + r7;
		r8 = r3 - r6;
		r9 = r3 + r6;

		r2 = r4 - r9;
		r3 = r4 + r9;
		r6 = r5 - r8;
		r7 = r5 + r8;

		r2 = MULSHIFT32(SQRT1_2, r2);	/* can use r4, r5, r8, or r9 for constant and lo32 scratch reg */
		r3 = MULSHIFT32(SQRT1_2, r3);
		r6 = MULSHIFT32(SQRT1_2, r6);
		r7 = MULSHIFT32(SQRT1_2, r7);

		r4 = r10 + r12;
		r5 = r10 - r12;
		r8 = r11 + r14;
		r9 = r11 - r14;
	
		r10 = r0[0];
		r11 = r0[2];
		r12 = r0[4];
		r14 = r0[6];

		r10 += r11;
		r12 += r14;

		r4 >>= 1;
		r10 += r12;
		r4 += (r10 >> 1);
		r0[ 0] = r4;
		r4 -= (r10 >> 1);
		r4 = (r10 >> 1) - r4;
		r0[ 8] = r4;

		r9 >>= 1;
		r10 -= 2*r12;
		r4 = (r10 >> 1) + r9;
		r0[ 4] = r4;
		r4 = (r10 >> 1) - r9;
		r0[12] = r4;
		r10 += r12;

		r10 -= 2*r11;
		r12 -= 2*r14;

		r4 =  r0[1];
		r9 =  r0[3];
		r11 = r0[5];
		r14 = r0[7];

		r4 += r9;
		r11 += r14;

		r8 >>= 1;
		r4 += r11;
		r8 += (r4 >> 1);
		r0[ 1] = r8;
		r8 -= (r4 >> 1);
		r8 = (r4 >> 1) - r8;
		r0[ 9] = r8;

		r5 >>= 1;
		r4 -= 2*r11;
		r8 = (r4 >> 1) - r5;
		r0[ 5] = r8;
		r8 = (r4 >> 1) + r5;
		r0[13] = r8;
		r4 += r11;

		r4 -= 2*r9;
		r11 -= 2*r14;

		r9 = r10 - r11;
		r10 += r11;
		r14 = r4 + r12;
		r4 -= r12;

		r5 = (r10 >> 1) + r7;
		r8 = (r4 >> 1) - r6;
		r0[ 2] = r5;
		r0[ 3] = r8;

		r5 = (r9 >> 1) - r2;
		r8 = (r14 >> 1) - r3;
		r0[ 6] = r5;
		r0[ 7] = r8;

		r5 = (r10 >> 1) - r7;
		r8 = (r4 >> 1) + r6;
		r0[10] = r5;
		r0[11] = r8;
		
		r5 = (r9 >> 1) + r2;
		r8 = (r14 >> 1) + r3;
		r0[14] = r5;
		r0[15] = r8;

		r0 += 16;
		r1--;
	} while (r1 != 0);
}

/**************************************************************************************
 * Function:    R4Core32
 *
 * Description: radix-4 pass for 32-point decimation-in-time FFT
 *
 * Inputs:      buffer of samples
 *
 * Outputs:     processed samples in same buffer
 *
 * Return:      none
 *
 * Notes:       gain 2 integer bits
 *              guard bits out = guard bits in - 1 (if inputs are full scale)
 *              see scaling comments in fft.c for base AAC
 *              uses 3-mul, 3-add butterflies instead of 4-mul, 2-add
 *              should compile with no stack spills on ARM (verify compiled output)
 *              current instruction count (per pass): 16 LDR, 16 STR, 4 SMULL, 61 ALU
 **************************************************************************************/
static void R4Core32(int *r0)
{
	int r2, r3, r4, r5, r6, r7;
	int r8, r9, r10, r12, r14;
	int *r1;

	r1 = (int *)twidTabOdd32;
	r10 = 8;
	do {
		/* can use r14 for lo32 scratch register in all MULSHIFT32 */
		r2 = r1[0];
		r3 = r1[1];
		r4 = r0[16];
		r5 = r0[17];
		r12 = r4 + r5;
		r12 = MULSHIFT32(r3, r12);
		r5  = MULSHIFT32(r2, r5) + r12;
		r2 += 2*r3;
		r4  = MULSHIFT32(r2, r4) - r12;	

		r2 = r1[2];
		r3 = r1[3];
		r6 = r0[32];
		r7 = r0[33];
		r12 = r6 + r7;
		r12 = MULSHIFT32(r3, r12);
		r7  = MULSHIFT32(r2, r7) + r12;
		r2 += 2*r3;
		r6  = MULSHIFT32(r2, r6) - r12;
		
		r2 = r1[4];
		r3 = r1[5];
		r8 = r0[48];
		r9 = r0[49];
		r12 = r8 + r9;
		r12 = MULSHIFT32(r3, r12);
		r9  = MULSHIFT32(r2, r9) + r12;
		r2 += 2*r3;
		r8  = MULSHIFT32(r2, r8) - r12;

		r2 = r0[0];
		r3 = r0[1];

		r12 = r6 + r8;
		r8  = r6 - r8;
		r14 = r9 - r7;
		r9  = r9 + r7;

		r6 = (r2 >> 2) - r4;
		r7 = (r3 >> 2) - r5;
		r4 += (r2 >> 2);
		r5 += (r3 >> 2);

		r2 = r4 + r12;
		r3 = r5 + r9;
		r0[0] = r2;
		r0[1] = r3;
		r2 = r6 - r14;
		r3 = r7 - r8;
		r0[16] = r2;
		r0[17] = r3;
		r2 = r4 - r12;
		r3 = r5 - r9;
		r0[32] = r2;
		r0[33] = r3;
		r2 = r6 + r14;
		r3 = r7 + r8;
		r0[48] = r2;
		r0[49] = r3;

		r0 += 2;
		r1 += 6;
		r10--;
	} while (r10 != 0);
}

/**************************************************************************************
 * Function:    FFT32C
 *
 * Description: Ken's very fast in-place radix-4 decimation-in-time FFT
 *
 * Inputs:      buffer of 32 complex samples (before bit-reversal)
 *
 * Outputs:     processed samples in same buffer
 *
 * Return:      none
 *
 * Notes:       assumes 3 guard bits in, gains 3 integer bits
 *              guard bits out = guard bits in - 2
 *              (guard bit analysis includes assumptions about steps immediately
 *               before and after, i.e. PreMul and PostMul for DCT)
 **************************************************************************************/
void FFT32C(int *x)
{
	/* decimation in time */
	BitReverse32(x);

	/* 32-point complex FFT */
	R8FirstPass32(x);	/* gain 1 int bit,  lose 2 GB (making assumptions about input) */
	R4Core32(x);		/* gain 2 int bits, lose 0 GB (making assumptions about input) */
}