fftreal.cpp

FFT源码

			/* First pass */
			for (pass = _nbr_bits - 1; pass >= 3; --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 sfr = sf + coef_index;
					const flt_t	*	const sfi = sfr + nbr_coef;
					flt_t *			const df1r = df + coef_index;
					flt_t *			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 */
					flt_t * const	df1i = df1r + h_nbr_coef;
					flt_t * const	df2i = df1i + nbr_coef;
					for (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 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);
						const flt_t		vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]
						const flt_t		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);

				/* Prepare to the next pass */
				if (pass < _nbr_bits - 1)
				{
					flt_t	* const	temp_ptr = df;
					df = sf;
					sf = temp_ptr;
				}
				else
				{
					sf = df;
					df = df_temp;
				}
			}

			/* Antepenultimate pass */
			{
				const flt_t		sqrt2_2 = _sqrt2_2;
				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 flt_t		vr = sf [coef_index + 1] - sf [coef_index + 3];
					const flt_t		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);
			}

			/* Penultimate and last pass at once */
			{
				coef_index = 0;
				const long *	bit_rev_lut_ptr = _bit_rev_lut.get_ptr ();
				const flt_t	*	sf2 = df;
				do
				{
					{
						const flt_t		b_0 = sf2 [0] + sf2 [2];
						const flt_t		b_2 = sf2 [0] - sf2 [2];
						const flt_t		b_1 = sf2 [1] * 2;
						const flt_t		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 flt_t		b_0 = sf2 [4] + sf2 [6];
						const flt_t		b_2 = sf2 [4] - sf2 [6];
						const flt_t		b_1 = sf2 [5] * 2;
						const flt_t		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);
			}
		}
	}

/*______________________________________________
 *
 * Special cases
 *______________________________________________
 */

	/* 4-point IFFT */
	else if (_nbr_bits == 2)
	{
		const flt_t		b_0 = f [0] + f [2];
		const flt_t		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/Output parameters:                                            */
/*        - x: pointer on array to rescale (time or frequency).             */
/*      Throws: Nothing                                                     */
/*==========================================================================*/

void	FFTReal::rescale (flt_t x []) const
{
	const flt_t		mul = flt_t (1.0 / _length);
	long				i = _length - 1;

	do
	{
		x [i] *= mul;
		--i;
	}
	while (i >= 0);
}



/*\\\ NESTED CLASS MEMBER FUNCTIONS \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/



/*==========================================================================*/
/*      Name: Constructor                                                   */
/*      Input parameters:                                                   */
/*        - nbr_bits: number of bits of the array on which we want to do a  */
/*                    FFT. Range: > 0                                       */
/*      Throws: std::bad_alloc                                              */
/*==========================================================================*/

FFTReal::BitReversedLUT::BitReversedLUT (const int nbr_bits)
{
	long				length;
	long				cnt;
	long				br_index;
	long				bit;

	length = 1L << nbr_bits;
	_ptr = new long [length];

	br_index = 0;
	_ptr [0] = 0;
	for (cnt = 1; cnt < length; ++cnt)
	{
		/* ++br_index (bit reversed) */
		bit = length >> 1;
		while (((br_index ^= bit) & bit) == 0)
		{
			bit >>= 1;
		}

		_ptr [cnt] = br_index;
	}
}



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

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



/*==========================================================================*/
/*      Name: Constructor                                                   */
/*      Input parameters:                                                   */
/*        - nbr_bits: number of bits of the array on which we want to do a  */
/*                    FFT. Range: > 0                                       */
/*      Throws: std::bad_alloc, anything                                    */
/*==========================================================================*/

FFTReal::TrigoLUT::TrigoLUT (const int nbr_bits)
{
	long		total_len;

	_ptr = 0;
	if (nbr_bits > 3)
	{
		total_len = (1L << (nbr_bits - 1)) - 4;
		_ptr = new flt_t [total_len];

		const double	PI = atan (1) * 4;
		for (int level = 3; level < nbr_bits; ++level)
		{
			const long		level_len = 1L << (level - 1);
			flt_t	* const	level_ptr = const_cast<flt_t *> (get_ptr (level));
			const double	mul = PI / (level_len << 1);

			for (long i = 0; i < level_len; ++ i)
			{
				level_ptr [i] = (flt_t) cos (i * mul);
			}
		}
	}
}



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

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



#if defined (_MSC_VER)
#pragma pack (pop)
#endif	// _MSC_VER



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

	LEGAL

	Source code may be freely used for any purpose, including commercial
	applications. Programs must display in their "About" dialog-box (or
	documentation) a text telling they use these routines by Laurent de Soras.
	Modified source code can be distributed, but modifications must be clearly
	indicated.

	CONTACT

	Laurent de Soras
	92 avenue Albert 1er
	92500 Rueil-Malmaison
	France

	[e-mail]

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



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