From 38c835feaf1ff5d30a8c2db94e2951f0839d8812 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sat, 9 Apr 2016 13:23:43 +1000 Subject: [PATCH 01/12] First crack at FFT: problem with iterator_adaptor. --- include/boost/math/tools/fft.hpp | 118 +++++++++++++++++++++++++++++++++++++++ test/test_fft.cpp | 29 ++++++++++ 2 files changed, 147 insertions(+) create mode 100644 include/boost/math/tools/fft.hpp create mode 100644 test/test_fft.cpp diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp new file mode 100644 index 000000000..2fc24178f --- /dev/null +++ b/include/boost/math/tools/fft.hpp @@ -0,0 +1,118 @@ +// (C) Copyright Jeremy William Murphy 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_TOOLS_FFT_HPP +#define BOOST_MATH_TOOLS_FFT_HPP + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace boost { namespace math { namespace tools { + +namespace detail +{ + template + class stride_iterator : public iterator_adaptor, Iterator> + { + typedef iterator_adaptor, Iterator> super_t; + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type base_distance; + friend class iterator_core_access; + + base_distance stride; + + public: + stride_iterator() {} + + stride_iterator(Iterator x, base_distance stride) : super_t(x), stride(stride) {} + + template + stride_iterator(stride_iterator const &x, typename enable_if_convertible::type* = 0) : super_t(x.base()) + {} + + template + stride_iterator(stride_iterator const &x, base_distance stride, typename enable_if_convertible::type* = 0) : super_t(x.base()), stride(x.stride * stride) + {} + + private: + void increment() { std::advance(this->base_reference(), stride); } + void decrement() { std::advance(this->base_reference(), -stride); } + + template + BOOST_DEDUCED_TYPENAME super_t::difference_type + distance_to(iterator_adaptor const& y) const + { + return (y.base() - this->base_reference()) / stride; + } + }; + + template + stride_iterator make_stride_iterator(stride_iterator x, N stride) + { + return stride_iterator(x, stride); + } + + // BOOST_DEDUCED_TYPENAME disable_if_c >::value, stride_iterator >::type + + template + stride_iterator make_stride_iterator(Iterator x, N stride) + { + return stride_iterator(x, stride); + } +} + +template +void inner_fft(I y_0, I y_1, J y, N n, CRU &w, CRU const &w_n) +{ + for (unsigned k = 0; k < n / 2 - 1; k++) + { + y[k] = y_0[k] + w * y_1[k]; + y[k + n / 2] = y_0[k] - w * y_1[k]; + w *= w_n; + } +} + +template +void fft(I first, I last, J y) +{ + using namespace boost::math::tools::detail; + + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type difference_type; + typedef std::complex complex_t; + difference_type const n = distance(first, last); + BOOST_ASSERT((n & (n - 1)) == 0); // We require that n is a power of 2. + if (n == 1) + return; + complex_t const i(0, 1); + complex_t const w_n = std::exp((2.0 * M_PI * i) / static_cast(n)); + complex_t w = 1; + std::vector y_0, y_1; + y_0.resize(n / 2); + y_1.resize(n / 2); + + fft(make_stride_iterator(first, 2), make_stride_iterator(last - 1, 2), y_0.begin()); + fft(make_stride_iterator(first + 1, 2), make_stride_iterator(last, 2), y_1.begin()); + + for (unsigned k = 0; k < n / 2 - 1; k++) + { + y[k] = y_0[k] + w * y_1[k]; + y[k + n / 2] = y_0[k] - w * y_1[k]; + w *= w_n; + } +} + +}}} + +#endif + + diff --git a/test/test_fft.cpp b/test/test_fft.cpp new file mode 100644 index 000000000..79864b3dc --- /dev/null +++ b/test/test_fft.cpp @@ -0,0 +1,29 @@ +// (C) Copyright Jeremy Murphy 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#define BOOST_TEST_MAIN + +#include +#include +#include +#include +#include +#include + +typedef boost::mpl::list signed_integral_test_types; +// , boost::multiprecision::cpp_int + +using namespace boost; +using namespace boost::math::tools; + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_basic, T, signed_integral_test_types) +{ + std::vector a; + std::copy(make_counting_iterator(1), make_counting_iterator(5), std::back_inserter(a)); + std::vector< std::complex > result; + result.resize(4); + fft(a.begin(), a.end(), result.begin()); +} From a38deb9789f025e54d4c7ec23774de675726a571 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sat, 9 Apr 2016 13:55:39 +1000 Subject: [PATCH 02/12] Work around the compile error by just making them public. :\ --- include/boost/math/tools/fft.hpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index 2fc24178f..b8ecf5220 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -44,16 +44,22 @@ namespace detail stride_iterator(stride_iterator const &x, base_distance stride, typename enable_if_convertible::type* = 0) : super_t(x.base()), stride(x.stride * stride) {} - private: - void increment() { std::advance(this->base_reference(), stride); } - void decrement() { std::advance(this->base_reference(), -stride); } - + // TODO: These *should* be private but compilation fails unexpectedly. template BOOST_DEDUCED_TYPENAME super_t::difference_type distance_to(iterator_adaptor const& y) const { return (y.base() - this->base_reference()) / stride; } + + void advance(BOOST_DEDUCED_TYPENAME super_t::difference_type n) + { + std::advance(this->base_reference(), stride * n); + } + private: + void increment() { std::advance(this->base_reference(), stride); } + void decrement() { std::advance(this->base_reference(), -stride); } + }; template @@ -62,8 +68,6 @@ namespace detail return stride_iterator(x, stride); } - // BOOST_DEDUCED_TYPENAME disable_if_c >::value, stride_iterator >::type - template stride_iterator make_stride_iterator(Iterator x, N stride) { From 718e491c73a5242add50f9e9e24d511986d3909b Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 10 Apr 2016 00:03:23 +1000 Subject: [PATCH 03/12] Clean up recursive_fft, add iterative_fft. --- include/boost/math/tools/fft.hpp | 97 ++++++++++++++++++++++++++++++++-------- test/test_fft.cpp | 60 ++++++++++++++++++++++--- 2 files changed, 134 insertions(+), 23 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index b8ecf5220..bcc01cf0b 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -11,7 +11,6 @@ #include #include #include -#include #include #include @@ -75,28 +74,22 @@ namespace detail } } -template -void inner_fft(I y_0, I y_1, J y, N n, CRU &w, CRU const &w_n) -{ - for (unsigned k = 0; k < n / 2 - 1; k++) - { - y[k] = y_0[k] + w * y_1[k]; - y[k + n / 2] = y_0[k] - w * y_1[k]; - w *= w_n; - } -} +// Recursive FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 template -void fft(I first, I last, J y) +void recursive_fft(I first, I last, J y) { using namespace boost::math::tools::detail; - typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type difference_type; - typedef std::complex complex_t; - difference_type const n = distance(first, last); + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type N; + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type complex_t; + N const n = distance(first, last); BOOST_ASSERT((n & (n - 1)) == 0); // We require that n is a power of 2. if (n == 1) + { + *y = *first; return; + } complex_t const i(0, 1); complex_t const w_n = std::exp((2.0 * M_PI * i) / static_cast(n)); complex_t w = 1; @@ -104,10 +97,10 @@ void fft(I first, I last, J y) y_0.resize(n / 2); y_1.resize(n / 2); - fft(make_stride_iterator(first, 2), make_stride_iterator(last - 1, 2), y_0.begin()); - fft(make_stride_iterator(first + 1, 2), make_stride_iterator(last, 2), y_1.begin()); + recursive_fft(make_stride_iterator(first, 2), make_stride_iterator(last, 2), y_0.begin()); + recursive_fft(make_stride_iterator(first + 1, 2), make_stride_iterator(last + 1, 2), y_1.begin()); - for (unsigned k = 0; k < n / 2 - 1; k++) + for (N k = 0; k < n / 2; k++) { y[k] = y_0[k] + w * y_1[k]; y[k + n / 2] = y_0[k] - w * y_1[k]; @@ -115,6 +108,74 @@ void fft(I first, I last, J y) } } + +// Based on code from Henry Warren Jr, Hacker's Delight, 2nd ed., 2013, p 138 +template +void increment_reversed_integer(N &x, N m) +{ + x ^= m; + if (x <= m) + { + do + { + m >>= 1; + x ^= m; + } + while (x < m); + } +} + + +// Presumably this can be done more elegantly (and efficiently) in combination +// with a permutation_iterator, etc. +template +void bit_reversed_permute(I first, I last) +{ + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type N; + N const n = std::distance(first, last); + BOOST_ASSERT((n & (n - 1)) == 0); // This might not be a requirement. + + N const m = n / 2; + N i_rev = m; + for (I i = first + 1; i != last - 1; i++) + { + if (i_rev > i - first) + std::iter_swap(i, first + i_rev); + increment_reversed_integer(i_rev, m); + } +} + + +// Iterative FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 +template +void iterative_fft(I first, I last, J A) +{ + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type complex_t; + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type N; + std::copy(first, last, A); + N const n = std::distance(first, last); + bit_reversed_permute(A, A + n); + complex_t const i(0, 1); + N log2_n = log2(n); // TODO: Err... C99? + for (N s = 1; s <= log2_n; s++) + { + N const m = N(1) << s; + complex_t const w_m = std::exp((2.0 * M_PI * i) / static_cast(m)); + for (N k = 0; k < n; k += m) + { + complex_t w = 1; + N const m_2 = m / 2; + for (N j = 0; j < m_2; j++) + { + complex_t const t = w * A[k + j + m_2]; + A[k + j + m_2] = A[k + j] - t; + A[k + j] += t; + w *= w_m; + } + } + } +} + }}} #endif diff --git a/test/test_fft.cpp b/test/test_fft.cpp index 79864b3dc..d74cf3eb8 100644 --- a/test/test_fft.cpp +++ b/test/test_fft.cpp @@ -6,24 +6,74 @@ #include #define BOOST_TEST_MAIN +#include #include #include #include #include #include +#include + +#include + #include +#include +#include -typedef boost::mpl::list signed_integral_test_types; +typedef boost::mpl::list signed_integral_test_types; // , boost::multiprecision::cpp_int using namespace boost; using namespace boost::math::tools; -BOOST_AUTO_TEST_CASE_TEMPLATE(test_basic, T, signed_integral_test_types) +typedef std::complex complex_t; + +// TODO: Fixtures. + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_recursive_fft, T, signed_integral_test_types) +{ + std::vector a; + std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); + std::vector< std::complex > result; + result.resize(16); + recursive_fft(a.begin(), a.end(), result.begin()); + // TODO: Compare against benchmark result (GSL). +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_bit_reversed_permute, T, signed_integral_test_types) +{ + std::vector a; + std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); + bit_reversed_permute(a.begin(), a.end()); + boost::array b = {0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15}; + BOOST_CHECK_EQUAL_COLLECTIONS(boost::begin(a), boost::end(a), boost::begin(b), boost::end(b)); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_iterative_fft, T, signed_integral_test_types) { std::vector a; - std::copy(make_counting_iterator(1), make_counting_iterator(5), std::back_inserter(a)); + std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); std::vector< std::complex > result; - result.resize(4); - fft(a.begin(), a.end(), result.begin()); + result.resize(16); + iterative_fft(boost::begin(a), boost::end(a), boost::begin(result)); + + gsl_fft_real_radix2_transform(a.data(), 1, a.size()); + + // TODO: Compare against benchmark result (GSL). +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_fft_agreement, T, signed_integral_test_types) +{ + // Tests whether the two implementations produce the same result. + std::vector a; + std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); + std::vector< std::complex > recursive_result; + std::vector< std::complex > iterative_result; + recursive_result.resize(16); + iterative_result.resize(16); + iterative_fft(boost::begin(a), boost::end(a), boost::begin(iterative_result)); + recursive_fft(a.begin(), a.end(), recursive_result.begin()); + BOOST_CHECK_EQUAL_COLLECTIONS(iterative_result.begin(), iterative_result.end(), recursive_result.begin(), recursive_result.end()); } From 57b0e6589477db02b007cf20be20d459783026b1 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 10 Apr 2016 00:30:10 +1000 Subject: [PATCH 04/12] Simplify bit-reversal incrementing. --- include/boost/math/tools/fft.hpp | 36 ++++++++++++------------------------ test/test_fft.cpp | 9 --------- 2 files changed, 12 insertions(+), 33 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index bcc01cf0b..175480aeb 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -114,7 +114,7 @@ template void increment_reversed_integer(N &x, N m) { x ^= m; - if (x <= m) + if (x < m) { do { @@ -126,37 +126,24 @@ void increment_reversed_integer(N &x, N m) } -// Presumably this can be done more elegantly (and efficiently) in combination -// with a permutation_iterator, etc. -template -void bit_reversed_permute(I first, I last) -{ - typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type N; - N const n = std::distance(first, last); - BOOST_ASSERT((n & (n - 1)) == 0); // This might not be a requirement. - - N const m = n / 2; - N i_rev = m; - for (I i = first + 1; i != last - 1; i++) - { - if (i_rev > i - first) - std::iter_swap(i, first + i_rev); - increment_reversed_integer(i_rev, m); - } -} - - // Iterative FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 template void iterative_fft(I first, I last, J A) { typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type complex_t; typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type N; - std::copy(first, last, A); + // Copy and permute the input. + // TODO: Presumably can be done more elegantly with permutation_iterator. N const n = std::distance(first, last); - bit_reversed_permute(A, A + n); + N i_rev = 0; + for (; first != last; first++) + { + A[i_rev] = *first; + increment_reversed_integer(i_rev, n / 2); + } + complex_t const i(0, 1); - N log2_n = log2(n); // TODO: Err... C99? + N const log2_n = log2(n); // TODO: Err... C99? for (N s = 1; s <= log2_n; s++) { N const m = N(1) << s; @@ -168,6 +155,7 @@ void iterative_fft(I first, I last, J A) for (N j = 0; j < m_2; j++) { complex_t const t = w * A[k + j + m_2]; + // TODO: Assignment into arrays probably not good for performance. A[k + j + m_2] = A[k + j] - t; A[k + j] += t; w *= w_m; diff --git a/test/test_fft.cpp b/test/test_fft.cpp index d74cf3eb8..7351f3755 100644 --- a/test/test_fft.cpp +++ b/test/test_fft.cpp @@ -41,15 +41,6 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_recursive_fft, T, signed_integral_test_types) } -BOOST_AUTO_TEST_CASE_TEMPLATE(test_bit_reversed_permute, T, signed_integral_test_types) -{ - std::vector a; - std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); - bit_reversed_permute(a.begin(), a.end()); - boost::array b = {0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15}; - BOOST_CHECK_EQUAL_COLLECTIONS(boost::begin(a), boost::end(a), boost::begin(b), boost::end(b)); -} - BOOST_AUTO_TEST_CASE_TEMPLATE(test_iterative_fft, T, signed_integral_test_types) { std::vector a; From 6ea11cfaf80a08096ffa54503e3cc128f159190e Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 10 Apr 2016 00:38:03 +1000 Subject: [PATCH 05/12] Make a note and a bugfix about the stride iterator. --- include/boost/math/tools/fft.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index 175480aeb..0c9d972b0 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -36,9 +36,13 @@ namespace detail stride_iterator(Iterator x, base_distance stride) : super_t(x), stride(stride) {} template - stride_iterator(stride_iterator const &x, typename enable_if_convertible::type* = 0) : super_t(x.base()) + stride_iterator(stride_iterator const &x, typename enable_if_convertible::type* = 0) : super_t(x.base()), stride(x.stride) {} + /* This constructor is the main reason for hand-rolling a stride + * iterator: when striding a stride iterator, I want to simply multiply + * the stride values, not make nested stride iterators. + */ template stride_iterator(stride_iterator const &x, base_distance stride, typename enable_if_convertible::type* = 0) : super_t(x.base()), stride(x.stride * stride) {} From 3041f51937071920a3ae45d2c65afdb7a301635e Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 11 Apr 2016 09:27:45 +1000 Subject: [PATCH 06/12] Simplify reversed integer increment. --- include/boost/math/tools/fft.hpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index 0c9d972b0..9528bc2de 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -118,14 +118,10 @@ template void increment_reversed_integer(N &x, N m) { x ^= m; - if (x < m) + while (x < m) { - do - { - m >>= 1; - x ^= m; - } - while (x < m); + m >>= 1; + x ^= m; } } From 667df2aa4cd76769caa0ba9f029156e306f4c88d Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 11 Apr 2016 09:28:52 +1000 Subject: [PATCH 07/12] fft: forward/inverse wrappers. --- include/boost/math/tools/fft.hpp | 47 ++++++++++++++++++++++++++++++++-------- test/test_fft.cpp | 23 +++++++++++++++++--- 2 files changed, 58 insertions(+), 12 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index 9528bc2de..1457f4128 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -11,7 +11,9 @@ #include #include #include +#include +#include #include #include #include @@ -76,6 +78,16 @@ namespace detail { return stride_iterator(x, stride); } + + + template + struct identity + { + T operator()(T const &x) const + { + return x; + } + }; } @@ -127,10 +139,10 @@ void increment_reversed_integer(N &x, N m) // Iterative FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 -template -void iterative_fft(I first, I last, J A) +template +void iterative_fft_radix2(I first, I last, J A, R1 sign) { - typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type complex_t; + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type output_type; typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type N; // Copy and permute the input. // TODO: Presumably can be done more elegantly with permutation_iterator. @@ -142,19 +154,19 @@ void iterative_fft(I first, I last, J A) increment_reversed_integer(i_rev, n / 2); } - complex_t const i(0, 1); + output_type const i(0, 1); N const log2_n = log2(n); // TODO: Err... C99? for (N s = 1; s <= log2_n; s++) { N const m = N(1) << s; - complex_t const w_m = std::exp((2.0 * M_PI * i) / static_cast(m)); + output_type const w_m = std::exp(sign(2.0 * M_PI * i) / static_cast(m)); for (N k = 0; k < n; k += m) { - complex_t w = 1; + output_type w = 1; N const m_2 = m / 2; for (N j = 0; j < m_2; j++) { - complex_t const t = w * A[k + j + m_2]; + output_type const t = w * A[k + j + m_2]; // TODO: Assignment into arrays probably not good for performance. A[k + j + m_2] = A[k + j] - t; A[k + j] += t; @@ -164,8 +176,25 @@ void iterative_fft(I first, I last, J A) } } -}}} -#endif +template +void iterative_fft_radix2_forward(I first, I last, J A) +{ + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type complex_t; + iterative_fft_radix2(first, last, A, std::negate()); +} +template +void iterative_fft_radix2_inverse(I first, I last, J A) +{ + using namespace boost::lambda; + typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type complex_t; + + iterative_fft_radix2(first, last, A, detail::identity()); + std::transform(A, A + (last - first), A, _1 / static_cast(last - first)); +} + +}}} + +#endif diff --git a/test/test_fft.cpp b/test/test_fft.cpp index 7351f3755..56c9be2f9 100644 --- a/test/test_fft.cpp +++ b/test/test_fft.cpp @@ -19,12 +19,14 @@ #include #include #include +#include typedef boost::mpl::list signed_integral_test_types; // , boost::multiprecision::cpp_int using namespace boost; using namespace boost::math::tools; +using boost::math::tools::detail::identity; typedef std::complex complex_t; @@ -41,13 +43,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_recursive_fft, T, signed_integral_test_types) } -BOOST_AUTO_TEST_CASE_TEMPLATE(test_iterative_fft, T, signed_integral_test_types) +BOOST_AUTO_TEST_CASE_TEMPLATE(test_iterative_fft_radix2, T, signed_integral_test_types) { std::vector a; std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); std::vector< std::complex > result; result.resize(16); - iterative_fft(boost::begin(a), boost::end(a), boost::begin(result)); + iterative_fft_radix2(boost::begin(a), boost::end(a), boost::begin(result), identity()); gsl_fft_real_radix2_transform(a.data(), 1, a.size()); @@ -64,7 +66,22 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_fft_agreement, T, signed_integral_test_types) std::vector< std::complex > iterative_result; recursive_result.resize(16); iterative_result.resize(16); - iterative_fft(boost::begin(a), boost::end(a), boost::begin(iterative_result)); + iterative_fft_radix2(boost::begin(a), boost::end(a), boost::begin(iterative_result), identity()); recursive_fft(a.begin(), a.end(), recursive_result.begin()); BOOST_CHECK_EQUAL_COLLECTIONS(iterative_result.begin(), iterative_result.end(), recursive_result.begin(), recursive_result.end()); } + + +BOOST_AUTO_TEST_CASE_TEMPLATE(test_fft_involution, T, signed_integral_test_types) +{ + std::vector a; + std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); + std::vector< std::complex > tmp; + tmp.resize(16); + iterative_fft_radix2_forward(a.begin(), a.end(), tmp.begin()); + std::vector b; + b.resize(16); + iterative_fft_radix2_inverse(tmp.begin(), tmp.end(), b.begin()); + for (std::size_t i = 0; i < a.size(); i++) + BOOST_CHECK_EQUAL(a[0], b[0].real()); +} From ccc9afcd25e2335c564b7fc23cb98c4da34d5ed9 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 11 Apr 2016 11:36:04 +1000 Subject: [PATCH 08/12] Add radix2 suffix to recursive_fft name. --- include/boost/math/tools/fft.hpp | 6 +++--- test/test_fft.cpp | 9 ++++++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index 1457f4128..4bb9ac497 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -93,7 +93,7 @@ namespace detail // Recursive FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 template -void recursive_fft(I first, I last, J y) +void recursive_fft_radix2(I first, I last, J y) { using namespace boost::math::tools::detail; @@ -113,8 +113,8 @@ void recursive_fft(I first, I last, J y) y_0.resize(n / 2); y_1.resize(n / 2); - recursive_fft(make_stride_iterator(first, 2), make_stride_iterator(last, 2), y_0.begin()); - recursive_fft(make_stride_iterator(first + 1, 2), make_stride_iterator(last + 1, 2), y_1.begin()); + recursive_fft_radix2(make_stride_iterator(first, 2), make_stride_iterator(last, 2), y_0.begin()); + recursive_fft_radix2(make_stride_iterator(first + 1, 2), make_stride_iterator(last + 1, 2), y_1.begin()); for (N k = 0; k < n / 2; k++) { diff --git a/test/test_fft.cpp b/test/test_fft.cpp index 56c9be2f9..bb81cda90 100644 --- a/test/test_fft.cpp +++ b/test/test_fft.cpp @@ -32,13 +32,13 @@ typedef std::complex complex_t; // TODO: Fixtures. -BOOST_AUTO_TEST_CASE_TEMPLATE(test_recursive_fft, T, signed_integral_test_types) +BOOST_AUTO_TEST_CASE_TEMPLATE(test_recursive_fft_radix2, T, signed_integral_test_types) { std::vector a; std::copy(make_counting_iterator(0), make_counting_iterator(16), std::back_inserter(a)); std::vector< std::complex > result; result.resize(16); - recursive_fft(a.begin(), a.end(), result.begin()); + recursive_fft_radix2(a.begin(), a.end(), result.begin()); // TODO: Compare against benchmark result (GSL). } @@ -67,11 +67,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_fft_agreement, T, signed_integral_test_types) recursive_result.resize(16); iterative_result.resize(16); iterative_fft_radix2(boost::begin(a), boost::end(a), boost::begin(iterative_result), identity()); - recursive_fft(a.begin(), a.end(), recursive_result.begin()); + recursive_fft_radix2(a.begin(), a.end(), recursive_result.begin()); BOOST_CHECK_EQUAL_COLLECTIONS(iterative_result.begin(), iterative_result.end(), recursive_result.begin(), recursive_result.end()); } +/* Check that fft_inverse(fft_forward(x)) == x. + * Might fail due to numerical error. + */ BOOST_AUTO_TEST_CASE_TEMPLATE(test_fft_involution, T, signed_integral_test_types) { std::vector a; From 029122bb46433a35243b35376101f921ec217268 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 14 Apr 2016 09:24:36 +1000 Subject: [PATCH 09/12] Minor loop variable simplification. --- include/boost/math/tools/fft.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index 4bb9ac497..a96958ba2 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -147,8 +147,7 @@ void iterative_fft_radix2(I first, I last, J A, R1 sign) // Copy and permute the input. // TODO: Presumably can be done more elegantly with permutation_iterator. N const n = std::distance(first, last); - N i_rev = 0; - for (; first != last; first++) + for (N i_rev = 0; first != last; first++) { A[i_rev] = *first; increment_reversed_integer(i_rev, n / 2); From 42ce81c7273e206a54902f3ec8b25b839881c59b Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 18 Apr 2016 09:21:08 +1000 Subject: [PATCH 10/12] [fft] Store w * y_1[k] for reuse in recursive algorithm. --- include/boost/math/tools/fft.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index a96958ba2..cfc8d5c4b 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -118,8 +118,9 @@ void recursive_fft_radix2(I first, I last, J y) for (N k = 0; k < n / 2; k++) { - y[k] = y_0[k] + w * y_1[k]; - y[k + n / 2] = y_0[k] - w * y_1[k]; + complex_t const t = w * y_1[k]; + y[k] = y_0[k] + t; + y[k + n / 2] = y_0[k] - t; w *= w_n; } } From 6c87601640abb0daf664c2d89b7210b8ed54c2c2 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Tue, 24 May 2016 09:15:48 +1000 Subject: [PATCH 11/12] Document parameters. --- include/boost/math/tools/fft.hpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index cfc8d5c4b..1e169301c 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -126,7 +126,11 @@ void recursive_fft_radix2(I first, I last, J y) } -// Based on code from Henry Warren Jr, Hacker's Delight, 2nd ed., 2013, p 138 +/** + * Based on code from Henry Warren Jr, Hacker's Delight, 2nd ed., 2013, p 138 + * @param x Reversed integer to increment. + * @param m Increment. + */ template void increment_reversed_integer(N &x, N m) { @@ -139,7 +143,14 @@ void increment_reversed_integer(N &x, N m) } -// Iterative FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 +/** + * Iterative FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 + * @tparam I InputIterator + * @tparam J Mutable RandomAccessIterator + * @tparam R1 Unary Transformation f(T) -> T + * @param first First element of data + * @param last One after last element of data + */ template void iterative_fft_radix2(I first, I last, J A, R1 sign) { From 8ea3d101f478d1563a17f7aac54a08ebab711d10 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 23 May 2016 20:01:35 +1000 Subject: [PATCH 12/12] Remove recursive implementation. It's probably not the implementation we'll actually use. --- include/boost/math/tools/fft.hpp | 105 --------------------------------------- 1 file changed, 105 deletions(-) diff --git a/include/boost/math/tools/fft.hpp b/include/boost/math/tools/fft.hpp index 1e169301c..ebe8dd034 100644 --- a/include/boost/math/tools/fft.hpp +++ b/include/boost/math/tools/fft.hpp @@ -20,111 +20,6 @@ #include namespace boost { namespace math { namespace tools { - -namespace detail -{ - template - class stride_iterator : public iterator_adaptor, Iterator> - { - typedef iterator_adaptor, Iterator> super_t; - typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type base_distance; - friend class iterator_core_access; - - base_distance stride; - - public: - stride_iterator() {} - - stride_iterator(Iterator x, base_distance stride) : super_t(x), stride(stride) {} - - template - stride_iterator(stride_iterator const &x, typename enable_if_convertible::type* = 0) : super_t(x.base()), stride(x.stride) - {} - - /* This constructor is the main reason for hand-rolling a stride - * iterator: when striding a stride iterator, I want to simply multiply - * the stride values, not make nested stride iterators. - */ - template - stride_iterator(stride_iterator const &x, base_distance stride, typename enable_if_convertible::type* = 0) : super_t(x.base()), stride(x.stride * stride) - {} - - // TODO: These *should* be private but compilation fails unexpectedly. - template - BOOST_DEDUCED_TYPENAME super_t::difference_type - distance_to(iterator_adaptor const& y) const - { - return (y.base() - this->base_reference()) / stride; - } - - void advance(BOOST_DEDUCED_TYPENAME super_t::difference_type n) - { - std::advance(this->base_reference(), stride * n); - } - private: - void increment() { std::advance(this->base_reference(), stride); } - void decrement() { std::advance(this->base_reference(), -stride); } - - }; - - template - stride_iterator make_stride_iterator(stride_iterator x, N stride) - { - return stride_iterator(x, stride); - } - - template - stride_iterator make_stride_iterator(Iterator x, N stride) - { - return stride_iterator(x, stride); - } - - - template - struct identity - { - T operator()(T const &x) const - { - return x; - } - }; -} - - -// Recursive FFT based on Cormen et al, Introduction to algorithms, 3rd ed., 2009 -template -void recursive_fft_radix2(I first, I last, J y) -{ - using namespace boost::math::tools::detail; - - typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::difference_type N; - typedef BOOST_DEDUCED_TYPENAME std::iterator_traits::value_type complex_t; - N const n = distance(first, last); - BOOST_ASSERT((n & (n - 1)) == 0); // We require that n is a power of 2. - if (n == 1) - { - *y = *first; - return; - } - complex_t const i(0, 1); - complex_t const w_n = std::exp((2.0 * M_PI * i) / static_cast(n)); - complex_t w = 1; - std::vector y_0, y_1; - y_0.resize(n / 2); - y_1.resize(n / 2); - - recursive_fft_radix2(make_stride_iterator(first, 2), make_stride_iterator(last, 2), y_0.begin()); - recursive_fft_radix2(make_stride_iterator(first + 1, 2), make_stride_iterator(last + 1, 2), y_1.begin()); - - for (N k = 0; k < n / 2; k++) - { - complex_t const t = w * y_1[k]; - y[k] = y_0[k] + t; - y[k + n / 2] = y_0[k] - t; - w *= w_n; - } -} - /** * Based on code from Henry Warren Jr, Hacker's Delight, 2nd ed., 2013, p 138