60941c2741
This is the beginning of revision history for this module. If you want to look at revision history older than this, please refer to the Qt Git wiki for how to use Git history grafting. At the time of writing, this wiki is located here: http://qt.gitorious.org/qt/pages/GitIntroductionWithQt If you have already performed the grafting and you don't see any history beyond this commit, try running "git log" with the "--follow" argument. Branched from the monolithic repo, Qt master branch, at commit 896db169ea224deb96c59ce8af800d019de63f12
917 lines
20 KiB
C++
917 lines
20 KiB
C++
/*****************************************************************************
|
|
|
|
FFTReal.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 (FFTReal_CURRENT_CODEHEADER)
|
|
#error Recursive inclusion of FFTReal code header.
|
|
#endif
|
|
#define FFTReal_CURRENT_CODEHEADER
|
|
|
|
#if ! defined (FFTReal_CODEHEADER_INCLUDED)
|
|
#define FFTReal_CODEHEADER_INCLUDED
|
|
|
|
|
|
|
|
/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
|
|
|
|
#include <cassert>
|
|
#include <cmath>
|
|
|
|
|
|
|
|
static inline bool FFTReal_is_pow2 (long x)
|
|
{
|
|
assert (x > 0);
|
|
|
|
return ((x & -x) == x);
|
|
}
|
|
|
|
|
|
|
|
static inline int FFTReal_get_next_pow2 (long x)
|
|
{
|
|
--x;
|
|
|
|
int p = 0;
|
|
while ((x & ~0xFFFFL) != 0)
|
|
{
|
|
p += 16;
|
|
x >>= 16;
|
|
}
|
|
while ((x & ~0xFL) != 0)
|
|
{
|
|
p += 4;
|
|
x >>= 4;
|
|
}
|
|
while (x > 0)
|
|
{
|
|
++p;
|
|
x >>= 1;
|
|
}
|
|
|
|
return (p);
|
|
}
|
|
|
|
|
|
|
|
/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
|
|
|
|
|
|
|
|
/*
|
|
==============================================================================
|
|
Name: ctor
|
|
Input parameters:
|
|
- length: length of the array on which we want to do a FFT. Range: power of
|
|
2 only, > 0.
|
|
Throws: std::bad_alloc
|
|
==============================================================================
|
|
*/
|
|
|
|
template <class DT>
|
|
FFTReal <DT>::FFTReal (long length)
|
|
: _length (length)
|
|
, _nbr_bits (FFTReal_get_next_pow2 (length))
|
|
, _br_lut ()
|
|
, _trigo_lut ()
|
|
, _buffer (length)
|
|
, _trigo_osc ()
|
|
{
|
|
assert (FFTReal_is_pow2 (length));
|
|
assert (_nbr_bits <= MAX_BIT_DEPTH);
|
|
|
|
init_br_lut ();
|
|
init_trigo_lut ();
|
|
init_trigo_osc ();
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
==============================================================================
|
|
Name: get_length
|
|
Description:
|
|
Returns the number of points processed by this FFT object.
|
|
Returns: The number of points, power of 2, > 0.
|
|
Throws: Nothing
|
|
==============================================================================
|
|
*/
|
|
|
|
template <class DT>
|
|
long FFTReal <DT>::get_length () const
|
|
{
|
|
return (_length);
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
==============================================================================
|
|
Name: do_fft
|
|
Description:
|
|
Compute the FFT of the array.
|
|
Input parameters:
|
|
- x: pointer on the source array (time).
|
|
Output parameters:
|
|
- f: pointer on the destination array (frequencies).
|
|
f [0...length(x)/2] = real values,
|
|
f [length(x)/2+1...length(x)-1] = negative imaginary values of
|
|
coefficents 1...length(x)/2-1.
|
|
Throws: Nothing
|
|
==============================================================================
|
|
*/
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::do_fft (DataType f [], const DataType x []) const
|
|
{
|
|
assert (f != 0);
|
|
assert (f != use_buffer ());
|
|
assert (x != 0);
|
|
assert (x != use_buffer ());
|
|
assert (x != f);
|
|
|
|
// General case
|
|
if (_nbr_bits > 2)
|
|
{
|
|
compute_fft_general (f, x);
|
|
}
|
|
|
|
// 4-point FFT
|
|
else if (_nbr_bits == 2)
|
|
{
|
|
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
|
|
else if (_nbr_bits == 1)
|
|
{
|
|
f [0] = x [0] + x [1];
|
|
f [1] = x [0] - x [1];
|
|
}
|
|
|
|
// 1-point FFT
|
|
else
|
|
{
|
|
f [0] = x [0];
|
|
}
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
==============================================================================
|
|
Name: do_ifft
|
|
Description:
|
|
Compute the inverse FFT of the array. Note that data must be post-scaled:
|
|
IFFT (FFT (x)) = x * length (x).
|
|
Input parameters:
|
|
- f: pointer on the source array (frequencies).
|
|
f [0...length(x)/2] = real values
|
|
f [length(x)/2+1...length(x)-1] = negative imaginary values of
|
|
coefficents 1...length(x)/2-1.
|
|
Output parameters:
|
|
- x: pointer on the destination array (time).
|
|
Throws: Nothing
|
|
==============================================================================
|
|
*/
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::do_ifft (const DataType f [], DataType x []) const
|
|
{
|
|
assert (f != 0);
|
|
assert (f != use_buffer ());
|
|
assert (x != 0);
|
|
assert (x != use_buffer ());
|
|
assert (x != f);
|
|
|
|
// General case
|
|
if (_nbr_bits > 2)
|
|
{
|
|
compute_ifft_general (f, x);
|
|
}
|
|
|
|
// 4-point IFFT
|
|
else if (_nbr_bits == 2)
|
|
{
|
|
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
|
|
else if (_nbr_bits == 1)
|
|
{
|
|
x [0] = f [0] + f [1];
|
|
x [1] = f [0] - f [1];
|
|
}
|
|
|
|
// 1-point IFFT
|
|
else
|
|
{
|
|
x [0] = f [0];
|
|
}
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
==============================================================================
|
|
Name: rescale
|
|
Description:
|
|
Scale an array by divide each element by its length. This function should
|
|
be called after FFT + IFFT.
|
|
Input parameters:
|
|
- x: pointer on array to rescale (time or frequency).
|
|
Throws: Nothing
|
|
==============================================================================
|
|
*/
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::rescale (DataType x []) const
|
|
{
|
|
const DataType mul = DataType (1.0 / _length);
|
|
|
|
if (_length < 4)
|
|
{
|
|
long i = _length - 1;
|
|
do
|
|
{
|
|
x [i] *= mul;
|
|
--i;
|
|
}
|
|
while (i >= 0);
|
|
}
|
|
|
|
else
|
|
{
|
|
assert ((_length & 3) == 0);
|
|
|
|
// Could be optimized with SIMD instruction sets (needs alignment check)
|
|
long i = _length - 4;
|
|
do
|
|
{
|
|
x [i + 0] *= mul;
|
|
x [i + 1] *= mul;
|
|
x [i + 2] *= mul;
|
|
x [i + 3] *= mul;
|
|
i -= 4;
|
|
}
|
|
while (i >= 0);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
/*
|
|
==============================================================================
|
|
Name: use_buffer
|
|
Description:
|
|
Access the internal buffer, whose length is the FFT one.
|
|
Buffer content will be erased at each do_fft() / do_ifft() call!
|
|
This buffer cannot be used as:
|
|
- source for FFT or IFFT done with this object
|
|
- destination for FFT or IFFT done with this object
|
|
Returns:
|
|
Buffer start address
|
|
Throws: Nothing
|
|
==============================================================================
|
|
*/
|
|
|
|
template <class DT>
|
|
typename FFTReal <DT>::DataType * FFTReal <DT>::use_buffer () const
|
|
{
|
|
return (&_buffer [0]);
|
|
}
|
|
|
|
|
|
|
|
/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
|
|
|
|
|
|
|
|
/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::init_br_lut ()
|
|
{
|
|
const long length = 1L << _nbr_bits;
|
|
_br_lut.resize (length);
|
|
|
|
_br_lut [0] = 0;
|
|
long br_index = 0;
|
|
for (long cnt = 1; cnt < length; ++cnt)
|
|
{
|
|
// ++br_index (bit reversed)
|
|
long bit = length >> 1;
|
|
while (((br_index ^= bit) & bit) == 0)
|
|
{
|
|
bit >>= 1;
|
|
}
|
|
|
|
_br_lut [cnt] = br_index;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::init_trigo_lut ()
|
|
{
|
|
using namespace std;
|
|
|
|
if (_nbr_bits > 3)
|
|
{
|
|
const long total_len = (1L << (_nbr_bits - 1)) - 4;
|
|
_trigo_lut.resize (total_len);
|
|
|
|
for (int level = 3; level < _nbr_bits; ++level)
|
|
{
|
|
const long level_len = 1L << (level - 1);
|
|
DataType * const level_ptr =
|
|
&_trigo_lut [get_trigo_level_index (level)];
|
|
const double mul = PI / (level_len << 1);
|
|
|
|
for (long i = 0; i < level_len; ++ i)
|
|
{
|
|
level_ptr [i] = static_cast <DataType> (cos (i * mul));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::init_trigo_osc ()
|
|
{
|
|
const int nbr_osc = _nbr_bits - TRIGO_BD_LIMIT;
|
|
if (nbr_osc > 0)
|
|
{
|
|
_trigo_osc.resize (nbr_osc);
|
|
|
|
for (int osc_cnt = 0; osc_cnt < nbr_osc; ++osc_cnt)
|
|
{
|
|
OscType & osc = _trigo_osc [osc_cnt];
|
|
|
|
const long len = 1L << (TRIGO_BD_LIMIT + osc_cnt);
|
|
const double mul = (0.5 * PI) / len;
|
|
osc.set_step (mul);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
const long * FFTReal <DT>::get_br_ptr () const
|
|
{
|
|
return (&_br_lut [0]);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
const typename FFTReal <DT>::DataType * FFTReal <DT>::get_trigo_ptr (int level) const
|
|
{
|
|
assert (level >= 3);
|
|
|
|
return (&_trigo_lut [get_trigo_level_index (level)]);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
long FFTReal <DT>::get_trigo_level_index (int level) const
|
|
{
|
|
assert (level >= 3);
|
|
|
|
return ((1L << (level - 1)) - 4);
|
|
}
|
|
|
|
|
|
|
|
// Transform in several passes
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_fft_general (DataType f [], const DataType x []) const
|
|
{
|
|
assert (f != 0);
|
|
assert (f != use_buffer ());
|
|
assert (x != 0);
|
|
assert (x != use_buffer ());
|
|
assert (x != f);
|
|
|
|
DataType * sf;
|
|
DataType * df;
|
|
|
|
if ((_nbr_bits & 1) != 0)
|
|
{
|
|
df = use_buffer ();
|
|
sf = f;
|
|
}
|
|
else
|
|
{
|
|
df = f;
|
|
sf = use_buffer ();
|
|
}
|
|
|
|
compute_direct_pass_1_2 (df, x);
|
|
compute_direct_pass_3 (sf, df);
|
|
|
|
for (int pass = 3; pass < _nbr_bits; ++ pass)
|
|
{
|
|
compute_direct_pass_n (df, sf, pass);
|
|
|
|
DataType * const temp_ptr = df;
|
|
df = sf;
|
|
sf = temp_ptr;
|
|
}
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_direct_pass_1_2 (DataType df [], const DataType x []) const
|
|
{
|
|
assert (df != 0);
|
|
assert (x != 0);
|
|
assert (df != x);
|
|
|
|
const long * const bit_rev_lut_ptr = get_br_ptr ();
|
|
long coef_index = 0;
|
|
do
|
|
{
|
|
const long rev_index_0 = bit_rev_lut_ptr [coef_index];
|
|
const long rev_index_1 = bit_rev_lut_ptr [coef_index + 1];
|
|
const long rev_index_2 = bit_rev_lut_ptr [coef_index + 2];
|
|
const long rev_index_3 = bit_rev_lut_ptr [coef_index + 3];
|
|
|
|
DataType * const df2 = df + coef_index;
|
|
df2 [1] = x [rev_index_0] - x [rev_index_1];
|
|
df2 [3] = x [rev_index_2] - x [rev_index_3];
|
|
|
|
const DataType sf_0 = x [rev_index_0] + x [rev_index_1];
|
|
const DataType sf_2 = x [rev_index_2] + x [rev_index_3];
|
|
|
|
df2 [0] = sf_0 + sf_2;
|
|
df2 [2] = sf_0 - sf_2;
|
|
|
|
coef_index += 4;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_direct_pass_3 (DataType df [], const DataType sf []) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
|
|
const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
|
|
long coef_index = 0;
|
|
do
|
|
{
|
|
DataType v;
|
|
|
|
df [coef_index] = sf [coef_index] + sf [coef_index + 4];
|
|
df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
|
|
df [coef_index + 2] = sf [coef_index + 2];
|
|
df [coef_index + 6] = sf [coef_index + 6];
|
|
|
|
v = (sf [coef_index + 5] - sf [coef_index + 7]) * sqrt2_2;
|
|
df [coef_index + 1] = sf [coef_index + 1] + v;
|
|
df [coef_index + 3] = sf [coef_index + 1] - v;
|
|
|
|
v = (sf [coef_index + 5] + sf [coef_index + 7]) * sqrt2_2;
|
|
df [coef_index + 5] = v + sf [coef_index + 3];
|
|
df [coef_index + 7] = v - sf [coef_index + 3];
|
|
|
|
coef_index += 8;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_direct_pass_n (DataType df [], const DataType sf [], int pass) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
assert (pass >= 3);
|
|
assert (pass < _nbr_bits);
|
|
|
|
if (pass <= TRIGO_BD_LIMIT)
|
|
{
|
|
compute_direct_pass_n_lut (df, sf, pass);
|
|
}
|
|
else
|
|
{
|
|
compute_direct_pass_n_osc (df, sf, pass);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_direct_pass_n_lut (DataType df [], const DataType sf [], int pass) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
assert (pass >= 3);
|
|
assert (pass < _nbr_bits);
|
|
|
|
const long nbr_coef = 1 << pass;
|
|
const long h_nbr_coef = nbr_coef >> 1;
|
|
const long d_nbr_coef = nbr_coef << 1;
|
|
long coef_index = 0;
|
|
const DataType * const cos_ptr = get_trigo_ptr (pass);
|
|
do
|
|
{
|
|
const DataType * const sf1r = sf + coef_index;
|
|
const DataType * const sf2r = sf1r + nbr_coef;
|
|
DataType * const dfr = df + coef_index;
|
|
DataType * const dfi = dfr + nbr_coef;
|
|
|
|
// Extreme coefficients are always real
|
|
dfr [0] = sf1r [0] + sf2r [0];
|
|
dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] =
|
|
dfr [h_nbr_coef] = sf1r [h_nbr_coef];
|
|
dfi [h_nbr_coef] = sf2r [h_nbr_coef];
|
|
|
|
// Others are conjugate complex numbers
|
|
const DataType * const sf1i = sf1r + h_nbr_coef;
|
|
const DataType * const sf2i = sf1i + nbr_coef;
|
|
for (long i = 1; i < h_nbr_coef; ++ i)
|
|
{
|
|
const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef);
|
|
const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
|
|
DataType v;
|
|
|
|
v = sf2r [i] * c - sf2i [i] * s;
|
|
dfr [i] = sf1r [i] + v;
|
|
dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] =
|
|
|
|
v = sf2r [i] * s + sf2i [i] * c;
|
|
dfi [i] = v + sf1i [i];
|
|
dfi [nbr_coef - i] = v - sf1i [i];
|
|
}
|
|
|
|
coef_index += d_nbr_coef;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_direct_pass_n_osc (DataType df [], const DataType sf [], int pass) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
assert (pass > TRIGO_BD_LIMIT);
|
|
assert (pass < _nbr_bits);
|
|
|
|
const long nbr_coef = 1 << pass;
|
|
const long h_nbr_coef = nbr_coef >> 1;
|
|
const long d_nbr_coef = nbr_coef << 1;
|
|
long coef_index = 0;
|
|
OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
|
|
do
|
|
{
|
|
const DataType * const sf1r = sf + coef_index;
|
|
const DataType * const sf2r = sf1r + nbr_coef;
|
|
DataType * const dfr = df + coef_index;
|
|
DataType * const dfi = dfr + nbr_coef;
|
|
|
|
osc.clear_buffers ();
|
|
|
|
// Extreme coefficients are always real
|
|
dfr [0] = sf1r [0] + sf2r [0];
|
|
dfi [0] = sf1r [0] - sf2r [0]; // dfr [nbr_coef] =
|
|
dfr [h_nbr_coef] = sf1r [h_nbr_coef];
|
|
dfi [h_nbr_coef] = sf2r [h_nbr_coef];
|
|
|
|
// Others are conjugate complex numbers
|
|
const DataType * const sf1i = sf1r + h_nbr_coef;
|
|
const DataType * const sf2i = sf1i + nbr_coef;
|
|
for (long i = 1; i < h_nbr_coef; ++ i)
|
|
{
|
|
osc.step ();
|
|
const DataType c = osc.get_cos ();
|
|
const DataType s = osc.get_sin ();
|
|
DataType v;
|
|
|
|
v = sf2r [i] * c - sf2i [i] * s;
|
|
dfr [i] = sf1r [i] + v;
|
|
dfi [-i] = sf1r [i] - v; // dfr [nbr_coef - i] =
|
|
|
|
v = sf2r [i] * s + sf2i [i] * c;
|
|
dfi [i] = v + sf1i [i];
|
|
dfi [nbr_coef - i] = v - sf1i [i];
|
|
}
|
|
|
|
coef_index += d_nbr_coef;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
// Transform in several pass
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_ifft_general (const DataType f [], DataType x []) const
|
|
{
|
|
assert (f != 0);
|
|
assert (f != use_buffer ());
|
|
assert (x != 0);
|
|
assert (x != use_buffer ());
|
|
assert (x != f);
|
|
|
|
DataType * sf = const_cast <DataType *> (f);
|
|
DataType * df;
|
|
DataType * df_temp;
|
|
|
|
if (_nbr_bits & 1)
|
|
{
|
|
df = use_buffer ();
|
|
df_temp = x;
|
|
}
|
|
else
|
|
{
|
|
df = x;
|
|
df_temp = use_buffer ();
|
|
}
|
|
|
|
for (int pass = _nbr_bits - 1; pass >= 3; -- pass)
|
|
{
|
|
compute_inverse_pass_n (df, sf, pass);
|
|
|
|
if (pass < _nbr_bits - 1)
|
|
{
|
|
DataType * const temp_ptr = df;
|
|
df = sf;
|
|
sf = temp_ptr;
|
|
}
|
|
else
|
|
{
|
|
sf = df;
|
|
df = df_temp;
|
|
}
|
|
}
|
|
|
|
compute_inverse_pass_3 (df, sf);
|
|
compute_inverse_pass_1_2 (x, df);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_inverse_pass_n (DataType df [], const DataType sf [], int pass) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
assert (pass >= 3);
|
|
assert (pass < _nbr_bits);
|
|
|
|
if (pass <= TRIGO_BD_LIMIT)
|
|
{
|
|
compute_inverse_pass_n_lut (df, sf, pass);
|
|
}
|
|
else
|
|
{
|
|
compute_inverse_pass_n_osc (df, sf, pass);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_inverse_pass_n_lut (DataType df [], const DataType sf [], int pass) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
assert (pass >= 3);
|
|
assert (pass < _nbr_bits);
|
|
|
|
const long nbr_coef = 1 << pass;
|
|
const long h_nbr_coef = nbr_coef >> 1;
|
|
const long d_nbr_coef = nbr_coef << 1;
|
|
long coef_index = 0;
|
|
const DataType * const cos_ptr = get_trigo_ptr (pass);
|
|
do
|
|
{
|
|
const DataType * const sfr = sf + coef_index;
|
|
const DataType * const sfi = sfr + nbr_coef;
|
|
DataType * const df1r = df + coef_index;
|
|
DataType * const df2r = df1r + nbr_coef;
|
|
|
|
// Extreme coefficients are always real
|
|
df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef]
|
|
df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef]
|
|
df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
|
|
df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
|
|
|
|
// Others are conjugate complex numbers
|
|
DataType * const df1i = df1r + h_nbr_coef;
|
|
DataType * const df2i = df1i + nbr_coef;
|
|
for (long i = 1; i < h_nbr_coef; ++ i)
|
|
{
|
|
df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i]
|
|
df1i [i] = sfi [i] - sfi [nbr_coef - i];
|
|
|
|
const DataType c = cos_ptr [i]; // cos (i*PI/nbr_coef);
|
|
const DataType s = cos_ptr [h_nbr_coef - i]; // sin (i*PI/nbr_coef);
|
|
const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i]
|
|
const DataType vi = sfi [i] + sfi [nbr_coef - i];
|
|
|
|
df2r [i] = vr * c + vi * s;
|
|
df2i [i] = vi * c - vr * s;
|
|
}
|
|
|
|
coef_index += d_nbr_coef;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_inverse_pass_n_osc (DataType df [], const DataType sf [], int pass) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
assert (pass > TRIGO_BD_LIMIT);
|
|
assert (pass < _nbr_bits);
|
|
|
|
const long nbr_coef = 1 << pass;
|
|
const long h_nbr_coef = nbr_coef >> 1;
|
|
const long d_nbr_coef = nbr_coef << 1;
|
|
long coef_index = 0;
|
|
OscType & osc = _trigo_osc [pass - (TRIGO_BD_LIMIT + 1)];
|
|
do
|
|
{
|
|
const DataType * const sfr = sf + coef_index;
|
|
const DataType * const sfi = sfr + nbr_coef;
|
|
DataType * const df1r = df + coef_index;
|
|
DataType * const df2r = df1r + nbr_coef;
|
|
|
|
osc.clear_buffers ();
|
|
|
|
// Extreme coefficients are always real
|
|
df1r [0] = sfr [0] + sfi [0]; // + sfr [nbr_coef]
|
|
df2r [0] = sfr [0] - sfi [0]; // - sfr [nbr_coef]
|
|
df1r [h_nbr_coef] = sfr [h_nbr_coef] * 2;
|
|
df2r [h_nbr_coef] = sfi [h_nbr_coef] * 2;
|
|
|
|
// Others are conjugate complex numbers
|
|
DataType * const df1i = df1r + h_nbr_coef;
|
|
DataType * const df2i = df1i + nbr_coef;
|
|
for (long i = 1; i < h_nbr_coef; ++ i)
|
|
{
|
|
df1r [i] = sfr [i] + sfi [-i]; // + sfr [nbr_coef - i]
|
|
df1i [i] = sfi [i] - sfi [nbr_coef - i];
|
|
|
|
osc.step ();
|
|
const DataType c = osc.get_cos ();
|
|
const DataType s = osc.get_sin ();
|
|
const DataType vr = sfr [i] - sfi [-i]; // - sfr [nbr_coef - i]
|
|
const DataType vi = sfi [i] + sfi [nbr_coef - i];
|
|
|
|
df2r [i] = vr * c + vi * s;
|
|
df2i [i] = vi * c - vr * s;
|
|
}
|
|
|
|
coef_index += d_nbr_coef;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_inverse_pass_3 (DataType df [], const DataType sf []) const
|
|
{
|
|
assert (df != 0);
|
|
assert (sf != 0);
|
|
assert (df != sf);
|
|
|
|
const DataType sqrt2_2 = DataType (SQRT2 * 0.5);
|
|
long coef_index = 0;
|
|
do
|
|
{
|
|
df [coef_index] = sf [coef_index] + sf [coef_index + 4];
|
|
df [coef_index + 4] = sf [coef_index] - sf [coef_index + 4];
|
|
df [coef_index + 2] = sf [coef_index + 2] * 2;
|
|
df [coef_index + 6] = sf [coef_index + 6] * 2;
|
|
|
|
df [coef_index + 1] = sf [coef_index + 1] + sf [coef_index + 3];
|
|
df [coef_index + 3] = sf [coef_index + 5] - sf [coef_index + 7];
|
|
|
|
const DataType vr = sf [coef_index + 1] - sf [coef_index + 3];
|
|
const DataType vi = sf [coef_index + 5] + sf [coef_index + 7];
|
|
|
|
df [coef_index + 5] = (vr + vi) * sqrt2_2;
|
|
df [coef_index + 7] = (vi - vr) * sqrt2_2;
|
|
|
|
coef_index += 8;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
template <class DT>
|
|
void FFTReal <DT>::compute_inverse_pass_1_2 (DataType x [], const DataType sf []) const
|
|
{
|
|
assert (x != 0);
|
|
assert (sf != 0);
|
|
assert (x != sf);
|
|
|
|
const long * bit_rev_lut_ptr = get_br_ptr ();
|
|
const DataType * sf2 = sf;
|
|
long coef_index = 0;
|
|
do
|
|
{
|
|
{
|
|
const DataType b_0 = sf2 [0] + sf2 [2];
|
|
const DataType b_2 = sf2 [0] - sf2 [2];
|
|
const DataType b_1 = sf2 [1] * 2;
|
|
const DataType b_3 = sf2 [3] * 2;
|
|
|
|
x [bit_rev_lut_ptr [0]] = b_0 + b_1;
|
|
x [bit_rev_lut_ptr [1]] = b_0 - b_1;
|
|
x [bit_rev_lut_ptr [2]] = b_2 + b_3;
|
|
x [bit_rev_lut_ptr [3]] = b_2 - b_3;
|
|
}
|
|
{
|
|
const DataType b_0 = sf2 [4] + sf2 [6];
|
|
const DataType b_2 = sf2 [4] - sf2 [6];
|
|
const DataType b_1 = sf2 [5] * 2;
|
|
const DataType b_3 = sf2 [7] * 2;
|
|
|
|
x [bit_rev_lut_ptr [4]] = b_0 + b_1;
|
|
x [bit_rev_lut_ptr [5]] = b_0 - b_1;
|
|
x [bit_rev_lut_ptr [6]] = b_2 + b_3;
|
|
x [bit_rev_lut_ptr [7]] = b_2 - b_3;
|
|
}
|
|
|
|
sf2 += 8;
|
|
coef_index += 8;
|
|
bit_rev_lut_ptr += 8;
|
|
}
|
|
while (coef_index < _length);
|
|
}
|
|
|
|
|
|
|
|
#endif // FFTReal_CODEHEADER_INCLUDED
|
|
|
|
#undef FFTReal_CURRENT_CODEHEADER
|
|
|
|
|
|
|
|
/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/
|