cxx11_tensor_fft.cpp
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2014 Jianwei Cui <thucjw@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #include "main.h"
11 #include <Eigen/CXX11/Tensor>
12 
13 using Eigen::Tensor;
14 
15 template <int DataLayout>
16 static void test_fft_2D_golden() {
17  Tensor<float, 2, DataLayout> input(2, 3);
18  input(0, 0) = 1;
19  input(0, 1) = 2;
20  input(0, 2) = 3;
21  input(1, 0) = 4;
22  input(1, 1) = 5;
23  input(1, 2) = 6;
24 
25  array<ptrdiff_t, 2> fft;
26  fft[0] = 0;
27  fft[1] = 1;
28 
29  Tensor<std::complex<float>, 2, DataLayout> output = input.template fft<Eigen::BothParts, Eigen::FFT_FORWARD>(fft);
30 
31  std::complex<float> output_golden[6]; // in ColMajor order
32  output_golden[0] = std::complex<float>(21, 0);
33  output_golden[1] = std::complex<float>(-9, 0);
34  output_golden[2] = std::complex<float>(-3, 1.73205);
35  output_golden[3] = std::complex<float>( 0, 0);
36  output_golden[4] = std::complex<float>(-3, -1.73205);
37  output_golden[5] = std::complex<float>(0 ,0);
38 
39  std::complex<float> c_offset = std::complex<float>(1.0, 1.0);
40 
41  if (DataLayout == ColMajor) {
42  VERIFY_IS_APPROX(output(0) + c_offset, output_golden[0] + c_offset);
43  VERIFY_IS_APPROX(output(1) + c_offset, output_golden[1] + c_offset);
44  VERIFY_IS_APPROX(output(2) + c_offset, output_golden[2] + c_offset);
45  VERIFY_IS_APPROX(output(3) + c_offset, output_golden[3] + c_offset);
46  VERIFY_IS_APPROX(output(4) + c_offset, output_golden[4] + c_offset);
47  VERIFY_IS_APPROX(output(5) + c_offset, output_golden[5] + c_offset);
48  }
49  else {
50  VERIFY_IS_APPROX(output(0)+ c_offset, output_golden[0]+ c_offset);
51  VERIFY_IS_APPROX(output(1)+ c_offset, output_golden[2]+ c_offset);
52  VERIFY_IS_APPROX(output(2)+ c_offset, output_golden[4]+ c_offset);
53  VERIFY_IS_APPROX(output(3)+ c_offset, output_golden[1]+ c_offset);
54  VERIFY_IS_APPROX(output(4)+ c_offset, output_golden[3]+ c_offset);
55  VERIFY_IS_APPROX(output(5)+ c_offset, output_golden[5]+ c_offset);
56  }
57 }
58 
59 static void test_fft_complex_input_golden() {
60  Tensor<std::complex<float>, 1, ColMajor> input(5);
61  input(0) = std::complex<float>(1, 1);
62  input(1) = std::complex<float>(2, 2);
63  input(2) = std::complex<float>(3, 3);
64  input(3) = std::complex<float>(4, 4);
65  input(4) = std::complex<float>(5, 5);
66 
67  array<ptrdiff_t, 1> fft;
68  fft[0] = 0;
69 
70  Tensor<std::complex<float>, 1, ColMajor> forward_output_both_parts = input.fft<BothParts, FFT_FORWARD>(fft);
71  Tensor<std::complex<float>, 1, ColMajor> reverse_output_both_parts = input.fft<BothParts, FFT_REVERSE>(fft);
72 
73  Tensor<float, 1, ColMajor> forward_output_real_part = input.fft<RealPart, FFT_FORWARD>(fft);
74  Tensor<float, 1, ColMajor> reverse_output_real_part = input.fft<RealPart, FFT_REVERSE>(fft);
75 
76  Tensor<float, 1, ColMajor> forward_output_imag_part = input.fft<ImagPart, FFT_FORWARD>(fft);
77  Tensor<float, 1, ColMajor> reverse_output_imag_part = input.fft<ImagPart, FFT_REVERSE>(fft);
78 
79  VERIFY_IS_EQUAL(forward_output_both_parts.dimension(0), input.dimension(0));
80  VERIFY_IS_EQUAL(reverse_output_both_parts.dimension(0), input.dimension(0));
81 
82  VERIFY_IS_EQUAL(forward_output_real_part.dimension(0), input.dimension(0));
83  VERIFY_IS_EQUAL(reverse_output_real_part.dimension(0), input.dimension(0));
84 
85  VERIFY_IS_EQUAL(forward_output_imag_part.dimension(0), input.dimension(0));
86  VERIFY_IS_EQUAL(reverse_output_imag_part.dimension(0), input.dimension(0));
87 
88  std::complex<float> forward_golden_result[5];
89  std::complex<float> reverse_golden_result[5];
90 
91  forward_golden_result[0] = std::complex<float>(15.000000000000000,+15.000000000000000);
92  forward_golden_result[1] = std::complex<float>(-5.940954801177935, +0.940954801177934);
93  forward_golden_result[2] = std::complex<float>(-3.312299240582266, -1.687700759417735);
94  forward_golden_result[3] = std::complex<float>(-1.687700759417735, -3.312299240582266);
95  forward_golden_result[4] = std::complex<float>( 0.940954801177934, -5.940954801177935);
96 
97  reverse_golden_result[0] = std::complex<float>( 3.000000000000000, + 3.000000000000000);
98  reverse_golden_result[1] = std::complex<float>( 0.188190960235587, - 1.188190960235587);
99  reverse_golden_result[2] = std::complex<float>(-0.337540151883547, - 0.662459848116453);
100  reverse_golden_result[3] = std::complex<float>(-0.662459848116453, - 0.337540151883547);
101  reverse_golden_result[4] = std::complex<float>(-1.188190960235587, + 0.188190960235587);
102 
103  for(int i = 0; i < 5; ++i) {
104  VERIFY_IS_APPROX(forward_output_both_parts(i), forward_golden_result[i]);
105  VERIFY_IS_APPROX(forward_output_real_part(i), forward_golden_result[i].real());
106  VERIFY_IS_APPROX(forward_output_imag_part(i), forward_golden_result[i].imag());
107  }
108 
109  for(int i = 0; i < 5; ++i) {
110  VERIFY_IS_APPROX(reverse_output_both_parts(i), reverse_golden_result[i]);
111  VERIFY_IS_APPROX(reverse_output_real_part(i), reverse_golden_result[i].real());
112  VERIFY_IS_APPROX(reverse_output_imag_part(i), reverse_golden_result[i].imag());
113  }
114 }
115 
116 static void test_fft_real_input_golden() {
117  Tensor<float, 1, ColMajor> input(5);
118  input(0) = 1.0;
119  input(1) = 2.0;
120  input(2) = 3.0;
121  input(3) = 4.0;
122  input(4) = 5.0;
123 
124  array<ptrdiff_t, 1> fft;
125  fft[0] = 0;
126 
127  Tensor<std::complex<float>, 1, ColMajor> forward_output_both_parts = input.fft<BothParts, FFT_FORWARD>(fft);
128  Tensor<std::complex<float>, 1, ColMajor> reverse_output_both_parts = input.fft<BothParts, FFT_REVERSE>(fft);
129 
130  Tensor<float, 1, ColMajor> forward_output_real_part = input.fft<RealPart, FFT_FORWARD>(fft);
131  Tensor<float, 1, ColMajor> reverse_output_real_part = input.fft<RealPart, FFT_REVERSE>(fft);
132 
133  Tensor<float, 1, ColMajor> forward_output_imag_part = input.fft<ImagPart, FFT_FORWARD>(fft);
134  Tensor<float, 1, ColMajor> reverse_output_imag_part = input.fft<ImagPart, FFT_REVERSE>(fft);
135 
136  VERIFY_IS_EQUAL(forward_output_both_parts.dimension(0), input.dimension(0));
137  VERIFY_IS_EQUAL(reverse_output_both_parts.dimension(0), input.dimension(0));
138 
139  VERIFY_IS_EQUAL(forward_output_real_part.dimension(0), input.dimension(0));
140  VERIFY_IS_EQUAL(reverse_output_real_part.dimension(0), input.dimension(0));
141 
142  VERIFY_IS_EQUAL(forward_output_imag_part.dimension(0), input.dimension(0));
143  VERIFY_IS_EQUAL(reverse_output_imag_part.dimension(0), input.dimension(0));
144 
145  std::complex<float> forward_golden_result[5];
146  std::complex<float> reverse_golden_result[5];
147 
148 
149  forward_golden_result[0] = std::complex<float>( 15, 0);
150  forward_golden_result[1] = std::complex<float>(-2.5, +3.44095480117793);
151  forward_golden_result[2] = std::complex<float>(-2.5, +0.81229924058227);
152  forward_golden_result[3] = std::complex<float>(-2.5, -0.81229924058227);
153  forward_golden_result[4] = std::complex<float>(-2.5, -3.44095480117793);
154 
155  reverse_golden_result[0] = std::complex<float>( 3.0, 0);
156  reverse_golden_result[1] = std::complex<float>(-0.5, -0.688190960235587);
157  reverse_golden_result[2] = std::complex<float>(-0.5, -0.162459848116453);
158  reverse_golden_result[3] = std::complex<float>(-0.5, +0.162459848116453);
159  reverse_golden_result[4] = std::complex<float>(-0.5, +0.688190960235587);
160 
161  std::complex<float> c_offset(1.0, 1.0);
162  float r_offset = 1.0;
163 
164  for(int i = 0; i < 5; ++i) {
165  VERIFY_IS_APPROX(forward_output_both_parts(i) + c_offset, forward_golden_result[i] + c_offset);
166  VERIFY_IS_APPROX(forward_output_real_part(i) + r_offset, forward_golden_result[i].real() + r_offset);
167  VERIFY_IS_APPROX(forward_output_imag_part(i) + r_offset, forward_golden_result[i].imag() + r_offset);
168  }
169 
170  for(int i = 0; i < 5; ++i) {
171  VERIFY_IS_APPROX(reverse_output_both_parts(i) + c_offset, reverse_golden_result[i] + c_offset);
172  VERIFY_IS_APPROX(reverse_output_real_part(i) + r_offset, reverse_golden_result[i].real() + r_offset);
173  VERIFY_IS_APPROX(reverse_output_imag_part(i) + r_offset, reverse_golden_result[i].imag() + r_offset);
174  }
175 }
176 
177 
178 template <int DataLayout, typename RealScalar, bool isComplexInput, int FFTResultType, int FFTDirection, int TensorRank>
179 static void test_fft_real_input_energy() {
180 
181  Eigen::DSizes<ptrdiff_t, TensorRank> dimensions;
182  ptrdiff_t total_size = 1;
183  for (int i = 0; i < TensorRank; ++i) {
184  dimensions[i] = rand() % 20 + 1;
185  total_size *= dimensions[i];
186  }
187  const DSizes<ptrdiff_t, TensorRank> arr = dimensions;
188 
189  typedef typename internal::conditional<isComplexInput == true, std::complex<RealScalar>, RealScalar>::type InputScalar;
190 
191  Tensor<InputScalar, TensorRank, DataLayout> input;
192  input.resize(arr);
193  input.setRandom();
194 
195  array<ptrdiff_t, TensorRank> fft;
196  for (int i = 0; i < TensorRank; ++i) {
197  fft[i] = i;
198  }
199 
200  typedef typename internal::conditional<FFTResultType == Eigen::BothParts, std::complex<RealScalar>, RealScalar>::type OutputScalar;
201  Tensor<OutputScalar, TensorRank, DataLayout> output;
202  output = input.template fft<FFTResultType, FFTDirection>(fft);
203 
204  for (int i = 0; i < TensorRank; ++i) {
205  VERIFY_IS_EQUAL(output.dimension(i), input.dimension(i));
206  }
207 
208  RealScalar energy_original = 0.0;
209  RealScalar energy_after_fft = 0.0;
210 
211  for (int i = 0; i < total_size; ++i) {
212  energy_original += numext::abs2(input(i));
213  }
214 
215  for (int i = 0; i < total_size; ++i) {
216  energy_after_fft += numext::abs2(output(i));
217  }
218 
219  if(FFTDirection == FFT_FORWARD) {
220  VERIFY_IS_APPROX(energy_original, energy_after_fft / total_size);
221  }
222  else {
223  VERIFY_IS_APPROX(energy_original, energy_after_fft * total_size);
224  }
225 }
226 
227 template <typename RealScalar>
228 static void test_fft_non_power_of_2_round_trip(int exponent) {
229  int n = (1 << exponent) + 1;
230 
231  Eigen::DSizes<ptrdiff_t, 1> dimensions;
232  dimensions[0] = n;
233  const DSizes<ptrdiff_t, 1> arr = dimensions;
234  Tensor<RealScalar, 1, ColMajor, ptrdiff_t> input;
235 
236  input.resize(arr);
237  input.setRandom();
238 
239  array<int, 1> fft;
240  fft[0] = 0;
241 
242  Tensor<std::complex<RealScalar>, 1, ColMajor> forward =
243  input.template fft<BothParts, FFT_FORWARD>(fft);
244 
245  Tensor<RealScalar, 1, ColMajor, ptrdiff_t> output =
246  forward.template fft<RealPart, FFT_REVERSE>(fft);
247 
248  for (int i = 0; i < n; ++i) {
249  RealScalar tol = test_precision<RealScalar>() *
250  (std::abs(input[i]) + std::abs(output[i]) + 1);
251  VERIFY_IS_APPROX_OR_LESS_THAN(std::abs(input[i] - output[i]), tol);
252  }
253 }
254 
255 EIGEN_DECLARE_TEST(cxx11_tensor_fft) {
256  test_fft_complex_input_golden();
257  test_fft_real_input_golden();
258 
259  test_fft_2D_golden<ColMajor>();
260  test_fft_2D_golden<RowMajor>();
261 
262  test_fft_real_input_energy<ColMajor, float, true, Eigen::BothParts, FFT_FORWARD, 1>();
263  test_fft_real_input_energy<ColMajor, double, true, Eigen::BothParts, FFT_FORWARD, 1>();
264  test_fft_real_input_energy<ColMajor, float, false, Eigen::BothParts, FFT_FORWARD, 1>();
265  test_fft_real_input_energy<ColMajor, double, false, Eigen::BothParts, FFT_FORWARD, 1>();
266 
267  test_fft_real_input_energy<ColMajor, float, true, Eigen::BothParts, FFT_FORWARD, 2>();
268  test_fft_real_input_energy<ColMajor, double, true, Eigen::BothParts, FFT_FORWARD, 2>();
269  test_fft_real_input_energy<ColMajor, float, false, Eigen::BothParts, FFT_FORWARD, 2>();
270  test_fft_real_input_energy<ColMajor, double, false, Eigen::BothParts, FFT_FORWARD, 2>();
271 
272  test_fft_real_input_energy<ColMajor, float, true, Eigen::BothParts, FFT_FORWARD, 3>();
273  test_fft_real_input_energy<ColMajor, double, true, Eigen::BothParts, FFT_FORWARD, 3>();
274  test_fft_real_input_energy<ColMajor, float, false, Eigen::BothParts, FFT_FORWARD, 3>();
275  test_fft_real_input_energy<ColMajor, double, false, Eigen::BothParts, FFT_FORWARD, 3>();
276 
277  test_fft_real_input_energy<ColMajor, float, true, Eigen::BothParts, FFT_FORWARD, 4>();
278  test_fft_real_input_energy<ColMajor, double, true, Eigen::BothParts, FFT_FORWARD, 4>();
279  test_fft_real_input_energy<ColMajor, float, false, Eigen::BothParts, FFT_FORWARD, 4>();
280  test_fft_real_input_energy<ColMajor, double, false, Eigen::BothParts, FFT_FORWARD, 4>();
281 
282  test_fft_real_input_energy<RowMajor, float, true, Eigen::BothParts, FFT_FORWARD, 1>();
283  test_fft_real_input_energy<RowMajor, double, true, Eigen::BothParts, FFT_FORWARD, 1>();
284  test_fft_real_input_energy<RowMajor, float, false, Eigen::BothParts, FFT_FORWARD, 1>();
285  test_fft_real_input_energy<RowMajor, double, false, Eigen::BothParts, FFT_FORWARD, 1>();
286 
287  test_fft_real_input_energy<RowMajor, float, true, Eigen::BothParts, FFT_FORWARD, 2>();
288  test_fft_real_input_energy<RowMajor, double, true, Eigen::BothParts, FFT_FORWARD, 2>();
289  test_fft_real_input_energy<RowMajor, float, false, Eigen::BothParts, FFT_FORWARD, 2>();
290  test_fft_real_input_energy<RowMajor, double, false, Eigen::BothParts, FFT_FORWARD, 2>();
291 
292  test_fft_real_input_energy<RowMajor, float, true, Eigen::BothParts, FFT_FORWARD, 3>();
293  test_fft_real_input_energy<RowMajor, double, true, Eigen::BothParts, FFT_FORWARD, 3>();
294  test_fft_real_input_energy<RowMajor, float, false, Eigen::BothParts, FFT_FORWARD, 3>();
295  test_fft_real_input_energy<RowMajor, double, false, Eigen::BothParts, FFT_FORWARD, 3>();
296 
297  test_fft_real_input_energy<RowMajor, float, true, Eigen::BothParts, FFT_FORWARD, 4>();
298  test_fft_real_input_energy<RowMajor, double, true, Eigen::BothParts, FFT_FORWARD, 4>();
299  test_fft_real_input_energy<RowMajor, float, false, Eigen::BothParts, FFT_FORWARD, 4>();
300  test_fft_real_input_energy<RowMajor, double, false, Eigen::BothParts, FFT_FORWARD, 4>();
301 
302  test_fft_non_power_of_2_round_trip<float>(7);
303  test_fft_non_power_of_2_round_trip<double>(7);
304 }
real
float real
Definition: datatypes.h:10
Eigen::ColMajor
Definition: Constants.h:319
n
int n
Definition: BiCGSTAB_simple.cpp:1
test_fft_real_input_energy
static void test_fft_real_input_energy()
Definition: cxx11_tensor_fft.cpp:179
Eigen::TensorBase< Tensor< Scalar_, NumIndices_, Options_, IndexType_ > >::setRandom
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Tensor< Scalar_, NumIndices_, Options_, IndexType_ > & setRandom()
Definition: TensorBase.h:996
arr
py::array arr
Definition: test_numpy_array.cpp:77
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(cxx11_tensor_fft)
Definition: cxx11_tensor_fft.cpp:255
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
abs2
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
Definition: ArrayCwiseUnaryOps.h:80
test_fft_non_power_of_2_round_trip
static void test_fft_non_power_of_2_round_trip(int exponent)
Definition: cxx11_tensor_fft.cpp:228
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::Tensor::resize
EIGEN_DEVICE_FUNC void resize(const array< Index, NumIndices > &dimensions)
Definition: Tensor.h:447
VERIFY_IS_APPROX_OR_LESS_THAN
#define VERIFY_IS_APPROX_OR_LESS_THAN(a, b)
Definition: main.h:392
imag
EIGEN_DEVICE_FUNC const ImagReturnType imag() const
Definition: CommonCwiseUnaryOps.h:109
Eigen::Tensor::dimension
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index dimension(std::size_t n) const
Definition: Tensor.h:101
gtsam::tol
const G double tol
Definition: Group.h:86
test_fft_real_input_golden
static void test_fft_real_input_golden()
Definition: cxx11_tensor_fft.cpp:116
test_fft_complex_input_golden
static void test_fft_complex_input_golden()
Definition: cxx11_tensor_fft.cpp:59
test_fft_2D_golden
static void test_fft_2D_golden()
Definition: cxx11_tensor_fft.cpp:16
dimensions
const std::vector< size_t > dimensions
Definition: testVerticalBlockMatrix.cpp:27
abs
#define abs(x)
Definition: datatypes.h:17
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
DataLayout
static const int DataLayout
Definition: cxx11_tensor_image_patch_sycl.cpp:24
Eigen::Tensor
The tensor class.
Definition: Tensor.h:63

gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:07