diff options
Diffstat (limited to 'examples/spectrum/3rdparty/fftreal/FFTRealFixLen.hpp')
| -rw-r--r-- | examples/spectrum/3rdparty/fftreal/FFTRealFixLen.hpp | 322 |
1 files changed, 322 insertions, 0 deletions
diff --git a/examples/spectrum/3rdparty/fftreal/FFTRealFixLen.hpp b/examples/spectrum/3rdparty/fftreal/FFTRealFixLen.hpp new file mode 100644 index 00000000..6defb009 --- /dev/null +++ b/examples/spectrum/3rdparty/fftreal/FFTRealFixLen.hpp @@ -0,0 +1,322 @@ +/***************************************************************************** + + FFTRealFixLen.hpp + Copyright (c) 2005 Laurent de Soras + +--- Legal stuff --- + +This library is free software; you can redistribute it and/or +modify it under the terms of the GNU Lesser General Public +License as published by the Free Software Foundation; either +version 2.1 of the License, or (at your option) any later version. + +This library is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public +License along with this library; if not, write to the Free Software +Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +*Tab=3***********************************************************************/ + + + +#if defined (FFTRealFixLen_CURRENT_CODEHEADER) + #error Recursive inclusion of FFTRealFixLen code header. +#endif +#define FFTRealFixLen_CURRENT_CODEHEADER + +#if ! defined (FFTRealFixLen_CODEHEADER_INCLUDED) +#define FFTRealFixLen_CODEHEADER_INCLUDED + + + +/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + +#include "def.h" +#include "FFTRealPassDirect.h" +#include "FFTRealPassInverse.h" +#include "FFTRealSelect.h" + +#include <cassert> +#include <cmath> + +namespace std { } + + + +/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +template <int LL2> +FFTRealFixLen <LL2>::FFTRealFixLen () +: _buffer (FFT_LEN) +, _br_data (BR_ARR_SIZE) +, _trigo_data (TRIGO_TABLE_ARR_SIZE) +, _trigo_osc () +{ + build_br_lut (); + build_trigo_lut (); + build_trigo_osc (); +} + + + +template <int LL2> +long FFTRealFixLen <LL2>::get_length () const +{ + return (FFT_LEN); +} + + + +// General case +template <int LL2> +void FFTRealFixLen <LL2>::do_fft (DataType f [], const DataType x []) +{ + assert (f != 0); + assert (x != 0); + assert (x != f); + assert (FFT_LEN_L2 >= 3); + + // Do the transform in several passes + const DataType * cos_ptr = &_trigo_data [0]; + const long * br_ptr = &_br_data [0]; + + FFTRealPassDirect <FFT_LEN_L2 - 1>::process ( + FFT_LEN, + f, + &_buffer [0], + x, + cos_ptr, + TRIGO_TABLE_ARR_SIZE, + br_ptr, + &_trigo_osc [0] + ); +} + +// 4-point FFT +template <> +void FFTRealFixLen <2>::do_fft (DataType f [], const DataType x []) +{ + assert (f != 0); + assert (x != 0); + assert (x != f); + + f [1] = x [0] - x [2]; + f [3] = x [1] - x [3]; + + const DataType b_0 = x [0] + x [2]; + const DataType b_2 = x [1] + x [3]; + + f [0] = b_0 + b_2; + f [2] = b_0 - b_2; +} + +// 2-point FFT +template <> +void FFTRealFixLen <1>::do_fft (DataType f [], const DataType x []) +{ + assert (f != 0); + assert (x != 0); + assert (x != f); + + f [0] = x [0] + x [1]; + f [1] = x [0] - x [1]; +} + +// 1-point FFT +template <> +void FFTRealFixLen <0>::do_fft (DataType f [], const DataType x []) +{ + assert (f != 0); + assert (x != 0); + + f [0] = x [0]; +} + + + +// General case +template <int LL2> +void FFTRealFixLen <LL2>::do_ifft (const DataType f [], DataType x []) +{ + assert (f != 0); + assert (x != 0); + assert (x != f); + assert (FFT_LEN_L2 >= 3); + + // Do the transform in several passes + DataType * s_ptr = + FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (&_buffer [0], x); + DataType * d_ptr = + FFTRealSelect <FFT_LEN_L2 & 1>::sel_bin (x, &_buffer [0]); + const DataType * cos_ptr = &_trigo_data [0]; + const long * br_ptr = &_br_data [0]; + + FFTRealPassInverse <FFT_LEN_L2 - 1>::process ( + FFT_LEN, + d_ptr, + s_ptr, + f, + cos_ptr, + TRIGO_TABLE_ARR_SIZE, + br_ptr, + &_trigo_osc [0] + ); +} + +// 4-point IFFT +template <> +void FFTRealFixLen <2>::do_ifft (const DataType f [], DataType x []) +{ + assert (f != 0); + assert (x != 0); + assert (x != f); + + const DataType b_0 = f [0] + f [2]; + const DataType 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 +template <> +void FFTRealFixLen <1>::do_ifft (const DataType f [], DataType x []) +{ + assert (f != 0); + assert (x != 0); + assert (x != f); + + x [0] = f [0] + f [1]; + x [1] = f [0] - f [1]; +} + +// 1-point IFFT +template <> +void FFTRealFixLen <0>::do_ifft (const DataType f [], DataType x []) +{ + assert (f != 0); + assert (x != 0); + assert (x != f); + + x [0] = f [0]; +} + + + + +template <int LL2> +void FFTRealFixLen <LL2>::rescale (DataType x []) const +{ + assert (x != 0); + + const DataType mul = DataType (1.0 / FFT_LEN); + + if (FFT_LEN < 4) + { + long i = FFT_LEN - 1; + do + { + x [i] *= mul; + --i; + } + while (i >= 0); + } + + else + { + assert ((FFT_LEN & 3) == 0); + + // Could be optimized with SIMD instruction sets (needs alignment check) + long i = FFT_LEN - 4; + do + { + x [i + 0] *= mul; + x [i + 1] *= mul; + x [i + 2] *= mul; + x [i + 3] *= mul; + i -= 4; + } + while (i >= 0); + } +} + + + +/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +template <int LL2> +void FFTRealFixLen <LL2>::build_br_lut () +{ + _br_data [0] = 0; + for (long cnt = 1; cnt < BR_ARR_SIZE; ++cnt) + { + long index = cnt << 2; + long br_index = 0; + + int bit_cnt = FFT_LEN_L2; + do + { + br_index <<= 1; + br_index += (index & 1); + index >>= 1; + + -- bit_cnt; + } + while (bit_cnt > 0); + + _br_data [cnt] = br_index; + } +} + + + +template <int LL2> +void FFTRealFixLen <LL2>::build_trigo_lut () +{ + const double mul = (0.5 * PI) / TRIGO_TABLE_ARR_SIZE; + for (long i = 0; i < TRIGO_TABLE_ARR_SIZE; ++ i) + { + using namespace std; + + _trigo_data [i] = DataType (cos (i * mul)); + } +} + + + +template <int LL2> +void FFTRealFixLen <LL2>::build_trigo_osc () +{ + for (int i = 0; i < NBR_TRIGO_OSC; ++i) + { + OscType & osc = _trigo_osc [i]; + + const long len = static_cast <long> (TRIGO_TABLE_ARR_SIZE) << (i + 1); + const double mul = (0.5 * PI) / len; + osc.set_step (mul); + } +} + + + +#endif // FFTRealFixLen_CODEHEADER_INCLUDED + +#undef FFTRealFixLen_CURRENT_CODEHEADER + + + +/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |