diff options
Diffstat (limited to 'demos/spectrum/3rdparty/fftreal/FFTRealPassInverse.hpp')
| -rw-r--r-- | demos/spectrum/3rdparty/fftreal/FFTRealPassInverse.hpp | 229 |
1 files changed, 229 insertions, 0 deletions
diff --git a/demos/spectrum/3rdparty/fftreal/FFTRealPassInverse.hpp b/demos/spectrum/3rdparty/fftreal/FFTRealPassInverse.hpp new file mode 100644 index 0000000..5737546 --- /dev/null +++ b/demos/spectrum/3rdparty/fftreal/FFTRealPassInverse.hpp @@ -0,0 +1,229 @@ +/***************************************************************************** + + FFTRealPassInverse.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 (FFTRealPassInverse_CURRENT_CODEHEADER) + #error Recursive inclusion of FFTRealPassInverse code header. +#endif +#define FFTRealPassInverse_CURRENT_CODEHEADER + +#if ! defined (FFTRealPassInverse_CODEHEADER_INCLUDED) +#define FFTRealPassInverse_CODEHEADER_INCLUDED + + + +/*\\\ INCLUDE FILES \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + +#include "FFTRealUseTrigo.h" + + + +/*\\\ PUBLIC \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +template <int PASS> +void FFTRealPassInverse <PASS>::process (long len, DataType dest_ptr [], DataType src_ptr [], const DataType f_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) +{ + process_internal ( + len, + dest_ptr, + f_ptr, + cos_ptr, + cos_len, + br_ptr, + osc_list + ); + FFTRealPassInverse <PASS - 1>::process_rec ( + len, + src_ptr, + dest_ptr, + cos_ptr, + cos_len, + br_ptr, + osc_list + ); +} + + + +template <int PASS> +void FFTRealPassInverse <PASS>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) +{ + process_internal ( + len, + dest_ptr, + src_ptr, + cos_ptr, + cos_len, + br_ptr, + osc_list + ); + FFTRealPassInverse <PASS - 1>::process_rec ( + len, + src_ptr, + dest_ptr, + cos_ptr, + cos_len, + br_ptr, + osc_list + ); +} + +template <> +void FFTRealPassInverse <0>::process_rec (long len, DataType dest_ptr [], DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) +{ + // Stops recursion +} + + + +template <int PASS> +void FFTRealPassInverse <PASS>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) +{ + const long dist = 1L << (PASS - 1); + const long c1_r = 0; + const long c1_i = dist; + const long c2_r = dist * 2; + const long c2_i = dist * 3; + const long cend = dist * 4; + const long table_step = cos_len >> (PASS - 1); + + enum { TRIGO_OSC = PASS - FFTRealFixLenParam::TRIGO_BD_LIMIT }; + enum { TRIGO_DIRECT = (TRIGO_OSC >= 0) ? 1 : 0 }; + + long coef_index = 0; + do + { + const DataType * const sf = src_ptr + coef_index; + DataType * const df = dest_ptr + coef_index; + + // Extreme coefficients are always real + df [c1_r] = sf [c1_r] + sf [c2_r]; + df [c2_r] = sf [c1_r] - sf [c2_r]; + df [c1_i] = sf [c1_i] * 2; + df [c2_i] = sf [c2_i] * 2; + + FFTRealUseTrigo <TRIGO_DIRECT>::prepare (osc_list [TRIGO_OSC]); + + // Others are conjugate complex numbers + for (long i = 1; i < dist; ++ i) + { + df [c1_r + i] = sf [c1_r + i] + sf [c2_r - i]; + df [c1_i + i] = sf [c2_r + i] - sf [cend - i]; + + DataType c; + DataType s; + FFTRealUseTrigo <TRIGO_DIRECT>::iterate ( + osc_list [TRIGO_OSC], + c, + s, + cos_ptr, + i * table_step, + (dist - i) * table_step + ); + + const DataType vr = sf [c1_r + i] - sf [c2_r - i]; + const DataType vi = sf [c2_r + i] + sf [cend - i]; + + df [c2_r + i] = vr * c + vi * s; + df [c2_i + i] = vi * c - vr * s; + } + + coef_index += cend; + } + while (coef_index < len); +} + +template <> +void FFTRealPassInverse <2>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) +{ + // Antepenultimate pass + const DataType sqrt2_2 = DataType (SQRT2 * 0.5); + + long coef_index = 0; + do + { + dest_ptr [coef_index ] = src_ptr [coef_index] + src_ptr [coef_index + 4]; + dest_ptr [coef_index + 4] = src_ptr [coef_index] - src_ptr [coef_index + 4]; + dest_ptr [coef_index + 2] = src_ptr [coef_index + 2] * 2; + dest_ptr [coef_index + 6] = src_ptr [coef_index + 6] * 2; + + dest_ptr [coef_index + 1] = src_ptr [coef_index + 1] + src_ptr [coef_index + 3]; + dest_ptr [coef_index + 3] = src_ptr [coef_index + 5] - src_ptr [coef_index + 7]; + + const DataType vr = src_ptr [coef_index + 1] - src_ptr [coef_index + 3]; + const DataType vi = src_ptr [coef_index + 5] + src_ptr [coef_index + 7]; + + dest_ptr [coef_index + 5] = (vr + vi) * sqrt2_2; + dest_ptr [coef_index + 7] = (vi - vr) * sqrt2_2; + + coef_index += 8; + } + while (coef_index < len); +} + +template <> +void FFTRealPassInverse <1>::process_internal (long len, DataType dest_ptr [], const DataType src_ptr [], const DataType cos_ptr [], long cos_len, const long br_ptr [], OscType osc_list []) +{ + // Penultimate and last pass at once + const long qlen = len >> 2; + + long coef_index = 0; + do + { + const long ri_0 = br_ptr [coef_index >> 2]; + + const DataType b_0 = src_ptr [coef_index ] + src_ptr [coef_index + 2]; + const DataType b_2 = src_ptr [coef_index ] - src_ptr [coef_index + 2]; + const DataType b_1 = src_ptr [coef_index + 1] * 2; + const DataType b_3 = src_ptr [coef_index + 3] * 2; + + dest_ptr [ri_0 ] = b_0 + b_1; + dest_ptr [ri_0 + 2 * qlen] = b_0 - b_1; + dest_ptr [ri_0 + 1 * qlen] = b_2 + b_3; + dest_ptr [ri_0 + 3 * qlen] = b_2 - b_3; + + coef_index += 4; + } + while (coef_index < len); +} + + + +/*\\\ PROTECTED \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +/*\\\ PRIVATE \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ + + + +#endif // FFTRealPassInverse_CODEHEADER_INCLUDED + +#undef FFTRealPassInverse_CURRENT_CODEHEADER + + + +/*\\\ EOF \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\*/ |