1 /*****************************************************************************

     2 

     3         FFTReal.hpp

     4         Copyright (c) 2005 Laurent de Soras

     5 

     6 --- Legal stuff ---

     7 

     8 This library is free software; you can redistribute it and/or

     9 modify it under the terms of the GNU Lesser General Public

    10 License as published by the Free Software Foundation; either

    11 version 2.1 of the License, or (at your option) any later version.

    12 

    13 This library is distributed in the hope that it will be useful,

    14 but WITHOUT ANY WARRANTY; without even the implied warranty of

    15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

    16 Lesser General Public License for more details.

    17 

    18 You should have received a copy of the GNU Lesser General Public

    19 License along with this library; if not, write to the Free Software

    20 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    21 

    22 *Tab=3***********************************************************************/

    23 

    24 

    25 

    26 #if defined (FFTReal_CURRENT_CODEHEADER)

    27 	#error Recursive inclusion of FFTReal code header.

    28 #endif

    29 #define	FFTReal_CURRENT_CODEHEADER

    30 

    31 #if ! defined (FFTReal_CODEHEADER_INCLUDED)

    32 #define	FFTReal_CODEHEADER_INCLUDED

    33 

    34 

    35 

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

    37 

    38 #include	<cassert>

    39 #include	<cmath>

    40 

    41 

    42 

    43 static inline bool	FFTReal_is_pow2 (long x)

    44 {

    45 	assert (x > 0);

    46 

    47 	return  ((x & -x) == x);

    48 }

    49 

    50 

    51 

    52 static inline int	FFTReal_get_next_pow2 (long x)

    53 {

    54 	--x;

    55 

    56 	int				p = 0;

    57 	while ((x & ~0xFFFFL) != 0)

    58 	{

    59 		p += 16;

    60 		x >>= 16;

    61 	}

    62 	while ((x & ~0xFL) != 0)

    63 	{

    64 		p += 4;

    65 		x >>= 4;

    66 	}

    67 	while (x > 0)

    68 	{

    69 		++p;

    70 		x >>= 1;

    71 	}

    72 

    73 	return (p);

    74 }

    75 

    76 

    77 

    78 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

    79 

    80 

    81 

    82 /*

    83 ==============================================================================

    84 Name: ctor

    85 Input parameters:

    86 	- length: length of the array on which we want to do a FFT. Range: power of

    87 		2 only, > 0.

    88 Throws: std::bad_alloc

    89 ==============================================================================

    90 */

    91 

    92 template <class DT>

    93 FFTReal <DT>::FFTReal (long length)

    94 :	_length (length)

    95 ,	_nbr_bits (FFTReal_get_next_pow2 (length))

    96 ,	_br_lut ()

    97 ,	_trigo_lut ()

    98 ,	_buffer (length)

    99 ,	_trigo_osc ()

   100 {

   101 	assert (FFTReal_is_pow2 (length));

   102 	assert (_nbr_bits <= MAX_BIT_DEPTH);

   103 

   104 	init_br_lut ();

   105 	init_trigo_lut ();

   106 	init_trigo_osc ();

   107 }

   108 

   109 

   110 

   111 /*

   112 ==============================================================================

   113 Name: get_length

   114 Description:

   115 	Returns the number of points processed by this FFT object.

   116 Returns: The number of points, power of 2, > 0.

   117 Throws: Nothing

   118 ==============================================================================

   119 */

   120 

   121 template <class DT>

   122 long	FFTReal <DT>::get_length () const

   123 {

   124 	return (_length);

   125 }

   126 

   127 

   128 

   129 /*

   130 ==============================================================================

   131 Name: do_fft

   132 Description:

   133 	Compute the FFT of the array.

   134 Input parameters:

   135 	- x: pointer on the source array (time).

   136 Output parameters:

   137 	- f: pointer on the destination array (frequencies).

   138 		f [0...length(x)/2] = real values,

   139 		f [length(x)/2+1...length(x)-1] = negative imaginary values of

   140 		coefficents 1...length(x)/2-1.

   141 Throws: Nothing

   142 ==============================================================================

   143 */

   144 

   145 template <class DT>

   146 void	FFTReal <DT>::do_fft (DataType f [], const DataType x []) const

   147 {

   148 	assert (f != 0);

   149 	assert (f != use_buffer ());

   150 	assert (x != 0);

   151 	assert (x != use_buffer ());

   152 	assert (x != f);

   153 

   154 	// General case

   155 	if (_nbr_bits > 2)

   156 	{

   157 		compute_fft_general (f, x);

   158 	}

   159 

   160 	// 4-point FFT

   161 	else if (_nbr_bits == 2)

   162 	{

   163 		f [1] = x [0] - x [2];

   164 		f [3] = x [1] - x [3];

   165 

   166 		const DataType	b_0 = x [0] + x [2];

   167 		const DataType	b_2 = x [1] + x [3];

   168 		

   169 		f [0] = b_0 + b_2;

   170 		f [2] = b_0 - b_2;

   171 	}

   172 

   173 	// 2-point FFT

   174 	else if (_nbr_bits == 1)

   175 	{

   176 		f [0] = x [0] + x [1];

   177 		f [1] = x [0] - x [1];

   178 	}

   179 

   180 	// 1-point FFT

   181 	else

   182 	{

   183 		f [0] = x [0];

   184 	}

   185 }

   186 

   187 

   188 

   189 /*

   190 ==============================================================================

   191 Name: do_ifft

   192 Description:

   193 	Compute the inverse FFT of the array. Note that data must be post-scaled:

   194 	IFFT (FFT (x)) = x * length (x).

   195 Input parameters:

   196 	- f: pointer on the source array (frequencies).

   197 		f [0...length(x)/2] = real values

   198 		f [length(x)/2+1...length(x)-1] = negative imaginary values of

   199 		coefficents 1...length(x)/2-1.

   200 Output parameters:

   201 	- x: pointer on the destination array (time).

   202 Throws: Nothing

   203 ==============================================================================

   204 */

   205 

   206 template <class DT>

   207 void	FFTReal <DT>::do_ifft (const DataType f [], DataType x []) const

   208 {

   209 	assert (f != 0);

   210 	assert (f != use_buffer ());

   211 	assert (x != 0);

   212 	assert (x != use_buffer ());

   213 	assert (x != f);

   214 

   215 	// General case

   216 	if (_nbr_bits > 2)

   217 	{

   218 		compute_ifft_general (f, x);

   219 	}

   220 

   221 	// 4-point IFFT

   222 	else if (_nbr_bits == 2)

   223 	{

   224 		const DataType	b_0 = f [0] + f [2];

   225 		const DataType	b_2 = f [0] - f [2];

   226 

   227 		x [0] = b_0 + f [1] * 2;

   228 		x [2] = b_0 - f [1] * 2;

   229 		x [1] = b_2 + f [3] * 2;

   230 		x [3] = b_2 - f [3] * 2;

   231 	}

   232 

   233 	// 2-point IFFT

   234 	else if (_nbr_bits == 1)

   235 	{

   236 		x [0] = f [0] + f [1];

   237 		x [1] = f [0] - f [1];

   238 	}

   239 

   240 	// 1-point IFFT

   241 	else

   242 	{

   243 		x [0] = f [0];

   244 	}

   245 }

   246 

   247 

   248 

   249 /*

   250 ==============================================================================

   251 Name: rescale

   252 Description:

   253 	Scale an array by divide each element by its length. This function should

   254 	be called after FFT + IFFT.

   255 Input parameters:

   256 	- x: pointer on array to rescale (time or frequency).

   257 Throws: Nothing

   258 ==============================================================================

   259 */

   260 

   261 template <class DT>

   262 void	FFTReal <DT>::rescale (DataType x []) const

   263 {

   264 	const DataType	mul = DataType (1.0 / _length);

   265 

   266 	if (_length < 4)

   267 	{

   268 		long				i = _length - 1;

   269 		do

   270 		{

   271 			x [i] *= mul;

   272 			--i;

   273 		}

   274 		while (i >= 0);

   275 	}

   276 

   277 	else

   278 	{

   279 		assert ((_length & 3) == 0);

   280 

   281 		// Could be optimized with SIMD instruction sets (needs alignment check)

   282 		long				i = _length - 4;

   283 		do

   284 		{

   285 			x [i + 0] *= mul;

   286 			x [i + 1] *= mul;

   287 			x [i + 2] *= mul;

   288 			x [i + 3] *= mul;

   289 			i -= 4;

   290 		}

   291 		while (i >= 0);

   292 	}

   293 }

   294 

   295 

   296 

   297 /*

   298 ==============================================================================

   299 Name: use_buffer

   300 Description:

   301 	Access the internal buffer, whose length is the FFT one.

   302 	Buffer content will be erased at each do_fft() / do_ifft() call!

   303 	This buffer cannot be used as:

   304 		- source for FFT or IFFT done with this object

   305 		- destination for FFT or IFFT done with this object

   306 Returns:

   307 	Buffer start address

   308 Throws: Nothing

   309 ==============================================================================

   310 */

   311 

   312 template <class DT>

   313 typename FFTReal <DT>::DataType *	FFTReal <DT>::use_buffer () const

   314 {

   315 	return (&_buffer [0]);

   316 }

   317 

   318 

   319 

   320 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

   321 

   322 

   323 

   324 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/

   325 

   326 

   327 

   328 template <class DT>

   329 void	FFTReal <DT>::init_br_lut ()

   330 {

   331 	const long		length = 1L << _nbr_bits;

   332 	_br_lut.resize (length);

   333 

   334 	_br_lut [0] = 0;

   335 	long				br_index = 0;

   336 	for (long cnt = 1; cnt < length; ++cnt)

   337 	{

   338 		// ++br_index (bit reversed)

   339 		long				bit = length >> 1;

   340 		while (((br_index ^= bit) & bit) == 0)

   341 		{

   342 			bit >>= 1;

   343 		}

   344 

   345 		_br_lut [cnt] = br_index;

   346 	}

   347 }

   348 

   349 

   350 

   351 template <class DT>

   352 void	FFTReal <DT>::init_trigo_lut ()

   353 {

   354 	using namespace std;

   355 

   356 	if (_nbr_bits > 3)

   357 	{

   358 		const long		total_len = (1L << (_nbr_bits - 1)) - 4;

   359 		_trigo_lut.resize (total_len);

   360 

   361 		for (int level = 3; level < _nbr_bits; ++level)

   362 		{

   363 			const long		level_len = 1L << (level - 1);

   364 			DataType	* const	level_ptr =

   365 				&_trigo_lut [get_trigo_level_index (level)];

   366 			const double	mul = PI / (level_len << 1);

   367 

   368 			for (long i = 0; i < level_len; ++ i)

   369 			{

   370 				level_ptr [i] = static_cast <DataType> (cos (i * mul));

   371 			}

   372 		}

   373 	}

   374 }

   375 

   376 

   377 

   378 template <class DT>

   379 void	FFTReal <DT>::init_trigo_osc ()

   380 {

   381 	const int		nbr_osc = _nbr_bits - TRIGO_BD_LIMIT;

   382 	if (nbr_osc > 0)

   383 	{

   384 		_trigo_osc.resize (nbr_osc);

   385 

   386 		for (int osc_cnt = 0; osc_cnt < nbr_osc; ++osc_cnt)

   387 		{

   388 			OscType &		osc = _trigo_osc [osc_cnt];

   389 

   390 			const long		len = 1L << (TRIGO_BD_LIMIT + osc_cnt);

   391 			const double	mul = (0.5 * PI) / len;

   392 			osc.set_step (mul);

   393 		}

   394 	}

   395 }

   396 

   397 

   398 

   399 template <class DT>

   400 const long *	FFTReal <DT>::get_br_ptr () const

   401 {

   402 	return (&_br_lut [0]);

   403 }

   404 

   405 

   406 

   407 template <class DT>

   408 const typename FFTReal <DT>::DataType *	FFTReal <DT>::get_trigo_ptr (int level) const

   409 {

   410 	assert (level >= 3);

   411 

   412 	return (&_trigo_lut [get_trigo_level_index (level)]);

   413 }

   414 

   415 

   416 

   417 template <class DT>

   418 long	FFTReal <DT>::get_trigo_level_index (int level) const

   419 {

   420 	assert (level >= 3);

   421 

   422 	return ((1L << (level - 1)) - 4);

   423 }

   424 

   425 

   426 

   427 // Transform in several passes

   428 template <class DT>

   429 void	FFTReal <DT>::compute_fft_general (DataType f [], const DataType x []) const

   430 {

   431 	assert (f != 0);

   432 	assert (f != use_buffer ());

   433 	assert (x != 0);

   434 	assert (x != use_buffer ());

   435 	assert (x != f);

   436 

   437 	DataType *		sf;

   438 	DataType *		df;

   439 

   440 	if ((_nbr_bits & 1) != 0)

   441 	{

   442 		df = use_buffer ();

   443 		sf = f;

   444 	}

   445 	else

   446 	{

   447 		df = f;

   448 		sf = use_buffer ();

   449 	}

   450 

   451 	compute_direct_pass_1_2 (df, x);

   452 	compute_direct_pass_3 (sf, df);

   453 

   454 	for (int pass = 3; pass < _nbr_bits; ++ pass)

   455 	{

   456 		compute_direct_pass_n (df, sf, pass);

   457 

   458 		DataType * const	temp_ptr = df;

   459 		df = sf;

   460 		sf = temp_ptr;

   461 	}

   462 }

   463 

   464 

   465 

   466 template <class DT>

   467 void	FFTReal <DT>::compute_direct_pass_1_2 (DataType df [], const DataType x []) const

   468 {

   469 	assert (df != 0);

   470 	assert (x != 0);

   471 	assert (df != x);

   472 

   473 	const long * const	bit_rev_lut_ptr = get_br_ptr ();

   474 	long				coef_index = 0;

   475 	do

   476 	{

   477 		const long		rev_index_0 = bit_rev_lut_ptr [coef_index];

   478 		const long		rev_index_1 = bit_rev_lut_ptr [coef_index + 1];

   479 		const long		rev_index_2 = bit_rev_lut_ptr [coef_index + 2];

   480 		const long		rev_index_3 = bit_rev_lut_ptr [coef_index + 3];

   481 

   482 		DataType	* const	df2 = df + coef_index;

   483 		df2 [1] = x [rev_index_0] - x [rev_index_1];

   484 		df2 [3] = x [rev_index_2] - x [rev_index_3];

   485 

   486 		const DataType	sf_0 = x [rev_index_0] + x [rev_index_1];

   487 		const DataType	sf_2 = x [rev_index_2] + x [rev_index_3];

   488 

   489 		df2 [0] = sf_0 + sf_2;

   490 		df2 [2] = sf_0 - sf_2;

   491 		

   492 		coef_index += 4;

   493 	}

   494 	while (coef_index < _length);

   495 }

   496 

   497 

   498 

   499 template <class DT>

   500 void	FFTReal <DT>::compute_direct_pass_3 (DataType df [], const DataType sf []) const

   501 {

   502 	assert (df != 0);

   503 	assert (sf != 0);

   504 	assert (df != sf);

   505 

   506 	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);

   507 	long				coef_index = 0;

   508 	do

   509 	{

   510 		DataType			v;

   511 

   512 		df [coef_index] = sf [coef_index] + sf [coef_index + 4];

   513 		df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];

   514 		df [coef_index + 2] = sf [coef_index + 2];

   515 		df [coef_index + 6] = sf [coef_index + 6];

   516 

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

   518 		df [coef_index + 1] = sf [coef_index + 1] + v;

   519 		df [coef_index + 3] = sf [coef_index + 1] - v;

   520 

   521 		v = (sf [coef_index + 5] + sf [coef_index + 7]) * sqrt2_2;

   522 		df [coef_index + 5] = v + sf [coef_index + 3];

   523 		df [coef_index + 7] = v - sf [coef_index + 3];

   524 

   525 		coef_index += 8;

   526 	}

   527 	while (coef_index < _length);

   528 }

   529 

   530 

   531 

   532 template <class DT>

   533 void	FFTReal <DT>::compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const

   534 {

   535 	assert (df != 0);

   536 	assert (sf != 0);

   537 	assert (df != sf);

   538 	assert (pass >= 3);

   539 	assert (pass < _nbr_bits);

   540 

   541 	if (pass <= TRIGO_BD_LIMIT)

   542 	{

   543 		compute_direct_pass_n_lut (df, sf, pass);

   544 	}

   545 	else

   546 	{

   547 		compute_direct_pass_n_osc (df, sf, pass);

   548 	}

   549 }

   550 

   551 

   552 

   553 template <class DT>

   554 void	FFTReal <DT>::compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const

   555 {

   556 	assert (df != 0);

   557 	assert (sf != 0);

   558 	assert (df != sf);

   559 	assert (pass >= 3);

   560 	assert (pass < _nbr_bits);

   561 

   562 	const long		nbr_coef = 1 << pass;

   563 	const long		h_nbr_coef = nbr_coef >> 1;

   564 	const long		d_nbr_coef = nbr_coef << 1;

   565 	long				coef_index = 0;

   566 	const DataType	* const	cos_ptr = get_trigo_ptr (pass);

   567 	do

   568 	{

   569 		const DataType	* const	sf1r = sf + coef_index;

   570 		const DataType	* const	sf2r = sf1r + nbr_coef;

   571 		DataType			* const	dfr = df + coef_index;

   572 		DataType			* const	dfi = dfr + nbr_coef;

   573 

   574 		// Extreme coefficients are always real

   575 		dfr [0] = sf1r [0] + sf2r [0];

   576 		dfi [0] = sf1r [0] - sf2r [0];	// dfr [nbr_coef] =

   577 		dfr [h_nbr_coef] = sf1r [h_nbr_coef];

   578 		dfi [h_nbr_coef] = sf2r [h_nbr_coef];

   579 

   580 		// Others are conjugate complex numbers

   581 		const DataType * const	sf1i = sf1r + h_nbr_coef;

   582 		const DataType * const	sf2i = sf1i + nbr_coef;

   583 		for (long i = 1; i < h_nbr_coef; ++ i)

   584 		{

   585 			const DataType	c = cos_ptr [i];					// cos (i*PI/nbr_coef);

   586 			const DataType	s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);

   587 			DataType	 		v;

   588 

   589 			v = sf2r [i] * c - sf2i [i] * s;

   590 			dfr [i] = sf1r [i] + v;

   591 			dfi [-i] = sf1r [i] - v;	// dfr [nbr_coef - i] =

   592 

   593 			v = sf2r [i] * s + sf2i [i] * c;

   594 			dfi [i] = v + sf1i [i];

   595 			dfi [nbr_coef - i] = v - sf1i [i];

   596 		}

   597 

   598 		coef_index += d_nbr_coef;

   599 	}

   600 	while (coef_index < _length);

   601 }

   602 

   603 

   604 

   605 template <class DT>

   606 void	FFTReal <DT>::compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const

   607 {

   608 	assert (df != 0);

   609 	assert (sf != 0);

   610 	assert (df != sf);

   611 	assert (pass > TRIGO_BD_LIMIT);

   612 	assert (pass < _nbr_bits);

   613 

   614 	const long		nbr_coef = 1 << pass;

   615 	const long		h_nbr_coef = nbr_coef >> 1;

   616 	const long		d_nbr_coef = nbr_coef << 1;

   617 	long				coef_index = 0;

   618 	OscType &		osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];

   619 	do

   620 	{

   621 		const DataType	* const	sf1r = sf + coef_index;

   622 		const DataType	* const	sf2r = sf1r + nbr_coef;

   623 		DataType			* const	dfr = df + coef_index;

   624 		DataType			* const	dfi = dfr + nbr_coef;

   625 

   626 		osc.clear_buffers ();

   627 

   628 		// Extreme coefficients are always real

   629 		dfr [0] = sf1r [0] + sf2r [0];

   630 		dfi [0] = sf1r [0] - sf2r [0];	// dfr [nbr_coef] =

   631 		dfr [h_nbr_coef] = sf1r [h_nbr_coef];

   632 		dfi [h_nbr_coef] = sf2r [h_nbr_coef];

   633 

   634 		// Others are conjugate complex numbers

   635 		const DataType * const	sf1i = sf1r + h_nbr_coef;

   636 		const DataType * const	sf2i = sf1i + nbr_coef;

   637 		for (long i = 1; i < h_nbr_coef; ++ i)

   638 		{

   639 			osc.step ();

   640 			const DataType	c = osc.get_cos ();

   641 			const DataType	s = osc.get_sin ();

   642 			DataType	 		v;

   643 

   644 			v = sf2r [i] * c - sf2i [i] * s;

   645 			dfr [i] = sf1r [i] + v;

   646 			dfi [-i] = sf1r [i] - v;	// dfr [nbr_coef - i] =

   647 

   648 			v = sf2r [i] * s + sf2i [i] * c;

   649 			dfi [i] = v + sf1i [i];

   650 			dfi [nbr_coef - i] = v - sf1i [i];

   651 		}

   652 

   653 		coef_index += d_nbr_coef;

   654 	}

   655 	while (coef_index < _length);

   656 }

   657 

   658 

   659 

   660 // Transform in several pass

   661 template <class DT>

   662 void	FFTReal <DT>::compute_ifft_general (const DataType f [], DataType x []) const

   663 {

   664 	assert (f != 0);

   665 	assert (f != use_buffer ());

   666 	assert (x != 0);

   667 	assert (x != use_buffer ());

   668 	assert (x != f);

   669 

   670 	DataType *		sf = const_cast <DataType *> (f);

   671 	DataType *		df;

   672 	DataType *		df_temp;

   673 

   674 	if (_nbr_bits & 1)

   675 	{

   676 		df = use_buffer ();

   677 		df_temp = x;

   678 	}

   679 	else

   680 	{

   681 		df = x;

   682 		df_temp = use_buffer ();

   683 	}

   684 

   685 	for (int pass = _nbr_bits - 1; pass >= 3; -- pass)

   686 	{

   687 		compute_inverse_pass_n (df, sf, pass);

   688 

   689 		if (pass < _nbr_bits - 1)

   690 		{

   691 			DataType	* const	temp_ptr = df;

   692 			df = sf;

   693 			sf = temp_ptr;

   694 		}

   695 		else

   696 		{

   697 			sf = df;

   698 			df = df_temp;

   699 		}

   700 	}

   701 

   702 	compute_inverse_pass_3 (df, sf);

   703 	compute_inverse_pass_1_2 (x, df);

   704 }

   705 

   706 

   707 

   708 template <class DT>

   709 void	FFTReal <DT>::compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const

   710 {

   711 	assert (df != 0);

   712 	assert (sf != 0);

   713 	assert (df != sf);

   714 	assert (pass >= 3);

   715 	assert (pass < _nbr_bits);

   716 

   717 	if (pass <= TRIGO_BD_LIMIT)

   718 	{

   719 		compute_inverse_pass_n_lut (df, sf, pass);

   720 	}

   721 	else

   722 	{

   723 		compute_inverse_pass_n_osc (df, sf, pass);

   724 	}

   725 }

   726 

   727 

   728 

   729 template <class DT>

   730 void	FFTReal <DT>::compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const

   731 {

   732 	assert (df != 0);

   733 	assert (sf != 0);

   734 	assert (df != sf);

   735 	assert (pass >= 3);

   736 	assert (pass < _nbr_bits);

   737 

   738 	const long		nbr_coef = 1 << pass;

   739 	const long		h_nbr_coef = nbr_coef >> 1;

   740 	const long		d_nbr_coef = nbr_coef << 1;

   741 	long				coef_index = 0;

   742 	const DataType * const	cos_ptr = get_trigo_ptr (pass);

   743 	do

   744 	{

   745 		const DataType	* const	sfr = sf + coef_index;

   746 		const DataType	* const	sfi = sfr + nbr_coef;

   747 		DataType			* const	df1r = df + coef_index;

   748 		DataType			* const	df2r = df1r + nbr_coef;

   749 

   750 		// Extreme coefficients are always real

   751 		df1r [0] = sfr [0] + sfi [0];		// + sfr [nbr_coef]

   752 		df2r [0] = sfr [0] - sfi [0];		// - sfr [nbr_coef]

   753 		df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;

   754 		df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;

   755 

   756 		// Others are conjugate complex numbers

   757 		DataType * const	df1i = df1r + h_nbr_coef;

   758 		DataType * const	df2i = df1i + nbr_coef;

   759 		for (long i = 1; i < h_nbr_coef; ++ i)

   760 		{

   761 			df1r [i] = sfr [i] + sfi [-i];		// + sfr [nbr_coef - i]

   762 			df1i [i] = sfi [i] - sfi [nbr_coef - i];

   763 

   764 			const DataType	c = cos_ptr [i];					// cos (i*PI/nbr_coef);

   765 			const DataType	s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);

   766 			const DataType	vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]

   767 			const DataType	vi = sfi [i] + sfi [nbr_coef - i];

   768 

   769 			df2r [i] = vr * c + vi * s;

   770 			df2i [i] = vi * c - vr * s;

   771 		}

   772 

   773 		coef_index += d_nbr_coef;

   774 	}

   775 	while (coef_index < _length);

   776 }

   777 

   778 

   779 

   780 template <class DT>

   781 void	FFTReal <DT>::compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const

   782 {

   783 	assert (df != 0);

   784 	assert (sf != 0);

   785 	assert (df != sf);

   786 	assert (pass > TRIGO_BD_LIMIT);

   787 	assert (pass < _nbr_bits);

   788 

   789 	const long		nbr_coef = 1 << pass;

   790 	const long		h_nbr_coef = nbr_coef >> 1;

   791 	const long		d_nbr_coef = nbr_coef << 1;

   792 	long				coef_index = 0;

   793 	OscType &		osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];

   794 	do

   795 	{

   796 		const DataType	* const	sfr = sf + coef_index;

   797 		const DataType	* const	sfi = sfr + nbr_coef;

   798 		DataType			* const	df1r = df + coef_index;

   799 		DataType			* const	df2r = df1r + nbr_coef;

   800 

   801 		osc.clear_buffers ();

   802 

   803 		// Extreme coefficients are always real

   804 		df1r [0] = sfr [0] + sfi [0];		// + sfr [nbr_coef]

   805 		df2r [0] = sfr [0] - sfi [0];		// - sfr [nbr_coef]

   806 		df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;

   807 		df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;

   808 

   809 		// Others are conjugate complex numbers

   810 		DataType * const	df1i = df1r + h_nbr_coef;

   811 		DataType * const	df2i = df1i + nbr_coef;

   812 		for (long i = 1; i < h_nbr_coef; ++ i)

   813 		{

   814 			df1r [i] = sfr [i] + sfi [-i];		// + sfr [nbr_coef - i]

   815 			df1i [i] = sfi [i] - sfi [nbr_coef - i];

   816 

   817 			osc.step ();

   818 			const DataType	c = osc.get_cos ();

   819 			const DataType	s = osc.get_sin ();

   820 			const DataType	vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]

   821 			const DataType	vi = sfi [i] + sfi [nbr_coef - i];

   822 

   823 			df2r [i] = vr * c + vi * s;

   824 			df2i [i] = vi * c - vr * s;

   825 		}

   826 

   827 		coef_index += d_nbr_coef;

   828 	}

   829 	while (coef_index < _length);

   830 }

   831 

   832 

   833 

   834 template <class DT>

   835 void	FFTReal <DT>::compute_inverse_pass_3 (DataType df [], const DataType sf []) const

   836 {

   837 	assert (df != 0);

   838 	assert (sf != 0);

   839 	assert (df != sf);

   840 

   841 	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);

   842 	long				coef_index = 0;

   843 	do

   844 	{

   845 		df [coef_index] = sf [coef_index] + sf [coef_index + 4];

   846 		df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];

   847 		df [coef_index + 2] = sf [coef_index + 2] * 2;

   848 		df [coef_index + 6] = sf [coef_index + 6] * 2;

   849 

   850 		df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];

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

   852 

   853 		const DataType	vr = sf [coef_index + 1] - sf [coef_index + 3];

   854 		const DataType	vi = sf [coef_index + 5] + sf [coef_index + 7];

   855 

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

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

   858 

   859 		coef_index += 8;

   860 	}

   861 	while (coef_index < _length);

   862 }

   863 

   864 

   865 

   866 template <class DT>

   867 void	FFTReal <DT>::compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const

   868 {

   869 	assert (x != 0);

   870 	assert (sf != 0);

   871 	assert (x != sf);

   872 

   873 	const long *	bit_rev_lut_ptr = get_br_ptr ();

   874 	const DataType *	sf2 = sf;

   875 	long				coef_index = 0;

   876 	do

   877 	{

   878 		{

   879 			const DataType	b_0 = sf2 [0] + sf2 [2];

   880 			const DataType	b_2 = sf2 [0] - sf2 [2];

   881 			const DataType	b_1 = sf2 [1] * 2;

   882 			const DataType	b_3 = sf2 [3] * 2;

   883 

   884 			x [bit_rev_lut_ptr [0]] = b_0 + b_1;

   885 			x [bit_rev_lut_ptr [1]] = b_0 - b_1;

   886 			x [bit_rev_lut_ptr [2]] = b_2 + b_3;

   887 			x [bit_rev_lut_ptr [3]] = b_2 - b_3;

   888 		}

   889 		{

   890 			const DataType	b_0 = sf2 [4] + sf2 [6];

   891 			const DataType	b_2 = sf2 [4] - sf2 [6];

   892 			const DataType	b_1 = sf2 [5] * 2;

   893 			const DataType	b_3 = sf2 [7] * 2;

   894 

   895 			x [bit_rev_lut_ptr [4]] = b_0 + b_1;

   896 			x [bit_rev_lut_ptr [5]] = b_0 - b_1;

   897 			x [bit_rev_lut_ptr [6]] = b_2 + b_3;

   898 			x [bit_rev_lut_ptr [7]] = b_2 - b_3;

   899 		}

   900 

   901 		sf2 += 8;

   902 		coef_index += 8;

   903 		bit_rev_lut_ptr += 8;

   904 	}

   905 	while (coef_index < _length);

   906 }

   907 

   908 

   909 

   910 #endif	// FFTReal_CODEHEADER_INCLUDED

   911 

   912 #undef FFTReal_CURRENT_CODEHEADER

   913 

   914 

   915 

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