lbajardsilogic@79: /*****************************************************************************
lbajardsilogic@79: *                                                                            *
lbajardsilogic@79: *       DIGITAL SIGNAL PROCESSING TOOLS                                      *
lbajardsilogic@79: *       Version 1.01, 1999/11/07                                             *
lbajardsilogic@79: *       (c) 1999 - Laurent de Soras                                          *
lbajardsilogic@79: *                                                                            *
lbajardsilogic@79: *       FFTReal.cpp                                                          *
lbajardsilogic@79: *       Fourier transformation of real number arrays.                        *
lbajardsilogic@79: *       Portable ISO C++                                                     *
lbajardsilogic@79: *                                                                            *
lbajardsilogic@79: * Tab = 3                                                                    *
lbajardsilogic@79: *****************************************************************************/
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
lbajardsilogic@79: 
lbajardsilogic@79: #include	"FFTReal.h"
lbajardsilogic@79: 
lbajardsilogic@79: #include	<cassert>
lbajardsilogic@79: #include	<cmath>
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: #if defined (_MSC_VER)
lbajardsilogic@79: #pragma pack (push, 8)
lbajardsilogic@79: #endif	// _MSC_VER
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*\\\ PUBLIC MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: Constructor                                                   */
lbajardsilogic@79: /*      Input parameters:                                                   */
lbajardsilogic@79: /*        - length: length of the array on which we want to do a FFT.       */
lbajardsilogic@79: /*                  Range: power of 2 only, > 0.                            */
lbajardsilogic@79: /*      Throws: std::bad_alloc, anything                                    */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: FFTReal::FFTReal (const long length)
lbajardsilogic@79: :	_length (length)
lbajardsilogic@79: ,	_nbr_bits (int (floor (log ((float)length) / log (2.0f) + 0.5f)))
lbajardsilogic@79: ,	_bit_rev_lut (int (floor (log ((float)length) / log (2.0f) + 0.5f)))
lbajardsilogic@79: ,	_trigo_lut (int (floor (log ((float)length) / log (2.0f) + 0.5f)))
lbajardsilogic@79: ,	_sqrt2_2 (flt_t (sqrt (2.0f) * 0.5))
lbajardsilogic@79: {
lbajardsilogic@79: 	assert ((1L << _nbr_bits) == length);
lbajardsilogic@79: 
lbajardsilogic@79: 	_buffer_ptr = 0;
lbajardsilogic@79: 	if (_nbr_bits > 2)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		_buffer_ptr = new flt_t [_length];
lbajardsilogic@79: 	}
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: Destructor                                                    */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: FFTReal::~FFTReal (void)
lbajardsilogic@79: {
lbajardsilogic@79: 	delete [] _buffer_ptr;
lbajardsilogic@79: 	_buffer_ptr = 0;
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: do_fft                                                        */
lbajardsilogic@79: /*      Description: Compute the FFT of the array.                          */
lbajardsilogic@79: /*      Input parameters:                                                   */
lbajardsilogic@79: /*        - x: pointer on the source array (time).                          */
lbajardsilogic@79: /*      Output parameters:                                                  */
lbajardsilogic@79: /*        - f: pointer on the destination array (frequencies).              */
lbajardsilogic@79: /*             f [0...length(x)/2] = real values,                           */
lbajardsilogic@79: /*             f [length(x)/2+1...length(x)-1] = imaginary values of        */
lbajardsilogic@79: /*               coefficents 1...length(x)/2-1.                             */
lbajardsilogic@79: /*      Throws: Nothing                                                     */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: void	FFTReal::do_fft (flt_t f [], const flt_t x []) const
lbajardsilogic@79: {
lbajardsilogic@79: 
lbajardsilogic@79: /*______________________________________________
lbajardsilogic@79:  *
lbajardsilogic@79:  * General case
lbajardsilogic@79:  *______________________________________________
lbajardsilogic@79:  */
lbajardsilogic@79: 
lbajardsilogic@79: 	if (_nbr_bits > 2)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		flt_t *			sf;
lbajardsilogic@79: 		flt_t *			df;
lbajardsilogic@79: 
lbajardsilogic@79: 		if (_nbr_bits & 1)
lbajardsilogic@79: 		{
lbajardsilogic@79: 			df = _buffer_ptr;
lbajardsilogic@79: 			sf = f;
lbajardsilogic@79: 		}
lbajardsilogic@79: 		else
lbajardsilogic@79: 		{
lbajardsilogic@79: 			df = f;
lbajardsilogic@79: 			sf = _buffer_ptr;
lbajardsilogic@79: 		}
lbajardsilogic@79: 
lbajardsilogic@79: 		/* Do the transformation in several pass */
lbajardsilogic@79: 		{
lbajardsilogic@79: 			int		pass;
lbajardsilogic@79: 			long		nbr_coef;
lbajardsilogic@79: 			long		h_nbr_coef;
lbajardsilogic@79: 			long		d_nbr_coef;
lbajardsilogic@79: 			long		coef_index;
lbajardsilogic@79: 
lbajardsilogic@79: 			/* First and second pass at once */
lbajardsilogic@79: 			{
lbajardsilogic@79: 				const long * const	bit_rev_lut_ptr = _bit_rev_lut.get_ptr ();
lbajardsilogic@79: 				coef_index = 0;
lbajardsilogic@79: 				do
lbajardsilogic@79: 				{
lbajardsilogic@79: 					const long		rev_index_0 = bit_rev_lut_ptr [coef_index];
lbajardsilogic@79: 					const long		rev_index_1 = bit_rev_lut_ptr [coef_index + 1];
lbajardsilogic@79: 					const long		rev_index_2 = bit_rev_lut_ptr [coef_index + 2];
lbajardsilogic@79: 					const long		rev_index_3 = bit_rev_lut_ptr [coef_index + 3];
lbajardsilogic@79: 
lbajardsilogic@79: 					flt_t	* const	df2 = df + coef_index;
lbajardsilogic@79: 					df2 [1] = x [rev_index_0] - x [rev_index_1];
lbajardsilogic@79: 					df2 [3] = x [rev_index_2] - x [rev_index_3];
lbajardsilogic@79: 
lbajardsilogic@79: 					const flt_t		sf_0 = x [rev_index_0] + x [rev_index_1];
lbajardsilogic@79: 					const flt_t		sf_2 = x [rev_index_2] + x [rev_index_3];
lbajardsilogic@79: 
lbajardsilogic@79: 					df2 [0] = sf_0 + sf_2;
lbajardsilogic@79: 					df2 [2] = sf_0 - sf_2;
lbajardsilogic@79: 					
lbajardsilogic@79: 					coef_index += 4;
lbajardsilogic@79: 				}
lbajardsilogic@79: 				while (coef_index < _length);
lbajardsilogic@79: 			}
lbajardsilogic@79: 
lbajardsilogic@79: 			/* Third pass */
lbajardsilogic@79: 			{
lbajardsilogic@79: 				coef_index = 0;
lbajardsilogic@79: 				const flt_t		sqrt2_2 = _sqrt2_2;
lbajardsilogic@79: 				do
lbajardsilogic@79: 				{
lbajardsilogic@79: 					flt_t				v;
lbajardsilogic@79: 
lbajardsilogic@79: 					sf [coef_index] = df [coef_index] + df [coef_index + 4];
lbajardsilogic@79: 					sf [coef_index + 4] = df [coef_index] - df [coef_index + 4];
lbajardsilogic@79: 					sf [coef_index + 2] = df [coef_index + 2];
lbajardsilogic@79: 					sf [coef_index + 6] = df [coef_index + 6];
lbajardsilogic@79: 
lbajardsilogic@79: 					v = (df [coef_index + 5] - df [coef_index + 7]) * sqrt2_2;
lbajardsilogic@79: 					sf [coef_index + 1] = df [coef_index + 1] + v;
lbajardsilogic@79: 					sf [coef_index + 3] = df [coef_index + 1] - v;
lbajardsilogic@79: 
lbajardsilogic@79: 					v = (df [coef_index + 5] + df [coef_index + 7]) * sqrt2_2;
lbajardsilogic@79: 					sf [coef_index + 5] = v + df [coef_index + 3];
lbajardsilogic@79: 					sf [coef_index + 7] = v - df [coef_index + 3];
lbajardsilogic@79: 
lbajardsilogic@79: 					coef_index += 8;
lbajardsilogic@79: 				}
lbajardsilogic@79: 				while (coef_index < _length);
lbajardsilogic@79: 			}
lbajardsilogic@79: 
lbajardsilogic@79: 			/* Next pass */
lbajardsilogic@79: 			for (pass = 3; pass < _nbr_bits; ++pass)
lbajardsilogic@79: 			{
lbajardsilogic@79: 				coef_index = 0;
lbajardsilogic@79: 				nbr_coef = 1 << pass;
lbajardsilogic@79: 				h_nbr_coef = nbr_coef >> 1;
lbajardsilogic@79: 				d_nbr_coef = nbr_coef << 1;
lbajardsilogic@79: 				const flt_t	* const	cos_ptr = _trigo_lut.get_ptr (pass);
lbajardsilogic@79: 				do
lbajardsilogic@79: 				{
lbajardsilogic@79: 					long				i;
lbajardsilogic@79: 					const flt_t	*	const sf1r = sf + coef_index;
lbajardsilogic@79: 					const flt_t	*	const sf2r = sf1r + nbr_coef;
lbajardsilogic@79: 					flt_t *			const dfr = df + coef_index;
lbajardsilogic@79: 					flt_t *			const dfi = dfr + nbr_coef;
lbajardsilogic@79: 
lbajardsilogic@79: 					/* Extreme coefficients are always real */
lbajardsilogic@79: 					dfr [0] = sf1r [0] + sf2r [0];
lbajardsilogic@79: 					dfi [0] = sf1r [0] - sf2r [0];	// dfr [nbr_coef] =
lbajardsilogic@79: 					dfr [h_nbr_coef] = sf1r [h_nbr_coef];
lbajardsilogic@79: 					dfi [h_nbr_coef] = sf2r [h_nbr_coef];
lbajardsilogic@79: 
lbajardsilogic@79: 					/* Others are conjugate complex numbers */
lbajardsilogic@79: 					const flt_t	* const	sf1i = sf1r + h_nbr_coef;
lbajardsilogic@79: 					const flt_t	* const	sf2i = sf1i + nbr_coef;
lbajardsilogic@79: 					for (i = 1; i < h_nbr_coef; ++ i)
lbajardsilogic@79: 					{
lbajardsilogic@79: 						const flt_t		c = cos_ptr [i];					// cos (i*PI/nbr_coef);
lbajardsilogic@79: 						const flt_t		s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);
lbajardsilogic@79: 						flt_t				v;
lbajardsilogic@79: 
lbajardsilogic@79: 						v = sf2r [i] * c - sf2i [i] * s;
lbajardsilogic@79: 						dfr [i] = sf1r [i] + v;
lbajardsilogic@79: 						dfi [-i] = sf1r [i] - v;	// dfr [nbr_coef - i] =
lbajardsilogic@79: 
lbajardsilogic@79: 						v = sf2r [i] * s + sf2i [i] * c;
lbajardsilogic@79: 						dfi [i] = v + sf1i [i];
lbajardsilogic@79: 						dfi [nbr_coef - i] = v - sf1i [i];
lbajardsilogic@79: 					}
lbajardsilogic@79: 
lbajardsilogic@79: 					coef_index += d_nbr_coef;
lbajardsilogic@79: 				}
lbajardsilogic@79: 				while (coef_index < _length);
lbajardsilogic@79: 
lbajardsilogic@79: 				/* Prepare to the next pass */
lbajardsilogic@79: 				{
lbajardsilogic@79: 					flt_t	* const		temp_ptr = df;
lbajardsilogic@79: 					df = sf;
lbajardsilogic@79: 					sf = temp_ptr;
lbajardsilogic@79: 				}
lbajardsilogic@79: 			}
lbajardsilogic@79: 		}
lbajardsilogic@79: 	}
lbajardsilogic@79: 
lbajardsilogic@79: /*______________________________________________
lbajardsilogic@79:  *
lbajardsilogic@79:  * Special cases
lbajardsilogic@79:  *______________________________________________
lbajardsilogic@79:  */
lbajardsilogic@79: 
lbajardsilogic@79: 	/* 4-point FFT */
lbajardsilogic@79: 	else if (_nbr_bits == 2)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		f [1] = x [0] - x [2];
lbajardsilogic@79: 		f [3] = x [1] - x [3];
lbajardsilogic@79: 
lbajardsilogic@79: 		const flt_t			b_0 = x [0] + x [2];
lbajardsilogic@79: 		const flt_t			b_2 = x [1] + x [3];
lbajardsilogic@79: 		
lbajardsilogic@79: 		f [0] = b_0 + b_2;
lbajardsilogic@79: 		f [2] = b_0 - b_2;
lbajardsilogic@79: 	}
lbajardsilogic@79: 
lbajardsilogic@79: 	/* 2-point FFT */
lbajardsilogic@79: 	else if (_nbr_bits == 1)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		f [0] = x [0] + x [1];
lbajardsilogic@79: 		f [1] = x [0] - x [1];
lbajardsilogic@79: 	}
lbajardsilogic@79: 
lbajardsilogic@79: 	/* 1-point FFT */
lbajardsilogic@79: 	else
lbajardsilogic@79: 	{
lbajardsilogic@79: 		f [0] = x [0];
lbajardsilogic@79: 	}
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: do_ifft                                                       */
lbajardsilogic@79: /*      Description: Compute the inverse FFT of the array. Notice that      */
lbajardsilogic@79: /*                   IFFT (FFT (x)) = x * length (x). Data must be          */
lbajardsilogic@79: /*                   post-scaled.                                           */
lbajardsilogic@79: /*      Input parameters:                                                   */
lbajardsilogic@79: /*        - f: pointer on the source array (frequencies).                   */
lbajardsilogic@79: /*             f [0...length(x)/2] = real values,                           */
lbajardsilogic@79: /*             f [length(x)/2+1...length(x)] = imaginary values of          */
lbajardsilogic@79: /*               coefficents 1...length(x)-1.                               */
lbajardsilogic@79: /*      Output parameters:                                                  */
lbajardsilogic@79: /*        - x: pointer on the destination array (time).                     */
lbajardsilogic@79: /*      Throws: Nothing                                                     */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: void	FFTReal::do_ifft (const flt_t f [], flt_t x []) const
lbajardsilogic@79: {
lbajardsilogic@79: 
lbajardsilogic@79: /*______________________________________________
lbajardsilogic@79:  *
lbajardsilogic@79:  * General case
lbajardsilogic@79:  *______________________________________________
lbajardsilogic@79:  */
lbajardsilogic@79: 
lbajardsilogic@79: 	if (_nbr_bits > 2)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		flt_t *			sf = const_cast <flt_t *> (f);
lbajardsilogic@79: 		flt_t *			df;
lbajardsilogic@79: 		flt_t *			df_temp;
lbajardsilogic@79: 
lbajardsilogic@79: 		if (_nbr_bits & 1)
lbajardsilogic@79: 		{
lbajardsilogic@79: 			df = _buffer_ptr;
lbajardsilogic@79: 			df_temp = x;
lbajardsilogic@79: 		}
lbajardsilogic@79: 		else
lbajardsilogic@79: 		{
lbajardsilogic@79: 			df = x;
lbajardsilogic@79: 			df_temp = _buffer_ptr;
lbajardsilogic@79: 		}
lbajardsilogic@79: 
lbajardsilogic@79: 		/* Do the transformation in several pass */
lbajardsilogic@79: 		{
lbajardsilogic@79: 			int			pass;
lbajardsilogic@79: 			long			nbr_coef;
lbajardsilogic@79: 			long			h_nbr_coef;
lbajardsilogic@79: 			long			d_nbr_coef;
lbajardsilogic@79: 			long			coef_index;
lbajardsilogic@79: 
lbajardsilogic@79: 			/* First pass */
lbajardsilogic@79: 			for (pass = _nbr_bits - 1; pass >= 3; --pass)
lbajardsilogic@79: 			{
lbajardsilogic@79: 				coef_index = 0;
lbajardsilogic@79: 				nbr_coef = 1 << pass;
lbajardsilogic@79: 				h_nbr_coef = nbr_coef >> 1;
lbajardsilogic@79: 				d_nbr_coef = nbr_coef << 1;
lbajardsilogic@79: 				const flt_t	*const cos_ptr = _trigo_lut.get_ptr (pass);
lbajardsilogic@79: 				do
lbajardsilogic@79: 				{
lbajardsilogic@79: 					long				i;
lbajardsilogic@79: 					const flt_t	*	const sfr = sf + coef_index;
lbajardsilogic@79: 					const flt_t	*	const sfi = sfr + nbr_coef;
lbajardsilogic@79: 					flt_t *			const df1r = df + coef_index;
lbajardsilogic@79: 					flt_t *			const df2r = df1r + nbr_coef;
lbajardsilogic@79: 
lbajardsilogic@79: 					/* Extreme coefficients are always real */
lbajardsilogic@79: 					df1r [0] = sfr [0] + sfi [0];		// + sfr [nbr_coef]
lbajardsilogic@79: 					df2r [0] = sfr [0] - sfi [0];		// - sfr [nbr_coef]
lbajardsilogic@79: 					df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
lbajardsilogic@79: 					df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
lbajardsilogic@79: 
lbajardsilogic@79: 					/* Others are conjugate complex numbers */
lbajardsilogic@79: 					flt_t * const	df1i = df1r + h_nbr_coef;
lbajardsilogic@79: 					flt_t * const	df2i = df1i + nbr_coef;
lbajardsilogic@79: 					for (i = 1; i < h_nbr_coef; ++ i)
lbajardsilogic@79: 					{
lbajardsilogic@79: 						df1r [i] = sfr [i] + sfi [-i];		// + sfr [nbr_coef - i]
lbajardsilogic@79: 						df1i [i] = sfi [i] - sfi [nbr_coef - i];
lbajardsilogic@79: 
lbajardsilogic@79: 						const flt_t		c = cos_ptr [i];					// cos (i*PI/nbr_coef);
lbajardsilogic@79: 						const flt_t		s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);
lbajardsilogic@79: 						const flt_t		vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]
lbajardsilogic@79: 						const flt_t		vi = sfi [i] + sfi [nbr_coef - i];
lbajardsilogic@79: 
lbajardsilogic@79: 						df2r [i] = vr * c + vi * s;
lbajardsilogic@79: 						df2i [i] = vi * c - vr * s;
lbajardsilogic@79: 					}
lbajardsilogic@79: 
lbajardsilogic@79: 					coef_index += d_nbr_coef;
lbajardsilogic@79: 				}
lbajardsilogic@79: 				while (coef_index < _length);
lbajardsilogic@79: 
lbajardsilogic@79: 				/* Prepare to the next pass */
lbajardsilogic@79: 				if (pass < _nbr_bits - 1)
lbajardsilogic@79: 				{
lbajardsilogic@79: 					flt_t	* const	temp_ptr = df;
lbajardsilogic@79: 					df = sf;
lbajardsilogic@79: 					sf = temp_ptr;
lbajardsilogic@79: 				}
lbajardsilogic@79: 				else
lbajardsilogic@79: 				{
lbajardsilogic@79: 					sf = df;
lbajardsilogic@79: 					df = df_temp;
lbajardsilogic@79: 				}
lbajardsilogic@79: 			}
lbajardsilogic@79: 
lbajardsilogic@79: 			/* Antepenultimate pass */
lbajardsilogic@79: 			{
lbajardsilogic@79: 				const flt_t		sqrt2_2 = _sqrt2_2;
lbajardsilogic@79: 				coef_index = 0;
lbajardsilogic@79: 				do
lbajardsilogic@79: 				{
lbajardsilogic@79: 					df [coef_index] = sf [coef_index] + sf [coef_index + 4];
lbajardsilogic@79: 					df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
lbajardsilogic@79: 					df [coef_index + 2] = sf [coef_index + 2] * 2;
lbajardsilogic@79: 					df [coef_index + 6] = sf [coef_index + 6] * 2;
lbajardsilogic@79: 
lbajardsilogic@79: 					df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];
lbajardsilogic@79: 					df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7];
lbajardsilogic@79: 
lbajardsilogic@79: 					const flt_t		vr = sf [coef_index + 1] - sf [coef_index + 3];
lbajardsilogic@79: 					const flt_t		vi = sf [coef_index + 5] + sf [coef_index + 7];
lbajardsilogic@79: 
lbajardsilogic@79: 					df [coef_index + 5] = (vr + vi) * sqrt2_2;
lbajardsilogic@79: 					df [coef_index + 7] = (vi - vr) * sqrt2_2;
lbajardsilogic@79: 
lbajardsilogic@79: 					coef_index += 8;
lbajardsilogic@79: 				}
lbajardsilogic@79: 				while (coef_index < _length);
lbajardsilogic@79: 			}
lbajardsilogic@79: 
lbajardsilogic@79: 			/* Penultimate and last pass at once */
lbajardsilogic@79: 			{
lbajardsilogic@79: 				coef_index = 0;
lbajardsilogic@79: 				const long *	bit_rev_lut_ptr = _bit_rev_lut.get_ptr ();
lbajardsilogic@79: 				const flt_t	*	sf2 = df;
lbajardsilogic@79: 				do
lbajardsilogic@79: 				{
lbajardsilogic@79: 					{
lbajardsilogic@79: 						const flt_t		b_0 = sf2 [0] + sf2 [2];
lbajardsilogic@79: 						const flt_t		b_2 = sf2 [0] - sf2 [2];
lbajardsilogic@79: 						const flt_t		b_1 = sf2 [1] * 2;
lbajardsilogic@79: 						const flt_t		b_3 = sf2 [3] * 2;
lbajardsilogic@79: 
lbajardsilogic@79: 						x [bit_rev_lut_ptr [0]] = b_0 + b_1;
lbajardsilogic@79: 						x [bit_rev_lut_ptr [1]] = b_0 - b_1;
lbajardsilogic@79: 						x [bit_rev_lut_ptr [2]] = b_2 + b_3;
lbajardsilogic@79: 						x [bit_rev_lut_ptr [3]] = b_2 - b_3;
lbajardsilogic@79: 					}
lbajardsilogic@79: 					{
lbajardsilogic@79: 						const flt_t		b_0 = sf2 [4] + sf2 [6];
lbajardsilogic@79: 						const flt_t		b_2 = sf2 [4] - sf2 [6];
lbajardsilogic@79: 						const flt_t		b_1 = sf2 [5] * 2;
lbajardsilogic@79: 						const flt_t		b_3 = sf2 [7] * 2;
lbajardsilogic@79: 
lbajardsilogic@79: 						x [bit_rev_lut_ptr [4]] = b_0 + b_1;
lbajardsilogic@79: 						x [bit_rev_lut_ptr [5]] = b_0 - b_1;
lbajardsilogic@79: 						x [bit_rev_lut_ptr [6]] = b_2 + b_3;
lbajardsilogic@79: 						x [bit_rev_lut_ptr [7]] = b_2 - b_3;
lbajardsilogic@79: 					}
lbajardsilogic@79: 
lbajardsilogic@79: 					sf2 += 8;
lbajardsilogic@79: 					coef_index += 8;
lbajardsilogic@79: 					bit_rev_lut_ptr += 8;
lbajardsilogic@79: 				}
lbajardsilogic@79: 				while (coef_index < _length);
lbajardsilogic@79: 			}
lbajardsilogic@79: 		}
lbajardsilogic@79: 	}
lbajardsilogic@79: 
lbajardsilogic@79: /*______________________________________________
lbajardsilogic@79:  *
lbajardsilogic@79:  * Special cases
lbajardsilogic@79:  *______________________________________________
lbajardsilogic@79:  */
lbajardsilogic@79: 
lbajardsilogic@79: 	/* 4-point IFFT */
lbajardsilogic@79: 	else if (_nbr_bits == 2)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		const flt_t		b_0 = f [0] + f [2];
lbajardsilogic@79: 		const flt_t		b_2 = f [0] - f [2];
lbajardsilogic@79: 
lbajardsilogic@79: 		x [0] = b_0 + f [1] * 2;
lbajardsilogic@79: 		x [2] = b_0 - f [1] * 2;
lbajardsilogic@79: 		x [1] = b_2 + f [3] * 2;
lbajardsilogic@79: 		x [3] = b_2 - f [3] * 2;
lbajardsilogic@79: 	}
lbajardsilogic@79: 
lbajardsilogic@79: 	/* 2-point IFFT */
lbajardsilogic@79: 	else if (_nbr_bits == 1)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		x [0] = f [0] + f [1];
lbajardsilogic@79: 		x [1] = f [0] - f [1];
lbajardsilogic@79: 	}
lbajardsilogic@79: 
lbajardsilogic@79: 	/* 1-point IFFT */
lbajardsilogic@79: 	else
lbajardsilogic@79: 	{
lbajardsilogic@79: 		x [0] = f [0];
lbajardsilogic@79: 	}
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: rescale                                                       */
lbajardsilogic@79: /*      Description: Scale an array by divide each element by its length.   */
lbajardsilogic@79: /*                   This function should be called after FFT + IFFT.       */
lbajardsilogic@79: /*      Input/Output parameters:                                            */
lbajardsilogic@79: /*        - x: pointer on array to rescale (time or frequency).             */
lbajardsilogic@79: /*      Throws: Nothing                                                     */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: void	FFTReal::rescale (flt_t x []) const
lbajardsilogic@79: {
lbajardsilogic@79: 	const double dpi=3.14179;
lbajardsilogic@79: 	const flt_t		mul = flt_t (1.0 / _length);
lbajardsilogic@79: 	long				i = _length - 1;
lbajardsilogic@79: 
lbajardsilogic@79: 	do
lbajardsilogic@79: 	{
lbajardsilogic@79: 		x [i] *= mul;
lbajardsilogic@79: 		//x [i] *= mul*(0.5*(1-cos(2*dpi*(i)/_length))); // dan made this change to do rewindowing
lbajardsilogic@79: 		--i;
lbajardsilogic@79: 	}
lbajardsilogic@79: 	while (i >= 0);
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*\\\ NESTED CLASS MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: Constructor                                                   */
lbajardsilogic@79: /*      Input parameters:                                                   */
lbajardsilogic@79: /*        - nbr_bits: number of bits of the array on which we want to do a  */
lbajardsilogic@79: /*                    FFT. Range: > 0                                       */
lbajardsilogic@79: /*      Throws: std::bad_alloc                                              */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: FFTReal::BitReversedLUT::BitReversedLUT (const int nbr_bits)
lbajardsilogic@79: {
lbajardsilogic@79: 	long				length;
lbajardsilogic@79: 	long				cnt;
lbajardsilogic@79: 	long				br_index;
lbajardsilogic@79: 	long				bit;
lbajardsilogic@79: 
lbajardsilogic@79: 	length = 1L << nbr_bits;
lbajardsilogic@79: 	_ptr = new long [length];
lbajardsilogic@79: 
lbajardsilogic@79: 	br_index = 0;
lbajardsilogic@79: 	_ptr [0] = 0;
lbajardsilogic@79: 	for (cnt = 1; cnt < length; ++cnt)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		/* ++br_index (bit reversed) */
lbajardsilogic@79: 		bit = length >> 1;
lbajardsilogic@79: 		while (((br_index ^= bit) & bit) == 0)
lbajardsilogic@79: 		{
lbajardsilogic@79: 			bit >>= 1;
lbajardsilogic@79: 		}
lbajardsilogic@79: 
lbajardsilogic@79: 		_ptr [cnt] = br_index;
lbajardsilogic@79: 	}
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: Destructor                                                    */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: FFTReal::BitReversedLUT::~BitReversedLUT (void)
lbajardsilogic@79: {
lbajardsilogic@79: 	delete [] _ptr;
lbajardsilogic@79: 	_ptr = 0;
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: Constructor                                                   */
lbajardsilogic@79: /*      Input parameters:                                                   */
lbajardsilogic@79: /*        - nbr_bits: number of bits of the array on which we want to do a  */
lbajardsilogic@79: /*                    FFT. Range: > 0                                       */
lbajardsilogic@79: /*      Throws: std::bad_alloc, anything                                    */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: FFTReal::TrigoLUT::TrigoLUT (const int nbr_bits)
lbajardsilogic@79: {
lbajardsilogic@79: 	long		total_len;
lbajardsilogic@79: 
lbajardsilogic@79: 	_ptr = 0;
lbajardsilogic@79: 	if (nbr_bits > 3)
lbajardsilogic@79: 	{
lbajardsilogic@79: 		total_len = (1L << (nbr_bits - 1)) - 4;
lbajardsilogic@79: 		_ptr = new flt_t [total_len];
lbajardsilogic@79: 
lbajardsilogic@79: 		const double	PI = atan (1.0f) * 4;
lbajardsilogic@79: 		for (int level = 3; level < nbr_bits; ++level)
lbajardsilogic@79: 		{
lbajardsilogic@79: 			const long		level_len = 1L << (level - 1);
lbajardsilogic@79: 			flt_t	* const	level_ptr = const_cast<flt_t *> (get_ptr (level));
lbajardsilogic@79: 			const double	mul = PI / (level_len << 1);
lbajardsilogic@79: 
lbajardsilogic@79: 			for (long i = 0; i < level_len; ++ i)
lbajardsilogic@79: 			{
lbajardsilogic@79: 				level_ptr [i] = (flt_t) cos (i * mul);
lbajardsilogic@79: 			}
lbajardsilogic@79: 		}
lbajardsilogic@79: 	}
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: /*      Name: Destructor                                                    */
lbajardsilogic@79: /*==========================================================================*/
lbajardsilogic@79: 
lbajardsilogic@79: FFTReal::TrigoLUT::~TrigoLUT (void)
lbajardsilogic@79: {
lbajardsilogic@79: 	delete [] _ptr;
lbajardsilogic@79: 	_ptr = 0;
lbajardsilogic@79: }
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: #if defined (_MSC_VER)
lbajardsilogic@79: #pragma pack (pop)
lbajardsilogic@79: #endif	// _MSC_VER
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*****************************************************************************
lbajardsilogic@79: 
lbajardsilogic@79: 	LEGAL
lbajardsilogic@79: 
lbajardsilogic@79: 	Source code may be freely used for any purpose, including commercial
lbajardsilogic@79: 	applications. Programs must display in their "About" dialog-box (or
lbajardsilogic@79: 	documentation) a text telling they use these routines by Laurent de Soras.
lbajardsilogic@79: 	Modified source code can be distributed, but modifications must be clearly
lbajardsilogic@79: 	indicated.
lbajardsilogic@79: 
lbajardsilogic@79: 	CONTACT
lbajardsilogic@79: 
lbajardsilogic@79: 	Laurent de Soras
lbajardsilogic@79: 	92 avenue Albert 1er
lbajardsilogic@79: 	92500 Rueil-Malmaison
lbajardsilogic@79: 	France
lbajardsilogic@79: 
lbajardsilogic@79: 	ldesoras@club-internet.fr
lbajardsilogic@79: 
lbajardsilogic@79: *****************************************************************************/
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: 
lbajardsilogic@79: /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
lbajardsilogic@79: