fftreal.hpp

Real FFT , fftw compiled component

		sf = temp_ptr;
	}
}



template <class DT>
void	FFTReal <DT>::compute_direct_pass_1_2 (DataType df [], const DataType x []) const
{
	assert (df != 0);
	assert (x != 0);
	assert (df != x);

	const long * const	bit_rev_lut_ptr = get_br_ptr ();
	long				coef_index = 0;
	do
	{
		const long		rev_index_0 = bit_rev_lut_ptr [coef_index];
		const long		rev_index_1 = bit_rev_lut_ptr [coef_index + 1];
		const long		rev_index_2 = bit_rev_lut_ptr [coef_index + 2];
		const long		rev_index_3 = bit_rev_lut_ptr [coef_index + 3];

		DataType	* const	df2 = df + coef_index;
		df2 [1] = x [rev_index_0] - x [rev_index_1];
		df2 [3] = x [rev_index_2] - x [rev_index_3];

		const DataType	sf_0 = x [rev_index_0] + x [rev_index_1];
		const DataType	sf_2 = x [rev_index_2] + x [rev_index_3];

		df2 [0] = sf_0 + sf_2;
		df2 [2] = sf_0 - sf_2;
		
		coef_index += 4;
	}
	while (coef_index < _length);
}



template <class DT>
void	FFTReal <DT>::compute_direct_pass_3 (DataType df [], const DataType sf []) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);

	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);
	long				coef_index = 0;
	do
	{
		DataType			v;

		df [coef_index] = sf [coef_index] + sf [coef_index + 4];
		df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
		df [coef_index + 2] = sf [coef_index + 2];
		df [coef_index + 6] = sf [coef_index + 6];

		v = (sf [coef_index + 5] - sf [coef_index + 7]) * sqrt2_2;
		df [coef_index + 1] = sf [coef_index + 1] + v;
		df [coef_index + 3] = sf [coef_index + 1] - v;

		v = (sf [coef_index + 5] + sf [coef_index + 7]) * sqrt2_2;
		df [coef_index + 5] = v + sf [coef_index + 3];
		df [coef_index + 7] = v - sf [coef_index + 3];

		coef_index += 8;
	}
	while (coef_index < _length);
}



template <class DT>
void	FFTReal <DT>::compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);
	assert (pass >= 3);
	assert (pass < _nbr_bits);

	if (pass <= TRIGO_BD_LIMIT)
	{
		compute_direct_pass_n_lut (df, sf, pass);
	}
	else
	{
		compute_direct_pass_n_osc (df, sf, pass);
	}
}



template <class DT>
void	FFTReal <DT>::compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);
	assert (pass >= 3);
	assert (pass < _nbr_bits);

	const long		nbr_coef = 1 << pass;
	const long		h_nbr_coef = nbr_coef >> 1;
	const long		d_nbr_coef = nbr_coef << 1;
	long				coef_index = 0;
	const DataType	* const	cos_ptr = get_trigo_ptr (pass);
	do
	{
		const DataType	* const	sf1r = sf + coef_index;
		const DataType	* const	sf2r = sf1r + nbr_coef;
		DataType			* const	dfr = df + coef_index;
		DataType			* const	dfi = dfr + nbr_coef;

		// Extreme coefficients are always real
		dfr [0] = sf1r [0] + sf2r [0];
		dfi [0] = sf1r [0] - sf2r [0];	// dfr [nbr_coef] =
		dfr [h_nbr_coef] = sf1r [h_nbr_coef];
		dfi [h_nbr_coef] = sf2r [h_nbr_coef];

		// Others are conjugate complex numbers
		const DataType * const	sf1i = sf1r + h_nbr_coef;
		const DataType * const	sf2i = sf1i + nbr_coef;
		for (long i = 1; i < h_nbr_coef; ++ i)
		{
			const DataType	c = cos_ptr [i];					// cos (i*PI/nbr_coef);
			const DataType	s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);
			DataType	 		v;

			v = sf2r [i] * c - sf2i [i] * s;
			dfr [i] = sf1r [i] + v;
			dfi [-i] = sf1r [i] - v;	// dfr [nbr_coef - i] =

			v = sf2r [i] * s + sf2i [i] * c;
			dfi [i] = v + sf1i [i];
			dfi [nbr_coef - i] = v - sf1i [i];
		}

		coef_index += d_nbr_coef;
	}
	while (coef_index < _length);
}



template <class DT>
void	FFTReal <DT>::compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);
	assert (pass > TRIGO_BD_LIMIT);
	assert (pass < _nbr_bits);

	const long		nbr_coef = 1 << pass;
	const long		h_nbr_coef = nbr_coef >> 1;
	const long		d_nbr_coef = nbr_coef << 1;
	long				coef_index = 0;
	OscType &		osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
	do
	{
		const DataType	* const	sf1r = sf + coef_index;
		const DataType	* const	sf2r = sf1r + nbr_coef;
		DataType			* const	dfr = df + coef_index;
		DataType			* const	dfi = dfr + nbr_coef;

		osc.clear_buffers ();

		// Extreme coefficients are always real
		dfr [0] = sf1r [0] + sf2r [0];
		dfi [0] = sf1r [0] - sf2r [0];	// dfr [nbr_coef] =
		dfr [h_nbr_coef] = sf1r [h_nbr_coef];
		dfi [h_nbr_coef] = sf2r [h_nbr_coef];

		// Others are conjugate complex numbers
		const DataType * const	sf1i = sf1r + h_nbr_coef;
		const DataType * const	sf2i = sf1i + nbr_coef;
		for (long i = 1; i < h_nbr_coef; ++ i)
		{
			osc.step ();
			const DataType	c = osc.get_cos ();
			const DataType	s = osc.get_sin ();
			DataType	 		v;

			v = sf2r [i] * c - sf2i [i] * s;
			dfr [i] = sf1r [i] + v;
			dfi [-i] = sf1r [i] - v;	// dfr [nbr_coef - i] =

			v = sf2r [i] * s + sf2i [i] * c;
			dfi [i] = v + sf1i [i];
			dfi [nbr_coef - i] = v - sf1i [i];
		}

		coef_index += d_nbr_coef;
	}
	while (coef_index < _length);
}



// Transform in several pass
template <class DT>
void	FFTReal <DT>::compute_ifft_general (const DataType f [], DataType x []) const
{
	assert (f != 0);
	assert (f != use_buffer ());
	assert (x != 0);
	assert (x != use_buffer ());
	assert (x != f);

	DataType *		sf = const_cast <DataType *> (f);
	DataType *		df;
	DataType *		df_temp;

	if (_nbr_bits & 1)
	{
		df = use_buffer ();
		df_temp = x;
	}
	else
	{
		df = x;
		df_temp = use_buffer ();
	}

	for (int pass = _nbr_bits - 1; pass >= 3; -- pass)
	{
		compute_inverse_pass_n (df, sf, pass);

		if (pass < _nbr_bits - 1)
		{
			DataType	* const	temp_ptr = df;
			df = sf;
			sf = temp_ptr;
		}
		else
		{
			sf = df;
			df = df_temp;
		}
	}

	compute_inverse_pass_3 (df, sf);
	compute_inverse_pass_1_2 (x, df);
}



template <class DT>
void	FFTReal <DT>::compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);
	assert (pass >= 3);
	assert (pass < _nbr_bits);

	if (pass <= TRIGO_BD_LIMIT)
	{
		compute_inverse_pass_n_lut (df, sf, pass);
	}
	else
	{
		compute_inverse_pass_n_osc (df, sf, pass);
	}
}



template <class DT>
void	FFTReal <DT>::compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);
	assert (pass >= 3);
	assert (pass < _nbr_bits);

	const long		nbr_coef = 1 << pass;
	const long		h_nbr_coef = nbr_coef >> 1;
	const long		d_nbr_coef = nbr_coef << 1;
	long				coef_index = 0;
	const DataType * const	cos_ptr = get_trigo_ptr (pass);
	do
	{
		const DataType	* const	sfr = sf + coef_index;
		const DataType	* const	sfi = sfr + nbr_coef;
		DataType			* const	df1r = df + coef_index;
		DataType			* const	df2r = df1r + nbr_coef;

		// Extreme coefficients are always real
		df1r [0] = sfr [0] + sfi [0];		// + sfr [nbr_coef]
		df2r [0] = sfr [0] - sfi [0];		// - sfr [nbr_coef]
		df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
		df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;

		// Others are conjugate complex numbers
		DataType * const	df1i = df1r + h_nbr_coef;
		DataType * const	df2i = df1i + nbr_coef;
		for (long i = 1; i < h_nbr_coef; ++ i)
		{
			df1r [i] = sfr [i] + sfi [-i];		// + sfr [nbr_coef - i]
			df1i [i] = sfi [i] - sfi [nbr_coef - i];

			const DataType	c = cos_ptr [i];					// cos (i*PI/nbr_coef);
			const DataType	s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);
			const DataType	vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]
			const DataType	vi = sfi [i] + sfi [nbr_coef - i];

			df2r [i] = vr * c + vi * s;
			df2i [i] = vi * c - vr * s;
		}

		coef_index += d_nbr_coef;
	}
	while (coef_index < _length);
}



template <class DT>
void	FFTReal <DT>::compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);
	assert (pass > TRIGO_BD_LIMIT);
	assert (pass < _nbr_bits);

	const long		nbr_coef = 1 << pass;
	const long		h_nbr_coef = nbr_coef >> 1;
	const long		d_nbr_coef = nbr_coef << 1;
	long				coef_index = 0;
	OscType &		osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
	do
	{
		const DataType	* const	sfr = sf + coef_index;
		const DataType	* const	sfi = sfr + nbr_coef;
		DataType			* const	df1r = df + coef_index;
		DataType			* const	df2r = df1r + nbr_coef;

		osc.clear_buffers ();

		// Extreme coefficients are always real
		df1r [0] = sfr [0] + sfi [0];		// + sfr [nbr_coef]
		df2r [0] = sfr [0] - sfi [0];		// - sfr [nbr_coef]
		df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
		df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;

		// Others are conjugate complex numbers
		DataType * const	df1i = df1r + h_nbr_coef;
		DataType * const	df2i = df1i + nbr_coef;
		for (long i = 1; i < h_nbr_coef; ++ i)
		{
			df1r [i] = sfr [i] + sfi [-i];		// + sfr [nbr_coef - i]
			df1i [i] = sfi [i] - sfi [nbr_coef - i];

			osc.step ();
			const DataType	c = osc.get_cos ();
			const DataType	s = osc.get_sin ();
			const DataType	vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]
			const DataType	vi = sfi [i] + sfi [nbr_coef - i];

			df2r [i] = vr * c + vi * s;
			df2i [i] = vi * c - vr * s;
		}

		coef_index += d_nbr_coef;
	}
	while (coef_index < _length);
}



template <class DT>
void	FFTReal <DT>::compute_inverse_pass_3 (DataType df [], const DataType sf []) const
{
	assert (df != 0);
	assert (sf != 0);
	assert (df != sf);

	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);
	long				coef_index = 0;
	do
	{
		df [coef_index] = sf [coef_index] + sf [coef_index + 4];
		df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
		df [coef_index + 2] = sf [coef_index + 2] * 2;
		df [coef_index + 6] = sf [coef_index + 6] * 2;

		df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];
		df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7];

		const DataType	vr = sf [coef_index + 1] - sf [coef_index + 3];
		const DataType	vi = sf [coef_index + 5] + sf [coef_index + 7];

		df [coef_index + 5] = (vr + vi) * sqrt2_2;
		df [coef_index + 7] = (vi - vr) * sqrt2_2;

		coef_index += 8;
	}
	while (coef_index < _length);
}



template <class DT>
void	FFTReal <DT>::compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const
{
	assert (x != 0);
	assert (sf != 0);
	assert (x != sf);

	const long *	bit_rev_lut_ptr = get_br_ptr ();
	const DataType *	sf2 = sf;
	long				coef_index = 0;
	do
	{
		{
			const DataType	b_0 = sf2 [0] + sf2 [2];
			const DataType	b_2 = sf2 [0] - sf2 [2];
			const DataType	b_1 = sf2 [1] * 2;
			const DataType	b_3 = sf2 [3] * 2;

			x [bit_rev_lut_ptr [0]] = b_0 + b_1;
			x [bit_rev_lut_ptr [1]] = b_0 - b_1;
			x [bit_rev_lut_ptr [2]] = b_2 + b_3;
			x [bit_rev_lut_ptr [3]] = b_2 - b_3;
		}
		{
			const DataType	b_0 = sf2 [4] + sf2 [6];
			const DataType	b_2 = sf2 [4] - sf2 [6];
			const DataType	b_1 = sf2 [5] * 2;
			const DataType	b_3 = sf2 [7] * 2;

			x [bit_rev_lut_ptr [4]] = b_0 + b_1;
			x [bit_rev_lut_ptr [5]] = b_0 - b_1;
			x [bit_rev_lut_ptr [6]] = b_2 + b_3;
			x [bit_rev_lut_ptr [7]] = b_2 - b_3;
		}

		sf2 += 8;
		coef_index += 8;
		bit_rev_lut_ptr += 8;
	}
	while (coef_index < _length);
}



#endif	// FFTReal_CODEHEADER_INCLUDED

#undef FFTReal_CURRENT_CODEHEADER



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