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

        FFTRealFixLen.hpp

        Copyright (c) 2005 Laurent de Soras

--- Legal stuff ---

This library is free software; you can redistribute it and/or

modify it under the terms of the GNU Lesser General Public

License as published by the Free Software Foundation; either

version 2.1 of the License, or (at your option) any later version.

This library 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

Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public

License along with this library; if not, write to the Free Software

Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

*Tab=3***********************************************************************/

#if defined (FFTRealFixLen_CURRENT_CODEHEADER)

	#error Recursive inclusion of FFTRealFixLen code header.

#endif

#define	FFTRealFixLen_CURRENT_CODEHEADER

#if ! defined (FFTRealFixLen_CODEHEADER_INCLUDED)

#define	FFTRealFixLen_CODEHEADER_INCLUDED

/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

#include	"def.h"

#include	"FFTRealPassDirect.h"

#include	"FFTRealPassInverse.h"

#include	"FFTRealSelect.h"

#include	<cassert>

#include	<cmath>

namespace std { }

/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

template <int LL2>

FFTRealFixLen <LL2>::FFTRealFixLen ()

:	_buffer (FFT_LEN)

,	_br_data (BR_ARR_SIZE)

,	_trigo_data (TRIGO_TABLE_ARR_SIZE)

,	_trigo_osc ()

{

	build_br_lut ();

	build_trigo_lut ();

	build_trigo_osc ();

}

template <int LL2>

long	FFTRealFixLen <LL2>::get_length () const

{

	return (FFT_LEN);

}

// General case

template <int LL2>

void	FFTRealFixLen <LL2>::do_fft (DataType f [], const DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	assert (x != f);

	assert (FFT_LEN_L2 >= 3);

	// Do the transform in several passes

	const DataType	*	cos_ptr = &_trigo_data [0];

	const long *	br_ptr = &_br_data [0];

	FFTRealPassDirect <FFT_LEN_L2 - 1>::process (

		FFT_LEN,

		f,

		&_buffer [0],

		x,

		cos_ptr,

		TRIGO_TABLE_ARR_SIZE,

		br_ptr,

		&_trigo_osc [0]

	);

}

// 4-point FFT

template <>

void	FFTRealFixLen <2>::do_fft (DataType f [], const DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	assert (x != f);

	f [1] = x [0] - x [2];

	f [3] = x [1] - x [3];

	const DataType	b_0 = x [0] + x [2];

	const DataType	b_2 = x [1] + x [3];

	f [0] = b_0 + b_2;

	f [2] = b_0 - b_2;

}

// 2-point FFT

template <>

void	FFTRealFixLen <1>::do_fft (DataType f [], const DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	assert (x != f);

	f [0] = x [0] + x [1];

	f [1] = x [0] - x [1];

}

// 1-point FFT

template <>

void	FFTRealFixLen <0>::do_fft (DataType f [], const DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	f [0] = x [0];

}

// General case

template <int LL2>

void	FFTRealFixLen <LL2>::do_ifft (const DataType f [], DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	assert (x != f);

	assert (FFT_LEN_L2 >= 3);

	// Do the transform in several passes

	DataType *		s_ptr =

		FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (&_buffer [0], x);

	DataType *		d_ptr =

		FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (x, &_buffer [0]);

	const DataType	*	cos_ptr = &_trigo_data [0];

	const long *	br_ptr = &_br_data [0];

	FFTRealPassInverse <FFT_LEN_L2 - 1>::process (

		FFT_LEN,

		d_ptr,

		s_ptr,

		f,

		cos_ptr,

		TRIGO_TABLE_ARR_SIZE,

		br_ptr,

		&_trigo_osc [0]

	);

}

// 4-point IFFT

template <>

void	FFTRealFixLen <2>::do_ifft (const DataType f [], DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	assert (x != f);

	const DataType	b_0 = f [0] + f [2];

	const DataType	b_2 = f [0] - f [2];

	x [0] = b_0 + f [1] * 2;

	x [2] = b_0 - f [1] * 2;

	x [1] = b_2 + f [3] * 2;

	x [3] = b_2 - f [3] * 2;

}

// 2-point IFFT

template <>

void	FFTRealFixLen <1>::do_ifft (const DataType f [], DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	assert (x != f);

	x [0] = f [0] + f [1];

	x [1] = f [0] - f [1];

}

// 1-point IFFT

template <>

void	FFTRealFixLen <0>::do_ifft (const DataType f [], DataType x [])

{

	assert (f != 0);

	assert (x != 0);

	assert (x != f);

	x [0] = f [0];

}

template <int LL2>

void	FFTRealFixLen <LL2>::rescale (DataType x []) const

{

	assert (x != 0);

	const DataType	mul = DataType (1.0 / FFT_LEN);

	if (FFT_LEN < 4)

	{

		long				i = FFT_LEN - 1;

		do

		{

			x [i] *= mul;

			--i;

		}

		while (i >= 0);

	}

	else

	{

		assert ((FFT_LEN & 3) == 0);

		// Could be optimized with SIMD instruction sets (needs alignment check)

		long				i = FFT_LEN - 4;

		do

		{

			x [i + 0] *= mul;

			x [i + 1] *= mul;

			x [i + 2] *= mul;

			x [i + 3] *= mul;

			i -= 4;

		}

		while (i >= 0);

	}

}

/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

template <int LL2>

void	FFTRealFixLen <LL2>::build_br_lut ()

{

	_br_data [0] = 0;

	for (long cnt = 1; cnt < BR_ARR_SIZE; ++cnt)

	{

		long				index = cnt << 2;

		long				br_index = 0;

		int				bit_cnt = FFT_LEN_L2;

		do

		{

			br_index <<= 1;

			br_index += (index & 1);

			index >>= 1;

			-- bit_cnt;

		}

		while (bit_cnt > 0);

		_br_data [cnt] = br_index;

	}

}

template <int LL2>

void	FFTRealFixLen <LL2>::build_trigo_lut ()

{

	const double	mul = (0.5 * PI) / TRIGO_TABLE_ARR_SIZE;

	for (long i = 0; i < TRIGO_TABLE_ARR_SIZE; ++ i)

	{

		using namespace std;

		_trigo_data [i] = DataType (cos (i * mul));

	}

}

template <int LL2>

void	FFTRealFixLen <LL2>::build_trigo_osc ()

{

	for (int i = 0; i < NBR_TRIGO_OSC; ++i)

	{

		OscType &		osc = _trigo_osc [i];

		const long		len = static_cast <long> (TRIGO_TABLE_ARR_SIZE) << (i + 1);

		const double	mul = (0.5 * PI) / len;

		osc.set_step (mul);

	}

}

#endif	// FFTRealFixLen_CODEHEADER_INCLUDED

#undef FFTRealFixLen_CURRENT_CODEHEADER

/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/