fftreal.cpp

FFT源码

/*****************************************************************************
*                                                                            *
*       DIGITAL SIGNAL PROCESSING TOOLS                                      *
*       Version 1.01, 1999/11/07                                             *
*       (c) 1999 - Laurent de Soras                                          *
*                                                                            *
*       FFTReal.cpp                                                          *
*       Fourier transformation of real number arrays.                        *
*       Portable ISO C++                                                     *
*                                                                            *
* Tab = 3                                                                    *
*****************************************************************************/



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

#include	"FFTReal.h"

#include	<cassert>
#include	<cmath>



#if defined (_MSC_VER)
#pragma pack (push, 8)
#endif	// _MSC_VER



/*\\\ PUBLIC MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/



/*==========================================================================*/
/*      Name: Constructor                                                   */
/*      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, anything                                    */
/*==========================================================================*/

FFTReal::FFTReal (const long length)
:	_length (length)
,	_nbr_bits (int (floor (log (length) / log (2) + 0.5)))
,	_bit_rev_lut (int (floor (log (length) / log (2) + 0.5)))
,	_trigo_lut (int (floor (log (length) / log (2) + 0.5)))
,	_sqrt2_2 (flt_t (sqrt (2) * 0.5))
{
	assert ((1L << _nbr_bits) == length);

	_buffer_ptr = 0;
	if (_nbr_bits > 2)
	{
		_buffer_ptr = new flt_t [_length];
	}
}



/*==========================================================================*/
/*      Name: Destructor                                                    */
/*==========================================================================*/

FFTReal::~FFTReal (void)
{
	delete [] _buffer_ptr;
	_buffer_ptr = 0;
}



/*==========================================================================*/
/*      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] = imaginary values of        */
/*               coefficents 1...length(x)/2-1.                             */
/*      Throws: Nothing                                                     */
/*==========================================================================*/

void	FFTReal::do_fft (flt_t f [], const flt_t x []) const
{

/*______________________________________________
 *
 * General case
 *______________________________________________
 */

	if (_nbr_bits > 2)
	{
		flt_t *			sf;
		flt_t *			df;

		if (_nbr_bits & 1)
		{
			df = _buffer_ptr;
			sf = f;
		}
		else
		{
			df = f;
			sf = _buffer_ptr;
		}

		/* Do the transformation in several pass */
		{
			int		pass;
			long		nbr_coef;
			long		h_nbr_coef;
			long		d_nbr_coef;
			long		coef_index;

			/* First and second pass at once */
			{
				const long * const	bit_rev_lut_ptr = _bit_rev_lut.get_ptr ();
				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];

					flt_t	* 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 flt_t		sf_0 = x [rev_index_0] + x [rev_index_1];
					const flt_t		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);
			}

			/* Third pass */
			{
				coef_index = 0;
				const flt_t		sqrt2_2 = _sqrt2_2;
				do
				{
					flt_t				v;

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

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

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

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

			/* Next pass */
			for (pass = 3; pass < _nbr_bits; ++pass)
			{
				coef_index = 0;
				nbr_coef = 1 << pass;
				h_nbr_coef = nbr_coef >> 1;
				d_nbr_coef = nbr_coef << 1;
				const flt_t	* const	cos_ptr = _trigo_lut.get_ptr (pass);
				do
				{
					long				i;
					const flt_t	*	const sf1r = sf + coef_index;
					const flt_t	*	const sf2r = sf1r + nbr_coef;
					flt_t *			const dfr = df + coef_index;
					flt_t *			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 flt_t	* const	sf1i = sf1r + h_nbr_coef;
					const flt_t	* const	sf2i = sf1i + nbr_coef;
					for (i = 1; i < h_nbr_coef; ++ i)
					{
						const flt_t		c = cos_ptr [i];					// cos (i*PI/nbr_coef);
						const flt_t		s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);
						flt_t				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);

				/* Prepare to the next pass */
				{
					flt_t	* const		temp_ptr = df;
					df = sf;
					sf = temp_ptr;
				}
			}
		}
	}

/*______________________________________________
 *
 * Special cases
 *______________________________________________
 */

	/* 4-point FFT */
	else if (_nbr_bits == 2)
	{
		f [1] = x [0] - x [2];
		f [3] = x [1] - x [3];

		const flt_t			b_0 = x [0] + x [2];
		const flt_t			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. Notice that      */
/*                   IFFT (FFT (x)) = x * length (x). Data must be          */
/*                   post-scaled.                                           */
/*      Input parameters:                                                   */
/*        - f: pointer on the source array (frequencies).                   */
/*             f [0...length(x)/2] = real values,                           */
/*             f [length(x)/2+1...length(x)] = imaginary values of          */
/*               coefficents 1...length(x)-1.                               */
/*      Output parameters:                                                  */
/*        - x: pointer on the destination array (time).                     */
/*      Throws: Nothing                                                     */
/*==========================================================================*/

void	FFTReal::do_ifft (const flt_t f [], flt_t x []) const
{

/*______________________________________________
 *
 * General case
 *______________________________________________
 */

	if (_nbr_bits > 2)
	{
		flt_t *			sf = const_cast <flt_t *> (f);
		flt_t *			df;
		flt_t *			df_temp;

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

		/* Do the transformation in several pass */
		{
			int			pass;
			long			nbr_coef;
			long			h_nbr_coef;
			long			d_nbr_coef;
			long			coef_index;