|
5 | 5 |
|
6 | 6 | #include <xtl/xcomplex.hpp> |
7 | 7 |
|
8 | | -#include "./xtl_concepts.hpp" |
9 | 8 | #include "../containers/xarray.hpp" |
10 | 9 | #include "../core/xmath.hpp" |
11 | 10 | #include "../core/xnoalias.hpp" |
12 | 11 | #include "../generators/xbuilder.hpp" |
13 | 12 | #include "../misc/xcomplex.hpp" |
14 | 13 | #include "../views/xaxis_slice_iterator.hpp" |
15 | 14 | #include "../views/xview.hpp" |
| 15 | +#include "./xtl_concepts.hpp" |
16 | 16 |
|
17 | 17 | namespace xt |
18 | 18 | { |
@@ -118,65 +118,65 @@ namespace xt |
118 | 118 | } |
119 | 119 | } // namespace detail |
120 | 120 |
|
121 | | - /** |
| 121 | + /** |
122 | 122 | * @brief 1D FFT of an Nd array along a specified axis |
123 | 123 | * @param e an Nd expression to be transformed to the fourier domain |
124 | 124 | * @param axis the axis along which to perform the 1D FFT |
125 | 125 | * @return a transformed xarray of the specified precision |
126 | 126 | */ |
127 | | - template<class E> |
| 127 | + template <class E> |
128 | 128 | inline auto fft(E&& e, std::ptrdiff_t axis = -1) |
129 | 129 | { |
130 | | - using value_type = typename std::decay<E>::type::value_type; |
131 | | - if constexpr (xtl::is_complex<typename std::decay<E>::type::value_type>::value) |
132 | | - { |
133 | | - using precision = typename value_type::value_type; |
134 | | - const auto saxis = xt::normalize_axis(e.dimension(), axis); |
135 | | - const size_t N = e.shape(saxis); |
136 | | - const bool powerOfTwo = !(N == 0) && !(N & (N - 1)); |
137 | | - xt::xarray<std::complex<precision>> out = xt::eval(e); |
138 | | - auto begin = xt::axis_slice_begin(out, saxis); |
139 | | - auto end = xt::axis_slice_end(out, saxis); |
140 | | - for (auto iter = begin; iter != end; iter++) |
| 130 | + using value_type = typename std::decay<E>::type::value_type; |
| 131 | + if constexpr (xtl::is_complex<typename std::decay<E>::type::value_type>::value) |
141 | 132 | { |
142 | | - if (powerOfTwo) |
143 | | - { |
144 | | - xt::noalias(*iter) = detail::radix2(*iter); |
145 | | - } |
146 | | - else |
| 133 | + using precision = typename value_type::value_type; |
| 134 | + const auto saxis = xt::normalize_axis(e.dimension(), axis); |
| 135 | + const size_t N = e.shape(saxis); |
| 136 | + const bool powerOfTwo = !(N == 0) && !(N & (N - 1)); |
| 137 | + xt::xarray<std::complex<precision>> out = xt::eval(e); |
| 138 | + auto begin = xt::axis_slice_begin(out, saxis); |
| 139 | + auto end = xt::axis_slice_end(out, saxis); |
| 140 | + for (auto iter = begin; iter != end; iter++) |
147 | 141 | { |
148 | | - xt::noalias(*iter) = detail::transform_bluestein(*iter); |
| 142 | + if (powerOfTwo) |
| 143 | + { |
| 144 | + xt::noalias(*iter) = detail::radix2(*iter); |
| 145 | + } |
| 146 | + else |
| 147 | + { |
| 148 | + xt::noalias(*iter) = detail::transform_bluestein(*iter); |
| 149 | + } |
149 | 150 | } |
| 151 | + return out; |
150 | 152 | } |
151 | | - return out; |
152 | | - } |
153 | | - else |
154 | | - { |
155 | | - return fft(xt::cast<std::complex<value_type>>(e), axis); |
156 | | - } |
157 | | - } |
| 153 | + else |
| 154 | + { |
| 155 | + return fft(xt::cast<std::complex<value_type>>(e), axis); |
| 156 | + } |
| 157 | + } |
158 | 158 |
|
159 | 159 | template <class E> |
160 | 160 | inline auto ifft(E&& e, std::ptrdiff_t axis = -1) |
161 | 161 | { |
162 | | - if constexpr (xtl::is_complex<typename std::decay<E>::type::value_type>::value) |
163 | | - { |
164 | | - // check the length of the data on that axis |
165 | | - const std::size_t n = e.shape(axis); |
166 | | - if (n == 0) |
| 162 | + if constexpr (xtl::is_complex<typename std::decay<E>::type::value_type>::value) |
| 163 | + { |
| 164 | + // check the length of the data on that axis |
| 165 | + const std::size_t n = e.shape(axis); |
| 166 | + if (n == 0) |
| 167 | + { |
| 168 | + XTENSOR_THROW(std::runtime_error, "Cannot take the iFFT along an empty dimention"); |
| 169 | + } |
| 170 | + auto complex_args = xt::conj(e); |
| 171 | + auto fft_res = xt::fft::fft(complex_args, axis); |
| 172 | + fft_res = xt::conj(fft_res); |
| 173 | + return fft_res; |
| 174 | + } |
| 175 | + else |
167 | 176 | { |
168 | | - XTENSOR_THROW(std::runtime_error, "Cannot take the iFFT along an empty dimention"); |
| 177 | + using value_type = typename std::decay<E>::type::value_type; |
| 178 | + return ifft(xt::cast<std::complex<value_type>>(e), axis); |
169 | 179 | } |
170 | | - auto complex_args = xt::conj(e); |
171 | | - auto fft_res = xt::fft::fft(complex_args, axis); |
172 | | - fft_res = xt::conj(fft_res); |
173 | | - return fft_res; |
174 | | - } |
175 | | - else |
176 | | - { |
177 | | - using value_type = typename std::decay<E>::type::value_type; |
178 | | - return ifft(xt::cast<std::complex<value_type>>(e), axis); |
179 | | - } |
180 | 180 | } |
181 | 181 |
|
182 | 182 | /* |
|
0 commit comments