From 07b5f69ebbdf8c00d9dd9034a7f24195508ed2d4 Mon Sep 17 00:00:00 2001 From: bolero-MURAKAMI Date: Fri, 18 Apr 2014 11:34:28 +0900 Subject: [PATCH] add C++14 constexpr fft/ifft --- sprout/complex/complex.hpp | 10 ++- sprout/complex/values.hpp | 14 +++- sprout/numeric/fft.hpp | 2 + sprout/numeric/fft/fft.hpp | 125 ++++++++++++++++++++++++++++++++++++ sprout/numeric/fft/ifft.hpp | 58 +++++++++++++++++ 5 files changed, 205 insertions(+), 4 deletions(-) create mode 100644 sprout/numeric/fft/fft.hpp create mode 100644 sprout/numeric/fft/ifft.hpp diff --git a/sprout/complex/complex.hpp b/sprout/complex/complex.hpp index 627a6929..0edd5509 100644 --- a/sprout/complex/complex.hpp +++ b/sprout/complex/complex.hpp @@ -41,13 +41,19 @@ namespace sprout { SPROUT_CONSTEXPR complex(complex const& other) SPROUT_NOEXCEPT : re_(other.real()), im_(other.imag()) {} - SPROUT_CONSTEXPR T real() const SPROUT_NOEXCEPT { + SPROUT_CONSTEXPR T const& real() const SPROUT_NOEXCEPT { + return re_; + } + SPROUT_CXX14_CONSTEXPR T& real() SPROUT_NOEXCEPT { return re_; } SPROUT_CXX14_CONSTEXPR void real(T re) SPROUT_NOEXCEPT { re_ = re; } - SPROUT_CONSTEXPR T imag() const SPROUT_NOEXCEPT { + SPROUT_CONSTEXPR T const& imag() const SPROUT_NOEXCEPT { + return im_; + } + SPROUT_CXX14_CONSTEXPR T& imag() SPROUT_NOEXCEPT { return im_; } SPROUT_CXX14_CONSTEXPR void imag(T im) SPROUT_NOEXCEPT { diff --git a/sprout/complex/values.hpp b/sprout/complex/values.hpp index 27e91298..7a41dc5c 100644 --- a/sprout/complex/values.hpp +++ b/sprout/complex/values.hpp @@ -23,16 +23,26 @@ namespace sprout { // 26.4.7, values: template - inline SPROUT_CONSTEXPR T + inline SPROUT_CONSTEXPR T const& real(sprout::complex const& x) { return x.real(); } template - inline SPROUT_CONSTEXPR T + inline SPROUT_CONSTEXPR T const& imag(sprout::complex const& x) { return x.imag(); } template + inline SPROUT_CXX14_CONSTEXPR T& + real(sprout::complex& x) { + return x.real(); + } + template + inline SPROUT_CXX14_CONSTEXPR T& + imag(sprout::complex& x) { + return x.imag(); + } + template inline SPROUT_CONSTEXPR T abs(sprout::complex const& x) { return sprout::sqrt(sprout::norm(x)); diff --git a/sprout/numeric/fft.hpp b/sprout/numeric/fft.hpp index 2ff8a461..93b6d481 100644 --- a/sprout/numeric/fft.hpp +++ b/sprout/numeric/fft.hpp @@ -9,6 +9,8 @@ #define SPROUT_NUMERIC_FFT_HPP #include +#include +#include #include #endif // #ifndef SPROUT_NUMERIC_FFT_HPP diff --git a/sprout/numeric/fft/fft.hpp b/sprout/numeric/fft/fft.hpp new file mode 100644 index 00000000..e4224dc4 --- /dev/null +++ b/sprout/numeric/fft/fft.hpp @@ -0,0 +1,125 @@ +/*============================================================================= + Copyright (c) 2011-2014 Bolero MURAKAMI + https://github.com/bolero-MURAKAMI/Sprout + + Distributed under 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 SPROUT_NUMERIC_FFT_FFT_HPP +#define SPROUT_NUMERIC_FFT_FFT_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace sprout { + // + // fft + // + template + inline SPROUT_CXX14_CONSTEXPR void + fft(RandomAccessIterator first, RandomAccessIterator last) { + typedef typename std::iterator_traits::difference_type difference_type; + typedef typename std::iterator_traits::value_type value_type; + typedef typename value_type::value_type elem_type; + using sprout::real; + using sprout::imag; + difference_type const size = last - first; + // scrambler + for (difference_type i = 0, j = 1; j < size - 1; ++j) { + for (difference_type k = size >> 1; k > (i ^= k); k >>= 1) + ; + if (j < i) { + sprout::swap(first[i], first[j]); + } + } + // L shaped butterflies + elem_type const theta = -(sprout::math::half_pi() / size); + for (difference_type m = 4; m <= size; m <<= 1) { + difference_type mq = m >> 2; + // W == 1 + for (difference_type k = mq; k >= 1; k >>= 2) { + for (difference_type j = mq - k; j < mq - (k >> 1); ++j) { + difference_type j1 = j + mq; + difference_type j2 = j1 + mq; + difference_type j3 = j2 + mq; + value_type x1 = first[j] - first[j1]; + first[j] += first[j1]; + value_type x3 = first[j3] - first[j2]; + first[j2] += first[j3]; + real(first[j1]) = real(x1) - imag(x3); + imag(first[j1]) = imag(x1) + real(x3); + real(first[j3]) = real(x1) + imag(x3); + imag(first[j3]) = imag(x1) - real(x3); + } + } + if (m == size) { + continue; + } + difference_type irev = size >> 1; + elem_type w1r = sprout::cos(theta * irev); + for (difference_type k = mq; k >= 1; k >>= 2) { + for (difference_type j = m + mq - k; j < m + mq - (k >> 1); ++j) { + difference_type j1 = j + mq; + difference_type j2 = j1 + mq; + difference_type j3 = j2 + mq; + value_type x1 = first[j] - first[j1]; + first[j] += first[j1]; + value_type x3 = first[j3] - first[j2]; + first[j2] += first[j3]; + elem_type x0r = real(x1) - imag(x3); + elem_type x0i = imag(x1) + real(x3); + real(first[j1]) = w1r * (x0r + x0i); + imag(first[j1]) = w1r * (x0i - x0r); + x0r = real(x1) + imag(x3); + x0i = imag(x1) - real(x3); + real(first[j3]) = w1r * (-x0r + x0i); + imag(first[j3]) = w1r * (-x0i - x0r); + } + } + for (difference_type i = 2 * m; i < size; i += m) { + for (difference_type k = size >> 1; k > (irev ^= k); k >>= 1) + ; + elem_type w1r = sprout::cos(theta * irev); + elem_type w1i = sprout::sin(theta * irev); + elem_type w3r = sprout::cos(theta * 3 * irev); + elem_type w3i = sprout::sin(theta * 3 * irev); + for (difference_type k = mq; k >= 1; k >>= 2) { + for (difference_type j = i + mq - k; j < i + mq - (k >> 1); ++j) { + difference_type j1 = j + mq; + difference_type j2 = j1 + mq; + difference_type j3 = j2 + mq; + value_type x1 = first[j] - first[j1]; + first[j] += first[j1]; + value_type x3 = first[j3] - first[j2]; + first[j2] += first[j3]; + elem_type x0r = real(x1) - imag(x3); + elem_type x0i = imag(x1) + real(x3); + real(first[j1]) = w1r * x0r - w1i * x0i; + imag(first[j1]) = w1r * x0i + w1i * x0r; + x0r = real(x1) + imag(x3); + x0i = imag(x1) - real(x3); + real(first[j3]) = w3r * x0r - w3i * x0i; + imag(first[j3]) = w3r * x0i + w3i * x0r; + } + } + } + } + // radix 2 butterflies + difference_type mq = size >> 1; + for (difference_type k = mq; k >= 1; k >>= 2) { + for (difference_type j = mq - k; j < mq - (k >> 1); ++j) { + difference_type j1 = mq + j; + value_type x0 = first[j] - first[j1]; + first[j] += first[j1]; + first[j1] = x0; + } + } + } +} // namespace sprout + +#endif // #ifndef SPROUT_NUMERIC_FFT_FFT_HPP diff --git a/sprout/numeric/fft/ifft.hpp b/sprout/numeric/fft/ifft.hpp new file mode 100644 index 00000000..89c7be96 --- /dev/null +++ b/sprout/numeric/fft/ifft.hpp @@ -0,0 +1,58 @@ +/*============================================================================= + Copyright (c) 2011-2014 Bolero MURAKAMI + https://github.com/bolero-MURAKAMI/Sprout + + Distributed under 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 SPROUT_NUMERIC_FFT_IFFT_HPP +#define SPROUT_NUMERIC_FFT_IFFT_HPP + +#include +#include +#include +#include +#include + +namespace sprout { + namespace detail { + struct conj_value { + public: + template + SPROUT_CONSTEXPR Complex operator()(Complex const& x) const { + return conj(x); + } + }; + + template + struct conj_div_value { + private: + Size n_; + public: + explicit SPROUT_CONSTEXPR conj_div_value(Size n) + : n_(n) + {} + template + SPROUT_CONSTEXPR Complex operator()(Complex const& x) const { + return conj(x) / typename Complex::value_type(n_); + } + }; + template + SPROUT_CONSTEXPR sprout::detail::conj_div_value + conj_div(Size n) { + return sprout::detail::conj_div_value(n); + } + } // namespace detail + // + // ifft + // + template + inline SPROUT_CXX14_CONSTEXPR void + ifft(RandomAccessIterator first, RandomAccessIterator last) { + sprout::transform(first, last, first, sprout::detail::conj_value()); + sprout::fft(first, last); + sprout::transform(first, last, first, sprout::detail::conj_div(last - first)); + } +} // namespace sprout + +#endif // #ifndef SPROUT_NUMERIC_FFT_IFFT_HPP