1 // Licensed GNU LGPL v3 or later: http://www.gnu.org/licenses/lgpl.html
2
3 #include "smfft.hh"
4 #include "smutils.hh"
5 #include <algorithm>
6 #include <map>
7 #include <mutex>
8 #include "config.h"
9
10 #include <glib.h>
11 #include <glib/gstdio.h>
12 #include <string.h>
13 #include <unistd.h>
14
15 using namespace SpectMorph;
16 using std::map;
17 using std::string;
18
19 static bool enable_gsl_fft = false;
20 static bool randomize_new_fft_arrays = false;
21
22 void
use_gsl_fft(bool new_enable_gsl_fft)23 FFT::use_gsl_fft (bool new_enable_gsl_fft)
24 {
25 enable_gsl_fft = new_enable_gsl_fft;
26 }
27
28 void
debug_randomize_new_arrays(bool b)29 FFT::debug_randomize_new_arrays (bool b)
30 {
31 randomize_new_fft_arrays = b;
32 }
33
34 #if SPECTMORPH_HAVE_FFTW
35
36 #include <fftw3.h>
37
38 static void save_wisdom();
39
40 /*
41 * we force that only one thread at a time will be in planning mode
42 * - fftw planner is not threadsafe
43 * - neither is our code that generates plans
44 */
45 static std::mutex fftw_plan_mutex;
46 static std::mutex plan_map_mutex;
47
48 static fftwf_plan&
read_plan_map_threadsafe(std::map<int,fftwf_plan> & plan_map,size_t N)49 read_plan_map_threadsafe (std::map<int, fftwf_plan>& plan_map, size_t N)
50 {
51 /* std::map access is not threadsafe */
52 std::lock_guard<std::mutex> lg (plan_map_mutex);
53 return plan_map[N];
54 }
55
56 float *
new_array_float(size_t N)57 FFT::new_array_float (size_t N)
58 {
59 const size_t N_2 = N + 2; /* extra space for r2c extra complex output */
60
61 float *result = (float *) fftwf_malloc (sizeof (float) * N_2);
62
63 if (randomize_new_fft_arrays)
64 {
65 for (size_t i = 0; i < N_2; i++)
66 result[i] = g_random_double_range (-1, 1);
67 }
68 return result;
69 }
70
71 void
free_array_float(float * f)72 FFT::free_array_float (float *f)
73 {
74 fftwf_free (f);
75 }
76
77 static map<int, fftwf_plan> fftar_float_plan;
78
79 static int
plan_flags(FFT::PlanMode plan_mode)80 plan_flags (FFT::PlanMode plan_mode)
81 {
82 switch (plan_mode)
83 {
84 case FFT::PLAN_PATIENT: return (FFTW_PATIENT | FFTW_PRESERVE_INPUT | FFTW_WISDOM_ONLY);
85 case FFT::PLAN_ESTIMATE: return (FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
86 default: g_assert_not_reached();
87 }
88 }
89
90 void
fftar_float(size_t N,float * in,float * out,PlanMode plan_mode)91 FFT::fftar_float (size_t N, float *in, float *out, PlanMode plan_mode)
92 {
93 fftwf_plan& plan = read_plan_map_threadsafe (fftar_float_plan, N);
94
95 if (!plan)
96 {
97 std::lock_guard<std::mutex> lg (fftw_plan_mutex);
98 float *plan_in = new_array_float (N);
99 float *plan_out = new_array_float (N);
100 plan = fftwf_plan_dft_r2c_1d (N, plan_in, (fftwf_complex *) plan_out, plan_flags (plan_mode));
101 if (!plan) /* missing from wisdom -> create plan and save it */
102 {
103 plan = fftwf_plan_dft_r2c_1d (N, plan_in, (fftwf_complex *) plan_out, plan_flags (plan_mode) & ~FFTW_WISDOM_ONLY);
104 save_wisdom();
105 }
106 free_array_float (plan_out);
107 free_array_float (plan_in);
108 }
109 fftwf_execute_dft_r2c (plan, in, (fftwf_complex *) out);
110
111 out[1] = out[N];
112 }
113
114 static map<int, fftwf_plan> fftsr_float_plan;
115
116 void
fftsr_float(size_t N,float * in,float * out,PlanMode plan_mode)117 FFT::fftsr_float (size_t N, float *in, float *out, PlanMode plan_mode)
118 {
119 fftwf_plan& plan = read_plan_map_threadsafe (fftsr_float_plan, N);
120
121 if (!plan)
122 {
123 std::lock_guard<std::mutex> lg (fftw_plan_mutex);
124 float *plan_in = new_array_float (N);
125 float *plan_out = new_array_float (N);
126 plan = fftwf_plan_dft_c2r_1d (N, (fftwf_complex *) plan_in, plan_out, plan_flags (plan_mode));
127 if (!plan) /* missing from wisdom -> create plan and save it */
128 {
129 plan = fftwf_plan_dft_c2r_1d (N, (fftwf_complex *) plan_in, plan_out, plan_flags (plan_mode) & ~FFTW_WISDOM_ONLY);
130 save_wisdom();
131 }
132 free_array_float (plan_out);
133 free_array_float (plan_in);
134 }
135 in[N] = in[1];
136 in[N+1] = 0;
137 in[1] = 0;
138
139 fftwf_execute_dft_c2r (plan, (fftwf_complex *)in, out);
140
141 in[1] = in[N]; // we need to preserve the input array
142 }
143
144 static map<int, fftwf_plan> fftsr_destructive_float_plan;
145
146 void
fftsr_destructive_float(size_t N,float * in,float * out,PlanMode plan_mode)147 FFT::fftsr_destructive_float (size_t N, float *in, float *out, PlanMode plan_mode)
148 {
149 fftwf_plan& plan = read_plan_map_threadsafe (fftsr_destructive_float_plan, N);
150
151 if (!plan)
152 {
153 std::lock_guard<std::mutex> lg (fftw_plan_mutex);
154 int xplan_flags = plan_flags (plan_mode) & ~FFTW_PRESERVE_INPUT;
155 float *plan_in = new_array_float (N);
156 float *plan_out = new_array_float (N);
157 plan = fftwf_plan_dft_c2r_1d (N, (fftwf_complex *) plan_in, plan_out, xplan_flags);
158 if (!plan) /* missing from wisdom -> create plan and save it */
159 {
160 plan = fftwf_plan_dft_c2r_1d (N, (fftwf_complex *) plan_in, plan_out,
161 xplan_flags & ~FFTW_WISDOM_ONLY);
162 save_wisdom();
163 }
164 free_array_float (plan_out);
165 free_array_float (plan_in);
166 }
167 in[N] = in[1];
168 in[N+1] = 0;
169 in[1] = 0;
170
171 fftwf_execute_dft_c2r (plan, (fftwf_complex *)in, out);
172 }
173
174 static map<int, fftwf_plan> fftac_float_plan;
175
176 void
fftac_float(size_t N,float * in,float * out,PlanMode plan_mode)177 FFT::fftac_float (size_t N, float *in, float *out, PlanMode plan_mode)
178 {
179 fftwf_plan& plan = read_plan_map_threadsafe (fftac_float_plan, N);
180 if (!plan)
181 {
182 std::lock_guard<std::mutex> lg (fftw_plan_mutex);
183 float *plan_in = new_array_float (N * 2);
184 float *plan_out = new_array_float (N * 2);
185
186 plan = fftwf_plan_dft_1d (N, (fftwf_complex *) plan_in, (fftwf_complex *) plan_out,
187 FFTW_FORWARD, plan_flags (plan_mode));
188 if (!plan) /* missing from wisdom -> create plan and save it */
189 {
190 plan = fftwf_plan_dft_1d (N, (fftwf_complex *) plan_in, (fftwf_complex *) plan_out,
191 FFTW_FORWARD, plan_flags (plan_mode) & ~FFTW_WISDOM_ONLY);
192 save_wisdom();
193 }
194 free_array_float (plan_out);
195 free_array_float (plan_in);
196 }
197
198 fftwf_execute_dft (plan, (fftwf_complex *)in, (fftwf_complex *)out);
199 }
200
201 static map<int, fftwf_plan> fftsc_float_plan;
202
203 void
fftsc_float(size_t N,float * in,float * out,PlanMode plan_mode)204 FFT::fftsc_float (size_t N, float *in, float *out, PlanMode plan_mode)
205 {
206 fftwf_plan& plan = read_plan_map_threadsafe (fftsc_float_plan, N);
207 if (!plan)
208 {
209 std::lock_guard<std::mutex> lg (fftw_plan_mutex);
210 float *plan_in = new_array_float (N * 2);
211 float *plan_out = new_array_float (N * 2);
212
213 plan = fftwf_plan_dft_1d (N, (fftwf_complex *) plan_in, (fftwf_complex *) plan_out,
214 FFTW_BACKWARD, plan_flags (plan_mode));
215 if (!plan) /* missing from wisdom -> create plan and save it */
216 {
217 plan = fftwf_plan_dft_1d (N, (fftwf_complex *) plan_in, (fftwf_complex *) plan_out,
218 FFTW_BACKWARD, plan_flags (plan_mode) & ~FFTW_WISDOM_ONLY);
219 save_wisdom();
220 }
221 free_array_float (plan_out);
222 free_array_float (plan_in);
223 }
224 fftwf_execute_dft (plan, (fftwf_complex *)in, (fftwf_complex *)out);
225 }
226
227 static string
wisdom_filename()228 wisdom_filename()
229 {
230 const char *hostname = g_get_host_name();
231 return sm_get_user_dir (USER_DIR_DATA) + string ("/.fftw_wisdom_") + hostname;
232 }
233
234 static void
save_wisdom()235 save_wisdom()
236 {
237 /* detect if we're running in valgrind - in this case newly accumulated wisdom is probably flawed */
238 bool valgrind = false;
239
240 FILE *maps = fopen (string_printf ("/proc/%d/maps", getpid()).c_str(), "r");
241 if (maps)
242 {
243 char buffer[1024];
244 while (fgets (buffer, 1024, maps))
245 {
246 if (strstr (buffer, "vgpreload"))
247 valgrind = true;
248 }
249 fclose (maps);
250 }
251 if (valgrind)
252 {
253 printf ("FFT::save_wisdom(): not saving fft wisdom (running under valgrind)\n");
254 return;
255 }
256 /* atomically replace old wisdom file with new wisdom file
257 *
258 * its theoretically possible (but highly unlikely) that we leak a *wisdom*.new.12345 file
259 */
260 string new_wisdom_filename = string_printf ("%s.new.%d", wisdom_filename().c_str(), getpid());
261 FILE *outfile = fopen (new_wisdom_filename.c_str(), "w");
262 if (outfile)
263 {
264 fftwf_export_wisdom_to_file (outfile);
265 fclose (outfile);
266 g_rename (new_wisdom_filename.c_str(), wisdom_filename().c_str());
267 }
268 }
269
270 static void
load_wisdom()271 load_wisdom()
272 {
273 FILE *infile = fopen (wisdom_filename().c_str(), "r");
274 if (infile)
275 {
276 fftwf_import_wisdom_from_file (infile);
277 fclose (infile);
278 }
279 }
280
281 void
init()282 FFT::init()
283 {
284 /* If an "application" is structured as a set of "plugins" which are unaware
285 * of each other, we need to force a thread-safe planner to avoid crashes.
286 *
287 * As SpectMorph is to be run as plugin in many situations, we need this:
288 */
289 #if SPECTMORPH_HAVE_FFTW_THREADSAFE
290 fftwf_make_planner_thread_safe();
291 #endif
292 load_wisdom();
293 }
294
295 void
cleanup()296 FFT::cleanup()
297 {
298 typedef map<int, fftwf_plan> PlanMap;
299
300 auto cleanup_plans = [](PlanMap& plan_map) {
301 for (auto& plan_entry : plan_map)
302 fftwf_destroy_plan (plan_entry.second);
303
304 plan_map.clear();
305 };
306 cleanup_plans (fftar_float_plan);
307 cleanup_plans (fftsr_float_plan);
308 cleanup_plans (fftsr_destructive_float_plan);
309 cleanup_plans (fftac_float_plan);
310 cleanup_plans (fftsc_float_plan);
311 }
312
313 #else
314
315 #error "building without FFTW is not supported currently"
316
317 #endif
318