From d1fc657c1fe7b4ad70452cb37b13a377fa831e56 Mon Sep 17 00:00:00 2001 From: bolero-MURAKAMI Date: Fri, 18 Apr 2014 19:33:48 +0900 Subject: [PATCH] fix fft implementation --- sprout/complex.hpp | 1 + sprout/complex/complex.hpp | 6 ++--- sprout/complex/type_traits.hpp | 42 ++++++++++++++++++++++++++++++++++ sprout/complex/values.hpp | 4 ++-- sprout/numeric/fft/fft.hpp | 36 +++++++++++++++++++---------- 5 files changed, 72 insertions(+), 17 deletions(-) create mode 100644 sprout/complex/type_traits.hpp diff --git a/sprout/complex.hpp b/sprout/complex.hpp index ebf8decf..2bc5f937 100644 --- a/sprout/complex.hpp +++ b/sprout/complex.hpp @@ -16,5 +16,6 @@ #include #include #include +#include #endif // #ifndef SPROUT_COMPLEX_HPP diff --git a/sprout/complex/complex.hpp b/sprout/complex/complex.hpp index 0edd5509..3c05d066 100644 --- a/sprout/complex/complex.hpp +++ b/sprout/complex/complex.hpp @@ -15,9 +15,9 @@ namespace sprout { class complex; namespace detail { - template - inline SPROUT_CONSTEXPR T - complex_norm(sprout::complex const& x) { + template + inline SPROUT_CONSTEXPR typename Complex::value_type + complex_norm(Complex const& x) { return x.real() * x.real() + x.imag() * x.imag(); } } // namespace detail diff --git a/sprout/complex/type_traits.hpp b/sprout/complex/type_traits.hpp new file mode 100644 index 00000000..0ad74b3e --- /dev/null +++ b/sprout/complex/type_traits.hpp @@ -0,0 +1,42 @@ +/*============================================================================= + 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_COMPLEX_TYPE_TRAITS_HPP +#define SPROUT_COMPLEX_TYPE_TRAITS_HPP + +#include +#include +#include + +namespace sprout { + // + // is_complex + // + template + struct is_complex + : public sprout::false_type + {}; + template + struct is_complex + : public sprout::is_complex + {}; + template + struct is_complex + : public sprout::is_complex + {}; + template + struct is_complex > + : public sprout::true_type + {}; + +#if SPROUT_USE_VARIABLE_TEMPLATES + template + SPROUT_STATIC_CONSTEXPR bool is_complex_v = sprout::is_complex::value; +#endif // #if SPROUT_USE_VARIABLE_TEMPLATES +} // namespace sprout + +#endif // #ifndef SPROUT_COMPLEX_TYPE_TRAITS_HPP diff --git a/sprout/complex/values.hpp b/sprout/complex/values.hpp index 7a41dc5c..2d8aaec2 100644 --- a/sprout/complex/values.hpp +++ b/sprout/complex/values.hpp @@ -33,12 +33,12 @@ namespace sprout { return x.imag(); } template - inline SPROUT_CXX14_CONSTEXPR T& + inline SPROUT_CONSTEXPR T& real(sprout::complex& x) { return x.real(); } template - inline SPROUT_CXX14_CONSTEXPR T& + inline SPROUT_CONSTEXPR T& imag(sprout::complex& x) { return x.imag(); } diff --git a/sprout/numeric/fft/fft.hpp b/sprout/numeric/fft/fft.hpp index e4224dc4..3448d757 100644 --- a/sprout/numeric/fft/fft.hpp +++ b/sprout/numeric/fft/fft.hpp @@ -51,10 +51,14 @@ namespace sprout { 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); + first[j1] = value_type( + real(x1) - imag(x3), + imag(x1) + real(x3) + ); + first[j3] = value_type( + real(x1) + imag(x3), + imag(x1) - real(x3) + ); } } if (m == size) { @@ -73,12 +77,16 @@ namespace sprout { 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); + first[j1] = value_type( + w1r * (x0r + x0i), + 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); + first[j3] = value_type( + w1r * (-x0r + x0i), + w1r * (-x0i - x0r) + ); } } for (difference_type i = 2 * m; i < size; i += m) { @@ -99,12 +107,16 @@ namespace sprout { 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; + first[j1] = value_type( + w1r * x0r - w1i * x0i, + 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; + first[j3] = value_type( + w3r * x0r - w3i * x0i, + w3r * x0i + w3i * x0r + ); } } }