SourceXtractorPlusPlus  0.16
Please provide a description of the project.
FFT.cpp
Go to the documentation of this file.
1 
17 /*
18  * @file FFT.cpp
19  * @date 11/09/18
20  * @author aalvarez
21  */
22 
23 #include "SEFramework/FFT/FFT.h"
24 #include <boost/thread/shared_mutex.hpp>
25 #include <cassert>
26 #include <fftw3.h>
27 #include <map>
28 
29 namespace SourceXtractor {
30 
35 boost::mutex fftw_global_plan_mutex{};
36 
42 
48 
49 int fftRoundDimension(int size) {
50  // Precomputed lookup table for the optimal dimension
51  static int optimal_dims[] = {
52  1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 18, 20, 20,
53  21, 22, 24, 24, 25, 26, 27, 28, 30, 30, 32, 32, 33, 35, 35, 36, 39, 39, 39, 40, 42,
54  42, 44, 44, 45, 48, 48, 48, 49, 50, 52, 52, 54, 54, 55, 56, 60, 60, 60, 60, 63, 63,
55  63, 64, 65, 66, 70, 70, 70, 70, 72, 72, 75, 75, 75, 77, 77, 78, 80, 80, 81, 84, 84,
56  84, 88, 88, 88, 88, 90, 90, 91, 96, 96, 96, 96, 96, 98, 98, 99, 100, 104, 104, 104, 104,
57  105, 108, 108, 108, 110, 110, 112, 112, 117, 117, 117, 117, 117, 120, 120, 120, 125, 125, 125, 125, 125,
58  126, 128, 128, 130, 130, 132, 132, 135, 135, 135, 140, 140, 140, 140, 140, 144, 144, 144, 144, 147, 147,
59  147, 150, 150, 150, 154, 154, 154, 154, 156, 156, 160, 160, 160, 160, 162, 162, 165, 165, 165, 168, 168,
60  168, 175, 175, 175, 175, 175, 175, 175, 176, 180, 180, 180, 180, 182, 182, 189, 189, 189, 189, 189, 189,
61  189, 192, 192, 192, 195, 195, 195, 196, 198, 198, 200, 200, 208, 208, 208, 208, 208, 208, 208, 208, 210,
62  210, 216, 216, 216, 216, 216, 216, 220, 220, 220, 220, 224, 224, 224, 224, 225, 231, 231, 231, 231, 231,
63  231, 234, 234, 234, 240, 240, 240, 240, 240, 240, 243, 243, 243, 245, 245, 250, 250, 250, 250, 250, 252,
64  252, 256, 256, 256, 256, 260, 260, 260, 260, 264, 264, 264, 264, 270, 270, 270, 270, 270, 270, 273, 273,
65  273, 275, 275, 280, 280, 280, 280, 280, 288, 288, 288, 288, 288, 288, 288, 288, 294, 294, 294, 294, 294,
66  294, 297, 297, 297, 300, 300, 300, 308, 308, 308, 308, 308, 308, 308, 308, 312, 312, 312, 312, 315, 315,
67  315, 320, 320, 320, 320, 320, 324, 324, 324, 324, 325, 330, 330, 330, 330, 330, 336, 336, 336, 336, 336,
68  336, 343, 343, 343, 343, 343, 343, 343, 350, 350, 350, 350, 350, 350, 350, 351, 352, 360, 360, 360, 360,
69  360, 360, 360, 360, 364, 364, 364, 364, 375, 375, 375, 375, 375, 375, 375, 375, 375, 375, 375, 378, 378,
70  378, 384, 384, 384, 384, 384, 384, 385, 390, 390, 390, 390, 390, 392, 392, 396, 396, 396, 396, 400, 400,
71  400, 400, 405, 405, 405, 405, 405, 416, 416, 416, 416, 416, 416, 416, 416, 416, 416, 416, 420, 420, 420,
72  420, 432, 432, 432, 432, 432, 432, 432, 432, 432, 432, 432, 432, 440, 440, 440, 440, 440, 440, 440, 440,
73  441, 448, 448, 448, 448, 448, 448, 448, 450, 450, 455, 455, 455, 455, 455, 462, 462, 462, 462, 462, 462,
74  462, 468, 468, 468, 468, 468, 468, 480, 480, 480, 480, 480, 480, 480, 480, 480, 480, 480, 480, 486, 486,
75  486, 486, 486, 486, 490, 490, 490, 490, 495, 495, 495, 495, 495, 500, 500, 500, 500, 500, 504, 504, 504,
76  504, 512, 512, 512, 512, 512, 512, 512, 512, 520, 520, 520, 520, 520, 520, 520, 520, 525, 525, 525, 525,
77  525, 528, 528, 528, 539, 539, 539, 539, 539, 539, 539, 539, 539, 539, 539, 540, 546, 546, 546, 546, 546,
78  546, 550, 550, 550, 550, 560, 560, 560, 560, 560, 560, 560, 560, 560, 560, 567, 567, 567, 567, 567, 567,
79  567, 576, 576, 576, 576, 576, 576, 576, 576, 576, 585, 585, 585, 585, 585, 585, 585, 585, 585, 588, 588,
80  588, 594, 594, 594, 594, 594, 594, 600, 600, 600, 600, 600, 600, 616, 616, 616, 616, 616, 616, 616, 616,
81  616, 616, 616, 616, 616, 616, 616, 616, 624, 624, 624, 624, 624, 624, 624, 624, 625, 630, 630, 630, 630,
82  630, 637, 637, 637, 637, 637, 637, 637, 640, 640, 640, 648, 648, 648, 648, 648, 648, 648, 648, 650, 650,
83  660, 660, 660, 660, 660, 660, 660, 660, 660, 660, 672, 672, 672, 672, 672, 672, 672, 672, 672, 672, 672,
84  672, 675, 675, 675, 686, 686, 686, 686, 686, 686, 686, 686, 686, 686, 686, 693, 693, 693, 693, 693, 693,
85  693, 700, 700, 700, 700, 700, 700, 700, 702, 702, 704, 704, 720, 720, 720, 720, 720, 720, 720, 720, 720,
86  720, 720, 720, 720, 720, 720, 720, 728, 728, 728, 728, 728, 728, 728, 728, 729, 735, 735, 735, 735, 735,
87  735, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 750, 756, 756, 756, 756, 756,
88  756, 768, 768, 768, 768, 768, 768, 768, 768, 768, 768, 768, 768, 770, 770, 780, 780, 780, 780, 780, 780,
89  780, 780, 780, 780, 784, 784, 784, 784, 792, 792, 792, 792, 792, 792, 792, 792, 800, 800, 800, 800, 800,
90  800, 800, 800, 810, 810, 810, 810, 810, 810, 810, 810, 810, 810, 819, 819, 819, 819, 819, 819, 819, 819,
91  819, 825, 825, 825, 825, 825, 825, 832, 832, 832, 832, 832, 832, 832, 840, 840, 840, 840, 840, 840, 840,
92  840, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864,
93  864, 864, 864, 864, 875, 875, 875, 875, 875, 875, 875, 875, 875, 875, 875, 880, 880, 880, 880, 880, 882,
94  882, 891, 891, 891, 891, 891, 891, 891, 891, 891, 896, 896, 896, 896, 896, 900, 900, 900, 900, 910, 910,
95  910, 910, 910, 910, 910, 910, 910, 910, 924, 924, 924, 924, 924, 924, 924, 924, 924, 924, 924, 924, 924,
96  924, 936, 936, 936, 936, 936, 936, 936, 936, 936, 936, 936, 936, 945, 945, 945, 945, 945, 945, 945, 945,
97  945, 960, 960, 960, 960, 960, 960, 960, 960, 960, 960, 960, 960, 960, 960, 960, 972, 972, 972, 972, 972,
98  972, 972, 972, 972, 972, 972, 972, 975, 975, 975, 980, 980, 980, 980, 980, 990, 990, 990, 990, 990, 990,
99  990, 990, 990, 990, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1008, 1008, 1008, 1008, 1008, 1008, 1008,
100  1008, 1024, 1024, 1024,
101  };
102  static int n_optimal_dims = sizeof(optimal_dims) / sizeof(int);
103 
104  if (size < n_optimal_dims) {
105  return optimal_dims[size];
106  }
107  return (size / 512 + (size % 512 != 0)) * 512;
108 }
109 
110 template <typename T>
111 auto FFT<T>::createForwardPlan(int width, int height, std::vector<T>& inout) -> plan_ptr_t {
112  size_t phy_height = height;
113  size_t phy_width = 2 * (width / 2 + 1);
114  size_t mem_size = phy_height * phy_width;
115 
116  // Make sure the buffers are big enough
117  if (inout.size() < mem_size) {
118  inout.resize(mem_size);
119  }
120 
121  // Cache plan, as they can be reused
122  static boost::shared_mutex mutex;
123  static std::map<std::tuple<int, int>, plan_ptr_t> plan_cache;
124 
125  boost::upgrade_lock<boost::shared_mutex> read_lock{mutex};
126 
127  auto pi = plan_cache.find(std::make_tuple(width, height));
128  if (pi != plan_cache.end()) {
129  return pi->second;
130  }
131 
132  // No available plan yet, so get one from FFTW
133  boost::upgrade_to_unique_lock<boost::shared_mutex> write_lock{read_lock};
134  boost::lock_guard<boost::mutex> lock_planner{fftw_global_plan_mutex};
135 
136  pi = plan_cache.emplace(
137  std::make_tuple(width, height),
138  plan_ptr_t{
139  fftw_traits::func_plan_fwd(
140  height, width, // n0, n1
141  inout.data(), reinterpret_cast<complex_t*>(inout.data()), // in, out
142  FFTW_ESTIMATE | FFTW_DESTROY_INPUT // flags
143  ),
144  fftw_traits::func_destroy_plan}
145  ).first;
146 
147  return pi->second;
148 }
149 
150 template <typename T>
151 auto FFT<T>::createInversePlan(int width, int height, std::vector<T>& inout) -> plan_ptr_t {
152  size_t phy_height = height;
153  size_t phy_width = 2 * (width / 2 + 1);
154  size_t mem_size = phy_height * phy_width;
155 
156  // Make sure the buffers are big enough
157  if (inout.size() < mem_size) {
158  inout.resize(mem_size);
159  }
160 
161  // Cache plan, as they can be reused
162  static boost::shared_mutex mutex;
163  static std::map<std::tuple<int, int>, plan_ptr_t> plan_cache;
164 
165  boost::upgrade_lock<boost::shared_mutex> read_lock{mutex};
166 
167  auto pi = plan_cache.find(std::make_tuple(width, height));
168  if (pi != plan_cache.end()) {
169  return pi->second;
170  }
171 
172  // No available plan yet, so get one from FFTW
173  boost::upgrade_to_unique_lock<boost::shared_mutex> write_lock{read_lock};
174  boost::lock_guard<boost::mutex> lock_planner{fftw_global_plan_mutex};
175 
176  pi = plan_cache.emplace(
177  std::make_tuple(width, height),
178  plan_ptr_t{
179  fftw_traits::func_plan_inv(
180  height, width, // n0, n1
181  reinterpret_cast<complex_t*>(inout.data()), inout.data(), // in, out
182  FFTW_ESTIMATE | FFTW_DESTROY_INPUT // flags
183  ),
184  fftw_traits::func_destroy_plan}
185  ).first;
186 
187  return pi->second;
188 }
189 
190 template <typename T>
192  fftw_traits::func_execute_fwd(plan.get(), inout.data(), reinterpret_cast<complex_t*>(inout.data()));
193 }
194 
195 template <typename T>
197  fftw_traits::func_execute_inv(plan.get(), reinterpret_cast<complex_t*>(inout.data()), inout.data());
198 }
199 
200 template struct FFT<float>;
201 template struct FFT<double>;
202 
203 } // namespace SourceXtractor
T data(T... args)
T emplace(T... args)
T end(T... args)
T find(T... args)
T get(T... args)
T make_tuple(T... args)
constexpr double pi
boost::mutex fftw_global_plan_mutex
Definition: FFT.cpp:35
int fftRoundDimension(int size)
Definition: FFT.cpp:49
decltype(fftw_plan_dft_c2r_2d) typedef func_plan_inv_t
Definition: FFT.h:67
static func_plan_inv_t * func_plan_inv
Definition: FFT.h:73
static func_destroy_plan_t * func_destroy_plan
Definition: FFT.h:74
decltype(fftw_destroy_plan) typedef func_destroy_plan_t
Definition: FFT.h:68
static func_execute_inv_t * func_execute_inv
Definition: FFT.h:76
static func_plan_fwd_t * func_plan_fwd
Definition: FFT.h:72
decltype(fftw_execute_dft_c2r) typedef func_execute_inv_t
Definition: FFT.h:70
decltype(fftw_plan_dft_r2c_2d) typedef func_plan_fwd_t
Definition: FFT.h:66
static func_execute_fwd_t * func_execute_fwd
Definition: FFT.h:75
decltype(fftw_execute_dft_r2c) typedef func_execute_fwd_t
Definition: FFT.h:69
static func_execute_fwd_t * func_execute_fwd
Definition: FFT.h:55
static func_plan_fwd_t * func_plan_fwd
Definition: FFT.h:52
decltype(fftwf_execute_dft_r2c) typedef func_execute_fwd_t
Definition: FFT.h:49
static func_destroy_plan_t * func_destroy_plan
Definition: FFT.h:54
static func_execute_inv_t * func_execute_inv
Definition: FFT.h:56
static func_plan_inv_t * func_plan_inv
Definition: FFT.h:53
decltype(fftwf_plan_dft_r2c_2d) typedef func_plan_fwd_t
Definition: FFT.h:46
decltype(fftwf_plan_dft_c2r_2d) typedef func_plan_inv_t
Definition: FFT.h:47
decltype(fftwf_execute_dft_c2r) typedef func_execute_inv_t
Definition: FFT.h:50
decltype(fftwf_destroy_plan) typedef func_destroy_plan_t
Definition: FFT.h:48
Wraps the FFTW API with a more C++ like one.
Definition: FFT.h:86
static void executeInverse(plan_ptr_t &plan, std::vector< T > &inout)
Definition: FFT.cpp:196
static plan_ptr_t createForwardPlan(int width, int height, std::vector< T > &inout)
Definition: FFT.cpp:111
fftw_traits::complex_t complex_t
Definition: FFT.h:92
static plan_ptr_t createInversePlan(int width, int height, std::vector< T > &inout)
Definition: FFT.cpp:151
static void executeForward(plan_ptr_t &plan, std::vector< T > &inout)
Definition: FFT.cpp:191