BogaudioModules

BogaudioModules for VCV Rack
Log | Files | Refs | README | LICENSE

FFTReal.hpp (21034B)


      1 /*****************************************************************************
      2 
      3         FFTReal.hpp
      4         By Laurent de Soras
      5 
      6 --- Legal stuff ---
      7 
      8 This program is free software. It comes without any warranty, to
      9 the extent permitted by applicable law. You can redistribute it
     10 and/or modify it under the terms of the Do What The Fuck You Want
     11 To Public License, Version 2, as published by Sam Hocevar. See
     12 http://sam.zoy.org/wtfpl/COPYING for more details.
     13 
     14 *Tab=3***********************************************************************/
     15 
     16 
     17 
     18 #if defined (ffft_FFTReal_CURRENT_CODEHEADER)
     19 	#error Recursive inclusion of FFTReal code header.
     20 #endif
     21 #define	ffft_FFTReal_CURRENT_CODEHEADER
     22 
     23 #if ! defined (ffft_FFTReal_CODEHEADER_INCLUDED)
     24 #define	ffft_FFTReal_CODEHEADER_INCLUDED
     25 
     26 
     27 
     28 /*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
     29 
     30 #include	<cassert>
     31 #include	<cmath>
     32 
     33 
     34 
     35 namespace ffft
     36 {
     37 
     38 
     39 
     40 static inline bool	FFTReal_is_pow2 (long x)
     41 {
     42 	assert (x > 0);
     43 
     44 	return  ((x & -x) == x);
     45 }
     46 
     47 
     48 
     49 static inline int	FFTReal_get_next_pow2 (long x)
     50 {
     51 	--x;
     52 
     53 	int				p = 0;
     54 	while ((x & ~0xFFFFL) != 0)
     55 	{
     56 		p += 16;
     57 		x >>= 16;
     58 	}
     59 	while ((x & ~0xFL) != 0)
     60 	{
     61 		p += 4;
     62 		x >>= 4;
     63 	}
     64 	while (x > 0)
     65 	{
     66 		++p;
     67 		x >>= 1;
     68 	}
     69 
     70 	return (p);
     71 }
     72 
     73 
     74 
     75 /*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
     76 
     77 
     78 
     79 /*
     80 ==============================================================================
     81 Name: ctor
     82 Input parameters:
     83 	- length: length of the array on which we want to do a FFT. Range: power of
     84 		2 only, > 0.
     85 Throws: std::bad_alloc
     86 ==============================================================================
     87 */
     88 
     89 template <class DT>
     90 FFTReal <DT>::FFTReal (long length)
     91 :	_length (length)
     92 ,	_nbr_bits (FFTReal_get_next_pow2 (length))
     93 ,	_br_lut ()
     94 ,	_trigo_lut ()
     95 ,	_buffer (length)
     96 ,	_trigo_osc ()
     97 {
     98 	assert (FFTReal_is_pow2 (length));
     99 	assert (_nbr_bits <= MAX_BIT_DEPTH);
    100 
    101 	init_br_lut ();
    102 	init_trigo_lut ();
    103 	init_trigo_osc ();
    104 }
    105 
    106 
    107 
    108 /*
    109 ==============================================================================
    110 Name: get_length
    111 Description:
    112 	Returns the number of points processed by this FFT object.
    113 Returns: The number of points, power of 2, > 0.
    114 Throws: Nothing
    115 ==============================================================================
    116 */
    117 
    118 template <class DT>
    119 long	FFTReal <DT>::get_length () const
    120 {
    121 	return (_length);
    122 }
    123 
    124 
    125 
    126 /*
    127 ==============================================================================
    128 Name: do_fft
    129 Description:
    130 	Compute the FFT of the array.
    131 Input parameters:
    132 	- x: pointer on the source array (time).
    133 Output parameters:
    134 	- f: pointer on the destination array (frequencies).
    135 		f [0...length(x)/2] = real values,
    136 		f [length(x)/2+1...length(x)-1] = negative imaginary values of
    137 		coefficents 1...length(x)/2-1.
    138 Throws: Nothing
    139 ==============================================================================
    140 */
    141 
    142 template <class DT>
    143 void	FFTReal <DT>::do_fft (DataType f [], const DataType x []) const
    144 {
    145 	assert (f != 0);
    146 	assert (f != use_buffer ());
    147 	assert (x != 0);
    148 	assert (x != use_buffer ());
    149 	assert (x != f);
    150 
    151 	// General case
    152 	if (_nbr_bits > 2)
    153 	{
    154 		compute_fft_general (f, x);
    155 	}
    156 
    157 	// 4-point FFT
    158 	else if (_nbr_bits == 2)
    159 	{
    160 		f [1] = x [0] - x [2];
    161 		f [3] = x [1] - x [3];
    162 
    163 		const DataType	b_0 = x [0] + x [2];
    164 		const DataType	b_2 = x [1] + x [3];
    165 		
    166 		f [0] = b_0 + b_2;
    167 		f [2] = b_0 - b_2;
    168 	}
    169 
    170 	// 2-point FFT
    171 	else if (_nbr_bits == 1)
    172 	{
    173 		f [0] = x [0] + x [1];
    174 		f [1] = x [0] - x [1];
    175 	}
    176 
    177 	// 1-point FFT
    178 	else
    179 	{
    180 		f [0] = x [0];
    181 	}
    182 }
    183 
    184 
    185 
    186 /*
    187 ==============================================================================
    188 Name: do_ifft
    189 Description:
    190 	Compute the inverse FFT of the array. Note that data must be post-scaled:
    191 	IFFT (FFT (x)) = x * length (x).
    192 Input parameters:
    193 	- f: pointer on the source array (frequencies).
    194 		f [0...length(x)/2] = real values
    195 		f [length(x)/2+1...length(x)-1] = negative imaginary values of
    196 		coefficents 1...length(x)/2-1.
    197 Output parameters:
    198 	- x: pointer on the destination array (time).
    199 Throws: Nothing
    200 ==============================================================================
    201 */
    202 
    203 template <class DT>
    204 void	FFTReal <DT>::do_ifft (const DataType f [], DataType x []) const
    205 {
    206 	assert (f != 0);
    207 	assert (f != use_buffer ());
    208 	assert (x != 0);
    209 	assert (x != use_buffer ());
    210 	assert (x != f);
    211 
    212 	// General case
    213 	if (_nbr_bits > 2)
    214 	{
    215 		compute_ifft_general (f, x);
    216 	}
    217 
    218 	// 4-point IFFT
    219 	else if (_nbr_bits == 2)
    220 	{
    221 		const DataType	b_0 = f [0] + f [2];
    222 		const DataType	b_2 = f [0] - f [2];
    223 
    224 		x [0] = b_0 + f [1] * 2;
    225 		x [2] = b_0 - f [1] * 2;
    226 		x [1] = b_2 + f [3] * 2;
    227 		x [3] = b_2 - f [3] * 2;
    228 	}
    229 
    230 	// 2-point IFFT
    231 	else if (_nbr_bits == 1)
    232 	{
    233 		x [0] = f [0] + f [1];
    234 		x [1] = f [0] - f [1];
    235 	}
    236 
    237 	// 1-point IFFT
    238 	else
    239 	{
    240 		x [0] = f [0];
    241 	}
    242 }
    243 
    244 
    245 
    246 /*
    247 ==============================================================================
    248 Name: rescale
    249 Description:
    250 	Scale an array by divide each element by its length. This function should
    251 	be called after FFT + IFFT.
    252 Input parameters:
    253 	- x: pointer on array to rescale (time or frequency).
    254 Throws: Nothing
    255 ==============================================================================
    256 */
    257 
    258 template <class DT>
    259 void	FFTReal <DT>::rescale (DataType x []) const
    260 {
    261 	const DataType	mul = DataType (1.0 / _length);
    262 
    263 	if (_length < 4)
    264 	{
    265 		long				i = _length - 1;
    266 		do
    267 		{
    268 			x [i] *= mul;
    269 			--i;
    270 		}
    271 		while (i >= 0);
    272 	}
    273 
    274 	else
    275 	{
    276 		assert ((_length & 3) == 0);
    277 
    278 		// Could be optimized with SIMD instruction sets (needs alignment check)
    279 		long				i = _length - 4;
    280 		do
    281 		{
    282 			x [i + 0] *= mul;
    283 			x [i + 1] *= mul;
    284 			x [i + 2] *= mul;
    285 			x [i + 3] *= mul;
    286 			i -= 4;
    287 		}
    288 		while (i >= 0);
    289 	}
    290 }
    291 
    292 
    293 
    294 /*
    295 ==============================================================================
    296 Name: use_buffer
    297 Description:
    298 	Access the internal buffer, whose length is the FFT one.
    299 	Buffer content will be erased at each do_fft() / do_ifft() call!
    300 	This buffer cannot be used as:
    301 		- source for FFT or IFFT done with this object
    302 		- destination for FFT or IFFT done with this object
    303 Returns:
    304 	Buffer start address
    305 Throws: Nothing
    306 ==============================================================================
    307 */
    308 
    309 template <class DT>
    310 typename FFTReal <DT>::DataType *	FFTReal <DT>::use_buffer () const
    311 {
    312 	return (&_buffer [0]);
    313 }
    314 
    315 
    316 
    317 /*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
    318 
    319 
    320 
    321 /*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
    322 
    323 
    324 
    325 template <class DT>
    326 void	FFTReal <DT>::init_br_lut ()
    327 {
    328 	const long		length = 1L << _nbr_bits;
    329 	_br_lut.resize (length);
    330 
    331 	_br_lut [0] = 0;
    332 	long				br_index = 0;
    333 	for (long cnt = 1; cnt < length; ++cnt)
    334 	{
    335 		// ++br_index (bit reversed)
    336 		long				bit = length >> 1;
    337 		while (((br_index ^= bit) & bit) == 0)
    338 		{
    339 			bit >>= 1;
    340 		}
    341 
    342 		_br_lut [cnt] = br_index;
    343 	}
    344 }
    345 
    346 
    347 
    348 template <class DT>
    349 void	FFTReal <DT>::init_trigo_lut ()
    350 {
    351 	using namespace std;
    352 
    353 	if (_nbr_bits > 3)
    354 	{
    355 		const long		total_len = (1L << (_nbr_bits - 1)) - 4;
    356 		_trigo_lut.resize (total_len);
    357 
    358 		for (int level = 3; level < _nbr_bits; ++level)
    359 		{
    360 			const long		level_len = 1L << (level - 1);
    361 			DataType	* const	level_ptr =
    362 				&_trigo_lut [get_trigo_level_index (level)];
    363 			const double	mul = PI / (level_len << 1);
    364 
    365 			for (long i = 0; i < level_len; ++ i)
    366 			{
    367 				level_ptr [i] = static_cast <DataType> (cos (i * mul));
    368 			}
    369 		}
    370 	}
    371 }
    372 
    373 
    374 
    375 template <class DT>
    376 void	FFTReal <DT>::init_trigo_osc ()
    377 {
    378 	const int		nbr_osc = _nbr_bits - TRIGO_BD_LIMIT;
    379 	if (nbr_osc > 0)
    380 	{
    381 		_trigo_osc.resize (nbr_osc);
    382 
    383 		for (int osc_cnt = 0; osc_cnt < nbr_osc; ++osc_cnt)
    384 		{
    385 			OscType &		osc = _trigo_osc [osc_cnt];
    386 
    387 			const long		len = 1L << (TRIGO_BD_LIMIT + osc_cnt);
    388 			const double	mul = (0.5 * PI) / len;
    389 			osc.set_step (mul);
    390 		}
    391 	}
    392 }
    393 
    394 
    395 
    396 template <class DT>
    397 const long *	FFTReal <DT>::get_br_ptr () const
    398 {
    399 	return (&_br_lut [0]);
    400 }
    401 
    402 
    403 
    404 template <class DT>
    405 const typename FFTReal <DT>::DataType *	FFTReal <DT>::get_trigo_ptr (int level) const
    406 {
    407 	assert (level >= 3);
    408 
    409 	return (&_trigo_lut [get_trigo_level_index (level)]);
    410 }
    411 
    412 
    413 
    414 template <class DT>
    415 long	FFTReal <DT>::get_trigo_level_index (int level) const
    416 {
    417 	assert (level >= 3);
    418 
    419 	return ((1L << (level - 1)) - 4);
    420 }
    421 
    422 
    423 
    424 // Transform in several passes
    425 template <class DT>
    426 void	FFTReal <DT>::compute_fft_general (DataType f [], const DataType x []) const
    427 {
    428 	assert (f != 0);
    429 	assert (f != use_buffer ());
    430 	assert (x != 0);
    431 	assert (x != use_buffer ());
    432 	assert (x != f);
    433 
    434 	DataType *		sf;
    435 	DataType *		df;
    436 
    437 	if ((_nbr_bits & 1) != 0)
    438 	{
    439 		df = use_buffer ();
    440 		sf = f;
    441 	}
    442 	else
    443 	{
    444 		df = f;
    445 		sf = use_buffer ();
    446 	}
    447 
    448 	compute_direct_pass_1_2 (df, x);
    449 	compute_direct_pass_3 (sf, df);
    450 
    451 	for (int pass = 3; pass < _nbr_bits; ++ pass)
    452 	{
    453 		compute_direct_pass_n (df, sf, pass);
    454 
    455 		DataType * const	temp_ptr = df;
    456 		df = sf;
    457 		sf = temp_ptr;
    458 	}
    459 }
    460 
    461 
    462 
    463 template <class DT>
    464 void	FFTReal <DT>::compute_direct_pass_1_2 (DataType df [], const DataType x []) const
    465 {
    466 	assert (df != 0);
    467 	assert (x != 0);
    468 	assert (df != x);
    469 
    470 	const long * const	bit_rev_lut_ptr = get_br_ptr ();
    471 	long				coef_index = 0;
    472 	do
    473 	{
    474 		const long		rev_index_0 = bit_rev_lut_ptr [coef_index];
    475 		const long		rev_index_1 = bit_rev_lut_ptr [coef_index + 1];
    476 		const long		rev_index_2 = bit_rev_lut_ptr [coef_index + 2];
    477 		const long		rev_index_3 = bit_rev_lut_ptr [coef_index + 3];
    478 
    479 		DataType	* const	df2 = df + coef_index;
    480 		df2 [1] = x [rev_index_0] - x [rev_index_1];
    481 		df2 [3] = x [rev_index_2] - x [rev_index_3];
    482 
    483 		const DataType	sf_0 = x [rev_index_0] + x [rev_index_1];
    484 		const DataType	sf_2 = x [rev_index_2] + x [rev_index_3];
    485 
    486 		df2 [0] = sf_0 + sf_2;
    487 		df2 [2] = sf_0 - sf_2;
    488 		
    489 		coef_index += 4;
    490 	}
    491 	while (coef_index < _length);
    492 }
    493 
    494 
    495 
    496 template <class DT>
    497 void	FFTReal <DT>::compute_direct_pass_3 (DataType df [], const DataType sf []) const
    498 {
    499 	assert (df != 0);
    500 	assert (sf != 0);
    501 	assert (df != sf);
    502 
    503 	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);
    504 	long				coef_index = 0;
    505 	do
    506 	{
    507 		DataType			v;
    508 
    509 		df [coef_index] = sf [coef_index] + sf [coef_index + 4];
    510 		df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
    511 		df [coef_index + 2] = sf [coef_index + 2];
    512 		df [coef_index + 6] = sf [coef_index + 6];
    513 
    514 		v = (sf [coef_index + 5] - sf [coef_index + 7]) * sqrt2_2;
    515 		df [coef_index + 1] = sf [coef_index + 1] + v;
    516 		df [coef_index + 3] = sf [coef_index + 1] - v;
    517 
    518 		v = (sf [coef_index + 5] + sf [coef_index + 7]) * sqrt2_2;
    519 		df [coef_index + 5] = v + sf [coef_index + 3];
    520 		df [coef_index + 7] = v - sf [coef_index + 3];
    521 
    522 		coef_index += 8;
    523 	}
    524 	while (coef_index < _length);
    525 }
    526 
    527 
    528 
    529 template <class DT>
    530 void	FFTReal <DT>::compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const
    531 {
    532 	assert (df != 0);
    533 	assert (sf != 0);
    534 	assert (df != sf);
    535 	assert (pass >= 3);
    536 	assert (pass < _nbr_bits);
    537 
    538 	if (pass <= TRIGO_BD_LIMIT)
    539 	{
    540 		compute_direct_pass_n_lut (df, sf, pass);
    541 	}
    542 	else
    543 	{
    544 		compute_direct_pass_n_osc (df, sf, pass);
    545 	}
    546 }
    547 
    548 
    549 
    550 template <class DT>
    551 void	FFTReal <DT>::compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const
    552 {
    553 	assert (df != 0);
    554 	assert (sf != 0);
    555 	assert (df != sf);
    556 	assert (pass >= 3);
    557 	assert (pass < _nbr_bits);
    558 
    559 	const long		nbr_coef = 1 << pass;
    560 	const long		h_nbr_coef = nbr_coef >> 1;
    561 	const long		d_nbr_coef = nbr_coef << 1;
    562 	long				coef_index = 0;
    563 	const DataType	* const	cos_ptr = get_trigo_ptr (pass);
    564 	do
    565 	{
    566 		const DataType	* const	sf1r = sf + coef_index;
    567 		const DataType	* const	sf2r = sf1r + nbr_coef;
    568 		DataType			* const	dfr = df + coef_index;
    569 		DataType			* const	dfi = dfr + nbr_coef;
    570 
    571 		// Extreme coefficients are always real
    572 		dfr [0] = sf1r [0] + sf2r [0];
    573 		dfi [0] = sf1r [0] - sf2r [0];	// dfr [nbr_coef] =
    574 		dfr [h_nbr_coef] = sf1r [h_nbr_coef];
    575 		dfi [h_nbr_coef] = sf2r [h_nbr_coef];
    576 
    577 		// Others are conjugate complex numbers
    578 		const DataType * const	sf1i = sf1r + h_nbr_coef;
    579 		const DataType * const	sf2i = sf1i + nbr_coef;
    580 		for (long i = 1; i < h_nbr_coef; ++ i)
    581 		{
    582 			const DataType	c = cos_ptr [i];					// cos (i*PI/nbr_coef);
    583 			const DataType	s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);
    584 			DataType	 		v;
    585 
    586 			v = sf2r [i] * c - sf2i [i] * s;
    587 			dfr [i] = sf1r [i] + v;
    588 			dfi [-i] = sf1r [i] - v;	// dfr [nbr_coef - i] =
    589 
    590 			v = sf2r [i] * s + sf2i [i] * c;
    591 			dfi [i] = v + sf1i [i];
    592 			dfi [nbr_coef - i] = v - sf1i [i];
    593 		}
    594 
    595 		coef_index += d_nbr_coef;
    596 	}
    597 	while (coef_index < _length);
    598 }
    599 
    600 
    601 
    602 template <class DT>
    603 void	FFTReal <DT>::compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const
    604 {
    605 	assert (df != 0);
    606 	assert (sf != 0);
    607 	assert (df != sf);
    608 	assert (pass > TRIGO_BD_LIMIT);
    609 	assert (pass < _nbr_bits);
    610 
    611 	const long		nbr_coef = 1 << pass;
    612 	const long		h_nbr_coef = nbr_coef >> 1;
    613 	const long		d_nbr_coef = nbr_coef << 1;
    614 	long				coef_index = 0;
    615 	OscType &		osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
    616 	do
    617 	{
    618 		const DataType	* const	sf1r = sf + coef_index;
    619 		const DataType	* const	sf2r = sf1r + nbr_coef;
    620 		DataType			* const	dfr = df + coef_index;
    621 		DataType			* const	dfi = dfr + nbr_coef;
    622 
    623 		osc.clear_buffers ();
    624 
    625 		// Extreme coefficients are always real
    626 		dfr [0] = sf1r [0] + sf2r [0];
    627 		dfi [0] = sf1r [0] - sf2r [0];	// dfr [nbr_coef] =
    628 		dfr [h_nbr_coef] = sf1r [h_nbr_coef];
    629 		dfi [h_nbr_coef] = sf2r [h_nbr_coef];
    630 
    631 		// Others are conjugate complex numbers
    632 		const DataType * const	sf1i = sf1r + h_nbr_coef;
    633 		const DataType * const	sf2i = sf1i + nbr_coef;
    634 		for (long i = 1; i < h_nbr_coef; ++ i)
    635 		{
    636 			osc.step ();
    637 			const DataType	c = osc.get_cos ();
    638 			const DataType	s = osc.get_sin ();
    639 			DataType	 		v;
    640 
    641 			v = sf2r [i] * c - sf2i [i] * s;
    642 			dfr [i] = sf1r [i] + v;
    643 			dfi [-i] = sf1r [i] - v;	// dfr [nbr_coef - i] =
    644 
    645 			v = sf2r [i] * s + sf2i [i] * c;
    646 			dfi [i] = v + sf1i [i];
    647 			dfi [nbr_coef - i] = v - sf1i [i];
    648 		}
    649 
    650 		coef_index += d_nbr_coef;
    651 	}
    652 	while (coef_index < _length);
    653 }
    654 
    655 
    656 
    657 // Transform in several pass
    658 template <class DT>
    659 void	FFTReal <DT>::compute_ifft_general (const DataType f [], DataType x []) const
    660 {
    661 	assert (f != 0);
    662 	assert (f != use_buffer ());
    663 	assert (x != 0);
    664 	assert (x != use_buffer ());
    665 	assert (x != f);
    666 
    667 	DataType *		sf = const_cast <DataType *> (f);
    668 	DataType *		df;
    669 	DataType *		df_temp;
    670 
    671 	if (_nbr_bits & 1)
    672 	{
    673 		df = use_buffer ();
    674 		df_temp = x;
    675 	}
    676 	else
    677 	{
    678 		df = x;
    679 		df_temp = use_buffer ();
    680 	}
    681 
    682 	for (int pass = _nbr_bits - 1; pass >= 3; -- pass)
    683 	{
    684 		compute_inverse_pass_n (df, sf, pass);
    685 
    686 		if (pass < _nbr_bits - 1)
    687 		{
    688 			DataType	* const	temp_ptr = df;
    689 			df = sf;
    690 			sf = temp_ptr;
    691 		}
    692 		else
    693 		{
    694 			sf = df;
    695 			df = df_temp;
    696 		}
    697 	}
    698 
    699 	compute_inverse_pass_3 (df, sf);
    700 	compute_inverse_pass_1_2 (x, df);
    701 }
    702 
    703 
    704 
    705 template <class DT>
    706 void	FFTReal <DT>::compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const
    707 {
    708 	assert (df != 0);
    709 	assert (sf != 0);
    710 	assert (df != sf);
    711 	assert (pass >= 3);
    712 	assert (pass < _nbr_bits);
    713 
    714 	if (pass <= TRIGO_BD_LIMIT)
    715 	{
    716 		compute_inverse_pass_n_lut (df, sf, pass);
    717 	}
    718 	else
    719 	{
    720 		compute_inverse_pass_n_osc (df, sf, pass);
    721 	}
    722 }
    723 
    724 
    725 
    726 template <class DT>
    727 void	FFTReal <DT>::compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const
    728 {
    729 	assert (df != 0);
    730 	assert (sf != 0);
    731 	assert (df != sf);
    732 	assert (pass >= 3);
    733 	assert (pass < _nbr_bits);
    734 
    735 	const long		nbr_coef = 1 << pass;
    736 	const long		h_nbr_coef = nbr_coef >> 1;
    737 	const long		d_nbr_coef = nbr_coef << 1;
    738 	long				coef_index = 0;
    739 	const DataType * const	cos_ptr = get_trigo_ptr (pass);
    740 	do
    741 	{
    742 		const DataType	* const	sfr = sf + coef_index;
    743 		const DataType	* const	sfi = sfr + nbr_coef;
    744 		DataType			* const	df1r = df + coef_index;
    745 		DataType			* const	df2r = df1r + nbr_coef;
    746 
    747 		// Extreme coefficients are always real
    748 		df1r [0] = sfr [0] + sfi [0];		// + sfr [nbr_coef]
    749 		df2r [0] = sfr [0] - sfi [0];		// - sfr [nbr_coef]
    750 		df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
    751 		df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
    752 
    753 		// Others are conjugate complex numbers
    754 		DataType * const	df1i = df1r + h_nbr_coef;
    755 		DataType * const	df2i = df1i + nbr_coef;
    756 		for (long i = 1; i < h_nbr_coef; ++ i)
    757 		{
    758 			df1r [i] = sfr [i] + sfi [-i];		// + sfr [nbr_coef - i]
    759 			df1i [i] = sfi [i] - sfi [nbr_coef - i];
    760 
    761 			const DataType	c = cos_ptr [i];					// cos (i*PI/nbr_coef);
    762 			const DataType	s = cos_ptr [h_nbr_coef - i];	// sin (i*PI/nbr_coef);
    763 			const DataType	vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]
    764 			const DataType	vi = sfi [i] + sfi [nbr_coef - i];
    765 
    766 			df2r [i] = vr * c + vi * s;
    767 			df2i [i] = vi * c - vr * s;
    768 		}
    769 
    770 		coef_index += d_nbr_coef;
    771 	}
    772 	while (coef_index < _length);
    773 }
    774 
    775 
    776 
    777 template <class DT>
    778 void	FFTReal <DT>::compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const
    779 {
    780 	assert (df != 0);
    781 	assert (sf != 0);
    782 	assert (df != sf);
    783 	assert (pass > TRIGO_BD_LIMIT);
    784 	assert (pass < _nbr_bits);
    785 
    786 	const long		nbr_coef = 1 << pass;
    787 	const long		h_nbr_coef = nbr_coef >> 1;
    788 	const long		d_nbr_coef = nbr_coef << 1;
    789 	long				coef_index = 0;
    790 	OscType &		osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
    791 	do
    792 	{
    793 		const DataType	* const	sfr = sf + coef_index;
    794 		const DataType	* const	sfi = sfr + nbr_coef;
    795 		DataType			* const	df1r = df + coef_index;
    796 		DataType			* const	df2r = df1r + nbr_coef;
    797 
    798 		osc.clear_buffers ();
    799 
    800 		// Extreme coefficients are always real
    801 		df1r [0] = sfr [0] + sfi [0];		// + sfr [nbr_coef]
    802 		df2r [0] = sfr [0] - sfi [0];		// - sfr [nbr_coef]
    803 		df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
    804 		df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
    805 
    806 		// Others are conjugate complex numbers
    807 		DataType * const	df1i = df1r + h_nbr_coef;
    808 		DataType * const	df2i = df1i + nbr_coef;
    809 		for (long i = 1; i < h_nbr_coef; ++ i)
    810 		{
    811 			df1r [i] = sfr [i] + sfi [-i];		// + sfr [nbr_coef - i]
    812 			df1i [i] = sfi [i] - sfi [nbr_coef - i];
    813 
    814 			osc.step ();
    815 			const DataType	c = osc.get_cos ();
    816 			const DataType	s = osc.get_sin ();
    817 			const DataType	vr = sfr [i] - sfi [-i];		// - sfr [nbr_coef - i]
    818 			const DataType	vi = sfi [i] + sfi [nbr_coef - i];
    819 
    820 			df2r [i] = vr * c + vi * s;
    821 			df2i [i] = vi * c - vr * s;
    822 		}
    823 
    824 		coef_index += d_nbr_coef;
    825 	}
    826 	while (coef_index < _length);
    827 }
    828 
    829 
    830 
    831 template <class DT>
    832 void	FFTReal <DT>::compute_inverse_pass_3 (DataType df [], const DataType sf []) const
    833 {
    834 	assert (df != 0);
    835 	assert (sf != 0);
    836 	assert (df != sf);
    837 
    838 	const DataType	sqrt2_2 = DataType (SQRT2 * 0.5);
    839 	long				coef_index = 0;
    840 	do
    841 	{
    842 		df [coef_index] = sf [coef_index] + sf [coef_index + 4];
    843 		df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
    844 		df [coef_index + 2] = sf [coef_index + 2] * 2;
    845 		df [coef_index + 6] = sf [coef_index + 6] * 2;
    846 
    847 		df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];
    848 		df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7];
    849 
    850 		const DataType	vr = sf [coef_index + 1] - sf [coef_index + 3];
    851 		const DataType	vi = sf [coef_index + 5] + sf [coef_index + 7];
    852 
    853 		df [coef_index + 5] = (vr + vi) * sqrt2_2;
    854 		df [coef_index + 7] = (vi - vr) * sqrt2_2;
    855 
    856 		coef_index += 8;
    857 	}
    858 	while (coef_index < _length);
    859 }
    860 
    861 
    862 
    863 template <class DT>
    864 void	FFTReal <DT>::compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const
    865 {
    866 	assert (x != 0);
    867 	assert (sf != 0);
    868 	assert (x != sf);
    869 
    870 	const long *	bit_rev_lut_ptr = get_br_ptr ();
    871 	const DataType *	sf2 = sf;
    872 	long				coef_index = 0;
    873 	do
    874 	{
    875 		{
    876 			const DataType	b_0 = sf2 [0] + sf2 [2];
    877 			const DataType	b_2 = sf2 [0] - sf2 [2];
    878 			const DataType	b_1 = sf2 [1] * 2;
    879 			const DataType	b_3 = sf2 [3] * 2;
    880 
    881 			x [bit_rev_lut_ptr [0]] = b_0 + b_1;
    882 			x [bit_rev_lut_ptr [1]] = b_0 - b_1;
    883 			x [bit_rev_lut_ptr [2]] = b_2 + b_3;
    884 			x [bit_rev_lut_ptr [3]] = b_2 - b_3;
    885 		}
    886 		{
    887 			const DataType	b_0 = sf2 [4] + sf2 [6];
    888 			const DataType	b_2 = sf2 [4] - sf2 [6];
    889 			const DataType	b_1 = sf2 [5] * 2;
    890 			const DataType	b_3 = sf2 [7] * 2;
    891 
    892 			x [bit_rev_lut_ptr [4]] = b_0 + b_1;
    893 			x [bit_rev_lut_ptr [5]] = b_0 - b_1;
    894 			x [bit_rev_lut_ptr [6]] = b_2 + b_3;
    895 			x [bit_rev_lut_ptr [7]] = b_2 - b_3;
    896 		}
    897 
    898 		sf2 += 8;
    899 		coef_index += 8;
    900 		bit_rev_lut_ptr += 8;
    901 	}
    902 	while (coef_index < _length);
    903 }
    904 
    905 
    906 
    907 }	// namespace ffft
    908 
    909 
    910 
    911 #endif	// ffft_FFTReal_CODEHEADER_INCLUDED
    912 
    913 #undef ffft_FFTReal_CURRENT_CODEHEADER
    914 
    915 
    916 
    917 /*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/