BogaudioModules

FFTReal.hpp (21034B)

---

 /*****************************************************************************
 
         FFTReal.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_FFTReal_CURRENT_CODEHEADER)
 	#error Recursive inclusion of FFTReal code header.
 #endif
 #define	ffft_FFTReal_CURRENT_CODEHEADER
 
 #if ! defined (ffft_FFTReal_CODEHEADER_INCLUDED)
 #define	ffft_FFTReal_CODEHEADER_INCLUDED
 
 
 
 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
 
 #include	<cassert>
 #include	<cmath>
 
 
 
 namespace ffft
 {
 
 
 
 static inline bool	FFTReal_is_pow2 (long x)
 {
 	assert (x > 0);
 
 	return  ((x & -x) == x);
 }
 
 
 
 static inline int	FFTReal_get_next_pow2 (long x)
 {
 	--x;
 
 	int				p = 0;
 	while ((x & ~0xFFFFL) != 0)
 	{
 		p += 16;
 		x >>= 16;
 	}
 	while ((x & ~0xFL) != 0)
 	{
 		p += 4;
 		x >>= 4;
 	}
 	while (x > 0)
 	{
 		++p;
 		x >>= 1;
 	}
 
 	return (p);
 }
 
 
 
 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
 
 
 
 /*
 ==============================================================================
 Name: ctor
 Input parameters:
 	- length: length of the array on which we want to do a FFT. Range: power of
 		2 only, > 0.
 Throws: std::bad_alloc
 ==============================================================================
 */
 
 template <class DT>
 FFTReal <DT>::FFTReal (long length)
 :	_length (length)
 ,	_nbr_bits (FFTReal_get_next_pow2 (length))
 ,	_br_lut ()
 ,	_trigo_lut ()
 ,	_buffer (length)
 ,	_trigo_osc ()
 {
 	assert (FFTReal_is_pow2 (length));
 	assert (_nbr_bits <= MAX_BIT_DEPTH);
 
 	init_br_lut ();
 	init_trigo_lut ();
 	init_trigo_osc ();
 }
 
 
 
 /*
 ==============================================================================
 Name: get_length
 Description:
 	Returns the number of points processed by this FFT object.
 Returns: The number of points, power of 2, > 0.
 Throws: Nothing
 ==============================================================================
 */
 
 template <class DT>
 long	FFTReal <DT>::get_length () const
 {
 	return (_length);
 }
 
 
 
 /*
 ==============================================================================
 Name: do_fft
 Description:
 	Compute the FFT of the array.
 Input parameters:
 	- x: pointer on the source array (time).
 Output parameters:
 	- f: pointer on the destination array (frequencies).
 		f [0...length(x)/2] = real values,
 		f [length(x)/2+1...length(x)-1] = negative imaginary values of
 		coefficents 1...length(x)/2-1.
 Throws: Nothing
 ==============================================================================
 */
 
 template <class DT>
 void	FFTReal <DT>::do_fft (DataType f [], const DataType x []) const
 {
 	assert (f != 0);
 	assert (f != use_buffer ());
 	assert (x != 0);
 	assert (x != use_buffer ());
 	assert (x != f);
 
 	// General case
 	if (_nbr_bits > 2)
 	{
 		compute_fft_general (f, x);
 	}
 
 	// 4-point FFT
 	else if (_nbr_bits == 2)
 	{
 		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
 	else if (_nbr_bits == 1)
 	{
 		f [0] = x [0] + x [1];
 		f [1] = x [0] - x [1];
 	}
 
 	// 1-point FFT
 	else
 	{
 		f [0] = x [0];
 	}
 }
 
 
 
 /*
 ==============================================================================
 Name: do_ifft
 Description:
 	Compute the inverse FFT of the array. Note that data must be post-scaled:
 	IFFT (FFT (x)) = x * length (x).
 Input parameters:
 	- f: pointer on the source array (frequencies).
 		f [0...length(x)/2] = real values
 		f [length(x)/2+1...length(x)-1] = negative imaginary values of
 		coefficents 1...length(x)/2-1.
 Output parameters:
 	- x: pointer on the destination array (time).
 Throws: Nothing
 ==============================================================================
 */
 
 template <class DT>
 void	FFTReal <DT>::do_ifft (const DataType f [], DataType x []) const
 {
 	assert (f != 0);
 	assert (f != use_buffer ());
 	assert (x != 0);
 	assert (x != use_buffer ());
 	assert (x != f);
 
 	// General case
 	if (_nbr_bits > 2)
 	{
 		compute_ifft_general (f, x);
 	}
 
 	// 4-point IFFT
 	else if (_nbr_bits == 2)
 	{
 		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
 	else if (_nbr_bits == 1)
 	{
 		x [0] = f [0] + f [1];
 		x [1] = f [0] - f [1];
 	}
 
 	// 1-point IFFT
 	else
 	{
 		x [0] = f [0];
 	}
 }
 
 
 
 /*
 ==============================================================================
 Name: rescale
 Description:
 	Scale an array by divide each element by its length. This function should
 	be called after FFT + IFFT.
 Input parameters:
 	- x: pointer on array to rescale (time or frequency).
 Throws: Nothing
 ==============================================================================
 */
 
 template <class DT>
 void	FFTReal <DT>::rescale (DataType x []) const
 {
 	const DataType	mul = DataType (1.0 / _length);
 
 	if (_length < 4)
 	{
 		long				i = _length - 1;
 		do
 		{
 			x [i] *= mul;
 			--i;
 		}
 		while (i >= 0);
 	}
 
 	else
 	{
 		assert ((_length & 3) == 0);
 
 		// Could be optimized with SIMD instruction sets (needs alignment check)
 		long				i = _length - 4;
 		do
 		{
 			x [i + 0] *= mul;
 			x [i + 1] *= mul;
 			x [i + 2] *= mul;
 			x [i + 3] *= mul;
 			i -= 4;
 		}
 		while (i >= 0);
 	}
 }
 
 
 
 /*
 ==============================================================================
 Name: use_buffer
 Description:
 	Access the internal buffer, whose length is the FFT one.
 	Buffer content will be erased at each do_fft() / do_ifft() call!
 	This buffer cannot be used as:
 		- source for FFT or IFFT done with this object
 		- destination for FFT or IFFT done with this object
 Returns:
 	Buffer start address
 Throws: Nothing
 ==============================================================================
 */
 
 template <class DT>
 typename FFTReal <DT>::DataType *	FFTReal <DT>::use_buffer () const
 {
 	return (&_buffer [0]);
 }
 
 
 
 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
 
 
 
 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
 
 
 
 template <class DT>
 void	FFTReal <DT>::init_br_lut ()
 {
 	const long		length = 1L << _nbr_bits;
 	_br_lut.resize (length);
 
 	_br_lut [0] = 0;
 	long				br_index = 0;
 	for (long cnt = 1; cnt < length; ++cnt)
 	{
 		// ++br_index (bit reversed)
 		long				bit = length >> 1;
 		while (((br_index ^= bit) & bit) == 0)
 		{
 			bit >>= 1;
 		}
 
 		_br_lut [cnt] = br_index;
 	}
 }
 
 
 
 template <class DT>
 void	FFTReal <DT>::init_trigo_lut ()
 {
 	using namespace std;
 
 	if (_nbr_bits > 3)
 	{
 		const long		total_len = (1L << (_nbr_bits - 1)) - 4;
 		_trigo_lut.resize (total_len);
 
 		for (int level = 3; level < _nbr_bits; ++level)
 		{
 			const long		level_len = 1L << (level - 1);
 			DataType	* const	level_ptr =
 				&_trigo_lut [get_trigo_level_index (level)];
 			const double	mul = PI / (level_len << 1);
 
 			for (long i = 0; i < level_len; ++ i)
 			{
 				level_ptr [i] = static_cast <DataType> (cos (i * mul));
 			}
 		}
 	}
 }
 
 
 
 template <class DT>
 void	FFTReal <DT>::init_trigo_osc ()
 {
 	const int		nbr_osc = _nbr_bits - TRIGO_BD_LIMIT;
 	if (nbr_osc > 0)
 	{
 		_trigo_osc.resize (nbr_osc);
 
 		for (int osc_cnt = 0; osc_cnt < nbr_osc; ++osc_cnt)
 		{
 			OscType &		osc = _trigo_osc [osc_cnt];
 
 			const long		len = 1L << (TRIGO_BD_LIMIT + osc_cnt);
 			const double	mul = (0.5 * PI) / len;
 			osc.set_step (mul);
 		}
 	}
 }
 
 
 
 template <class DT>
 const long *	FFTReal <DT>::get_br_ptr () const
 {
 	return (&_br_lut [0]);
 }
 
 
 
 template <class DT>
 const typename FFTReal <DT>::DataType *	FFTReal <DT>::get_trigo_ptr (int level) const
 {
 	assert (level >= 3);
 
 	return (&_trigo_lut [get_trigo_level_index (level)]);
 }
 
 
 
 template <class DT>
 long	FFTReal <DT>::get_trigo_level_index (int level) const
 {
 	assert (level >= 3);
 
 	return ((1L << (level - 1)) - 4);
 }
 
 
 
 // Transform in several passes
 template <class DT>
 void	FFTReal <DT>::compute_fft_general (DataType f [], const DataType x []) const
 {
 	assert (f != 0);
 	assert (f != use_buffer ());
 	assert (x != 0);
 	assert (x != use_buffer ());
 	assert (x != f);
 
 	DataType *		sf;
 	DataType *		df;
 
 	if ((_nbr_bits & 1) != 0)
 	{
 		df = use_buffer ();
 		sf = f;
 	}
 	else
 	{
 		df = f;
 		sf = use_buffer ();
 	}
 
 	compute_direct_pass_1_2 (df, x);
 	compute_direct_pass_3 (sf, df);
 
 	for (int pass = 3; pass < _nbr_bits; ++ pass)
 	{
 		compute_direct_pass_n (df, sf, pass);
 
 		DataType * const	temp_ptr = df;
 		df = sf;
 		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);
 }
 
 
 
 }	// namespace ffft
 
 
 
 #endif	// ffft_FFTReal_CODEHEADER_INCLUDED
 
 #undef ffft_FFTReal_CURRENT_CODEHEADER
 
 
 
 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/