blob: 6defb0095bdbb5bad4e0c777c881d850ef0ed4cf [file] [log] [blame]
/*****************************************************************************
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 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/