source: trunk/demos/spectrum/3rdparty/fftreal/FFTRealFixLen.hpp@ 822

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

        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 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/