Serial-Studio/lib/QRealFourier/fftreal/FFTRealFixLen.hpp

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

        FFTRealFixLen.hpp
        By Laurent de Soras

--- Legal stuff ---

This program is free software. It comes without any warranty, to
the extent permitted by applicable law. You can redistribute it
and/or modify it under the terms of the Do What The Fuck You Want
To Public License, Version 2, as published by Sam Hocevar. See
http://sam.zoy.org/wtfpl/COPYING for more details.

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

#if defined(ffft_FFTRealFixLen_CURRENT_CODEHEADER)
#  error Recursive inclusion of FFTRealFixLen code header.
#endif
#define ffft_FFTRealFixLen_CURRENT_CODEHEADER

#if !defined(ffft_FFTRealFixLen_CODEHEADER_INCLUDED)
#  define ffft_FFTRealFixLen_CODEHEADER_INCLUDED

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

#  include "def.h"
#  include "FFTRealPassDirect.h"
#  include "FFTRealPassInverse.h"
#  include "FFTRealSelect.h"

#  include <cassert>
#  include <cmath>

namespace std
{
}

namespace ffft
{

/*\\\ 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<>
inline 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<>
inline 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<>
inline 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<>
inline 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<>
inline 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<>
inline 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);
  }
}

} // namespace ffft

#endif // ffft_FFTRealFixLen_CODEHEADER_INCLUDED

#undef ffft_FFTRealFixLen_CURRENT_CODEHEADER

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