forked from degauden/Elements
-
Notifications
You must be signed in to change notification settings - Fork 2
Expand file tree
/
Copy pathFftw.cpp
More file actions
99 lines (79 loc) · 3.19 KB
/
Copy pathFftw.cpp
File metadata and controls
99 lines (79 loc) · 3.19 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
/**
* @file Fftw.cpp
* @date January 6th, 2015
* @author Pierre Dubath
*
* @copyright 2012-2020 Euclid Science Ground Segment
*
* This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either version 3.0 of the License, or (at your option)
* any later version.
*
* This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to
* the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*
*/
#include <fftw3.h> // for fftw_destroy_plan, fftw_execute, fftw_plan_dft_1d, fftw_cleanup, fftw_complex,
// FFTW_ESTIMATE, fftw_plan, fftw_plan_s, FFTW_BACKWARD, FFTW_FORWARD
#include <cmath> // for cos
#include <cstdio> // for size_t, printf
#include <map> // for map
#include <string> // for allocator, string
#include <boost/format.hpp> // for basic_format, operator<<, format
#include "ElementsKernel/Main.h" // for MAIN_FOR
#include "ElementsKernel/MathConstants.h" // for pi
#include "ElementsKernel/Program.h" // for Program
#include "ElementsKernel/Unused.h" // for ELEMENTS_UNUSED
using std::map;
using std::string;
constexpr std::size_t N = 32;
namespace Elements::Examples {
class Fftw final : public Program {
public:
ExitCode mainMethod(ELEMENTS_UNUSED map<string, VariableValue>& args) override {
const auto log = Logging::getLogger("FftwExample");
fftw_complex in[N]; /* double [2] */
fftw_complex out[N];
fftw_complex in2[N];
using std::cos;
/* prepare a cosine wave */
for (size_t i = 0; i < N; i++) {
in[i][0] = cos(3.0 * 2.0 * Units::pi * static_cast<double>(i) / static_cast<double>(N));
in[i][1] = 0;
}
/* forward Fourier transform, save the result in 'out' */
fftw_plan_s* const p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
for (size_t i = 0; i < N; i++) {
log.info() << boost::format("freq: %3d %+9.5f %+9.5f I") % i % out[i][0] % out[i][1];
}
fftw_destroy_plan(p);
/* backward Fourier transform, save the result in 'in2' */
printf("\nInverse transform:\n");
fftw_plan_s* const q = fftw_plan_dft_1d(N, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(q);
/* normalize */
for (auto& i : in2) {
i[0] *= 1. / N;
i[1] *= 1. / N;
}
for (size_t i = 0; i < N; i++) {
log.info() << boost::format("recover: %3d %+9.5f %+9.5f I vs. %+9.5f %+9.5f I") % i % in[i][0] % in[i][1] %
in2[i][0] % in2[i][1];
}
fftw_destroy_plan(q);
fftw_cleanup();
log.info() << "This is the end of the test";
return ExitCode::OK;
}
};
} // namespace Elements::Examples
/**
* Implementation of a main using a base class macro
* This must be present in all Elements programs
*/
MAIN_FOR(Elements::Examples::Fftw)