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

        FFTRealPassDirect.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 (FFTRealPassDirect_CURRENT_CODEHEADER)

	#error Recursive inclusion of FFTRealPassDirect code header.

#endif

#define	FFTRealPassDirect_CURRENT_CODEHEADER

#if ! defined (FFTRealPassDirect_CODEHEADER_INCLUDED)

#define	FFTRealPassDirect_CODEHEADER_INCLUDED

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

#include	"FFTRealUseTrigo.h"

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

template <>

void	FFTRealPassDirect <1>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])

{

	// First and second pass at once

	const long		qlen = len >> 2;

	long				coef_index = 0;

	do

	{

		// To do: unroll the loop (2x).

		const long		ri_0 = br_ptr [coef_index >> 2];

		const long		ri_1 = ri_0 + 2 * qlen;	// bit_rev_lut_ptr [coef_index + 1];

		const long		ri_2 = ri_0 + 1 * qlen;	// bit_rev_lut_ptr [coef_index + 2];

		const long		ri_3 = ri_0 + 3 * qlen;	// bit_rev_lut_ptr [coef_index + 3];

		DataType	* const	df2 = dest_ptr + coef_index;

		df2 [1] = x_ptr [ri_0] - x_ptr [ri_1];

		df2 [3] = x_ptr [ri_2] - x_ptr [ri_3];

		const DataType	sf_0 = x_ptr [ri_0] + x_ptr [ri_1];

		const DataType	sf_2 = x_ptr [ri_2] + x_ptr [ri_3];

		df2 [0] = sf_0 + sf_2;

		df2 [2] = sf_0 - sf_2;

		coef_index += 4;

	}

	while (coef_index < len);

}

template <>

void	FFTRealPassDirect <2>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])

{

	// Executes "previous" passes first. Inverts source and destination buffers

	FFTRealPassDirect <1>::process (

		len,

		src_ptr,

		dest_ptr,

		x_ptr,

		cos_ptr,

		cos_len,

		br_ptr,

		osc_list

	);

	// Third pass

	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);

	long				coef_index = 0;

	do

	{

		dest_ptr [coef_index    ] = src_ptr [coef_index] + src_ptr [coef_index + 4];

		dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4];

		dest_ptr [coef_index + 2] = src_ptr [coef_index + 2];

		dest_ptr [coef_index + 6] = src_ptr [coef_index + 6];

		DataType			v;

		v = (src_ptr [coef_index + 5] - src_ptr [coef_index + 7]) * sqrt2_2;

		dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + v;

		dest_ptr [coef_index + 3] = src_ptr [coef_index + 1] - v;

		v = (src_ptr [coef_index + 5] + src_ptr [coef_index + 7]) * sqrt2_2;

		dest_ptr [coef_index + 5] = v + src_ptr [coef_index + 3];

		dest_ptr [coef_index + 7] = v - src_ptr [coef_index + 3];

		coef_index += 8;

	}

	while (coef_index < len);

}

template <int PASS>

void	FFTRealPassDirect <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType x_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list [])

{

	// Executes "previous" passes first. Inverts source and destination buffers

	FFTRealPassDirect <PASS - 1>::process (

		len,

		src_ptr,

		dest_ptr,

		x_ptr,

		cos_ptr,

		cos_len,

		br_ptr,

		osc_list

	);

	const long		dist = 1L << (PASS - 1);

	const long		c1_r = 0;

	const long		c1_i = dist;

	const long		c2_r = dist * 2;

	const long		c2_i = dist * 3;

	const long		cend = dist * 4;

	const long		table_step = cos_len >> (PASS - 1);

   enum {	TRIGO_OSC		= PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT	};

	enum {	TRIGO_DIRECT	= (TRIGO_OSC >= 0) ? 1 : 0	};

	long				coef_index = 0;

	do

	{

		const DataType	* const	sf = src_ptr + coef_index;

		DataType			* const	df = dest_ptr + coef_index;

		// Extreme coefficients are always real

		df [c1_r] = sf [c1_r] + sf [c2_r];

		df [c2_r] = sf [c1_r] - sf [c2_r];

		df [c1_i] = sf [c1_i];

		df [c2_i] = sf [c2_i];

		FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]);

		// Others are conjugate complex numbers

		for (long i = 1; i < dist; ++ i)

		{

			DataType			c;

			DataType			s;

			FFTRealUseTrigo <TRIGO_DIRECT>::iterate (

				osc_list [TRIGO_OSC],

				c,

				s,

				cos_ptr,

				i * table_step,

				(dist - i) * table_step

			);

			const DataType	sf_r_i = sf [c1_r + i];

			const DataType	sf_i_i = sf [c1_i + i];

			const DataType	v1 = sf [c2_r + i] * c - sf [c2_i + i] * s;

			df [c1_r + i] = sf_r_i + v1;

			df [c2_r - i] = sf_r_i - v1;

			const DataType	v2 = sf [c2_r + i] * s + sf [c2_i + i] * c;

			df [c2_r + i] = v2 + sf_i_i;

			df [cend - i] = v2 - sf_i_i;

		}

		coef_index += cend;

	}

	while (coef_index < len);

}

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

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

#endif	// FFTRealPassDirect_CODEHEADER_INCLUDED

#undef FFTRealPassDirect_CURRENT_CODEHEADER

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