Serial-Studio/libs/QRealFourier/fftreal/FFTRealFixLen.hpp

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

        FFTRealFixLen.hpp
        By Laurent de Soras

--- Legal stuff ---

This program is free software. It comes without any warranty, to
the extent permitted by applicable law. You can redistribute it
and/or modify it under the terms of the Do What The Fuck You Want
To Public License, Version 2, as published by Sam Hocevar. See
http://sam.zoy.org/wtfpl/COPYING for more details.

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



#if defined (ffft_FFTRealFixLen_CURRENT_CODEHEADER)
	#error Recursive inclusion of FFTRealFixLen code header.
#endif
#define	ffft_FFTRealFixLen_CURRENT_CODEHEADER

#if ! defined (ffft_FFTRealFixLen_CODEHEADER_INCLUDED)
#define	ffft_FFTRealFixLen_CODEHEADER_INCLUDED



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

#include	"def.h"
#include	"FFTRealPassDirect.h"
#include	"FFTRealPassInverse.h"
#include	"FFTRealSelect.h"

#include	<cassert>
#include	<cmath>

namespace std { }



namespace ffft
{



/*\\\ 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 <>
inline 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 <>
inline 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 <>
inline 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 <>
inline 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 <>
inline 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 <>
inline 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);
	}
}



}	// namespace ffft



#endif	// ffft_FFTRealFixLen_CODEHEADER_INCLUDED

#undef ffft_FFTRealFixLen_CURRENT_CODEHEADER



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