GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-01 16:24:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

Disable use of recurrence for computing twiddle factors.

Browse Source
This commit is contained in:
Benoit Steiner 2018-01-09 18:32:52 +00:00
commit 22de74aa76
2 changed files with 54 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
unsupported/Eigen/CXX11/src/Tensor/TensorFFT.h
View File

@ -231,20 +231,32 @@ struct TensorEvaluator<const TensorFFTOp<FFT, ArgType, FFTResultType, FFTDir>, D
// t_n = exp(sqrt(-1) * pi * n^2 / line_len) // t_n = exp(sqrt(-1) * pi * n^2 / line_len)
// for n = 0, 1,..., line_len-1. // for n = 0, 1,..., line_len-1.
// For n > 2 we use the recurrence t_n = t_{n-1}^2 / t_{n-2} * t_1^2 // For n > 2 we use the recurrence t_n = t_{n-1}^2 / t_{n-2} * t_1^2
pos_j_base_powered[0] = ComplexScalar(1, 0);
if (line_len > 1) { // The recurrence is correct in exact arithmetic, but causes
const RealScalar pi_over_len(EIGEN_PI / line_len); // numerical issues for large transforms, especially in
const ComplexScalar pos_j_base = ComplexScalar( // single-precision floating point.
std::cos(pi_over_len), std::sin(pi_over_len)); //
pos_j_base_powered[1] = pos_j_base; // pos_j_base_powered[0] = ComplexScalar(1, 0);
if (line_len > 2) { // if (line_len > 1) {
const ComplexScalar pos_j_base_sq = pos_j_base * pos_j_base; // const ComplexScalar pos_j_base = ComplexScalar(
for (int j = 2; j < line_len + 1; ++j) { // numext::cos(M_PI / line_len), numext::sin(M_PI / line_len));
pos_j_base_powered[j] = pos_j_base_powered[j - 1] * // pos_j_base_powered[1] = pos_j_base;
pos_j_base_powered[j - 1] / // if (line_len > 2) {
pos_j_base_powered[j - 2] * pos_j_base_sq; // const ComplexScalar pos_j_base_sq = pos_j_base * pos_j_base;
} // for (int i = 2; i < line_len + 1; ++i) {
} // pos_j_base_powered[i] = pos_j_base_powered[i - 1] *
// pos_j_base_powered[i - 1] /
// pos_j_base_powered[i - 2] *
// pos_j_base_sq;
// }
// }
// }
// TODO(rmlarsen): Find a way to use Eigen's vectorized sin
// and cosine functions here.
for (int j = 0; j < line_len + 1; ++j) {
double arg = ((EIGEN_PI * j) * j) / line_len;
std::complex<double> tmp(numext::cos(arg), numext::sin(arg));
pos_j_base_powered[j] = static_cast<ComplexScalar>(tmp);
} }
} }

28
unsupported/test/cxx11_tensor_fft.cpp
View File

@ -224,6 +224,32 @@ static void test_fft_real_input_energy() {
} }
} }
template <typename RealScalar>
static void test_fft_non_power_of_2_round_trip(int exponent) {
int n = (1 << exponent) + 1;
Eigen::DSizes<long, 1> dimensions;
dimensions[0] = n;
const DSizes<long, 1> arr = dimensions;
Tensor<RealScalar, 1, ColMajor, long> input;
input.resize(arr);
input.setRandom();
array<int, 1> fft;
fft[0] = 0;
Tensor<std::complex<RealScalar>, 1, ColMajor> forward =
input.template fft<BothParts, FFT_FORWARD>(fft);
Tensor<RealScalar, 1, ColMajor, long> output =
forward.template fft<RealPart, FFT_REVERSE>(fft);
for (int i = 0; i < n; ++i) {
VERIFY_IS_APPROX(input[i], output[i]);
}
}
void test_cxx11_tensor_fft() { void test_cxx11_tensor_fft() {
test_fft_complex_input_golden(); test_fft_complex_input_golden();
test_fft_real_input_golden(); test_fft_real_input_golden();
@ -270,4 +296,6 @@ void test_cxx11_tensor_fft() {
test_fft_real_input_energy<RowMajor, double, true, Eigen::BothParts, FFT_FORWARD, 4>(); test_fft_real_input_energy<RowMajor, double, true, Eigen::BothParts, FFT_FORWARD, 4>();
test_fft_real_input_energy<RowMajor, float, false, Eigen::BothParts, FFT_FORWARD, 4>(); test_fft_real_input_energy<RowMajor, float, false, Eigen::BothParts, FFT_FORWARD, 4>();
test_fft_real_input_energy<RowMajor, double, false, Eigen::BothParts, FFT_FORWARD, 4>(); test_fft_real_input_energy<RowMajor, double, false, Eigen::BothParts, FFT_FORWARD, 4>();
test_fft_non_power_of_2_round_trip<float>(7);
} }