SigPack - the C++ signal processing library
fftw.h
Go to the documentation of this file.
1 // This Source Code Form is subject to the terms of the Mozilla Public
2 // License, v. 2.0. If a copy of the MPL was not distributed with this
3 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
4 #ifndef SP_FFTW_H
5 #define SP_FFTW_H
6 #include <fftw3.h>
7 
8 namespace sp
9 {
18 
25 class FFTW
26 {
27 private:
28  fftw_plan pl_fft;
29  fftw_plan pl_ifft;
30  fftw_plan pl_fft_cx;
31  fftw_plan pl_ifft_cx;
32  unsigned int N;
33  unsigned int R, C;
34  int alg;
35  int export_alg;
38 public:
44  FFTW(unsigned int _N, int _alg = FFTW_ESTIMATE)
45  {
46  N = _N;
47  R = 0;
48  C = 0;
49  alg = _alg;
50  export_alg = FFTW_PATIENT;
51  pl_fft = NULL;
52  pl_ifft = NULL;
53  pl_fft_cx = NULL;
54  pl_ifft_cx = NULL;
55  }
56 
63  FFTW(unsigned int _R, unsigned int _C, int _alg)
64  {
65  R = _R;
66  C = _C;
67  N = 0;
68  alg = _alg;
69  export_alg = FFTW_PATIENT;
70  pl_fft = NULL;
71  pl_ifft = NULL;
72  pl_fft_cx = NULL;
73  pl_ifft_cx = NULL;
74  }
75 
80  {
81  if (pl_fft != NULL)
82  fftw_destroy_plan(pl_fft);
83  if (pl_ifft != NULL)
84  fftw_destroy_plan(pl_ifft);
85  if (pl_fft_cx != NULL)
86  fftw_destroy_plan(pl_fft_cx);
87  if (pl_ifft_cx != NULL)
88  fftw_destroy_plan(pl_ifft_cx);
89  fftw_cleanup();
90  }
91 
97  void fft_cx(arma::cx_vec &x, arma::cx_vec &Pxx)
98  {
99  fftw_complex *in = reinterpret_cast<fftw_complex *>(x.memptr());
100  fftw_complex *out = reinterpret_cast<fftw_complex *>(Pxx.memptr());
101  if (pl_fft_cx == NULL)
102  {
103  pl_fft_cx = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, alg);
104  if (pl_fft_cx == NULL)
105  {
106  err_handler("Unable to create complex data FFTW plan");
107  }
108  }
109  fftw_execute_dft(pl_fft_cx, in, out);
110  }
111 
117  arma::cx_vec fft_cx(arma::cx_vec &x)
118  {
119  arma::cx_vec Pxx(N);
120  fft_cx(x, Pxx);
121  return Pxx;
122  }
123 
129  void ifft_cx(arma::cx_vec &Pxx, arma::cx_vec &x)
130  {
131  fftw_complex *in = reinterpret_cast<fftw_complex *>(Pxx.memptr());
132  fftw_complex *out = reinterpret_cast<fftw_complex *>(x.memptr());
133  if (pl_ifft_cx == NULL)
134  {
135  pl_ifft_cx = fftw_plan_dft_1d(N, in, out, FFTW_BACKWARD, alg);
136  if (pl_ifft_cx == NULL)
137  {
138  err_handler("Unable to create complex data IFFTW plan");
139  }
140  }
141  fftw_execute_dft(pl_ifft_cx, in, out);
142  x /= N;
143  }
144 
150  arma::cx_vec ifft_cx(arma::cx_vec &Pxx)
151  {
152  arma::cx_vec x(N);
153  ifft_cx(Pxx, x);
154  return x;
155  }
156 
162  void fft(arma::vec &x, arma::cx_vec &Pxx)
163  {
164  double *in = x.memptr();
165  fftw_complex *out = reinterpret_cast<fftw_complex *>(Pxx.memptr());
166  if (pl_fft == NULL)
167  {
168  pl_fft = fftw_plan_dft_r2c_1d(N, in, out, alg);
169  if (pl_fft == NULL)
170  {
171  err_handler("Unable to create real data FFTW plan");
172  }
173  }
174 
175  fftw_execute_dft_r2c(pl_fft, in, out);
176  int offset = static_cast<int>(ceil(N / 2.0));
177  int n_elem = N - offset;
178  for (int i = 0; i < n_elem; ++i)
179  {
180  Pxx(offset + i) = std::conj(Pxx(n_elem - i));
181  }
182  }
183 
188  arma::cx_vec fft(arma::vec &x)
189  {
190  arma::cx_vec Pxx(N);
191  fft(x, Pxx);
192  return Pxx;
193  }
194 
200  void ifft(arma::cx_vec &Pxx, arma::vec &x)
201  {
202  fftw_complex *in = reinterpret_cast<fftw_complex *>(Pxx.memptr());
203  double *out = x.memptr();
204  if (pl_ifft == NULL)
205  {
206  pl_ifft = fftw_plan_dft_c2r_1d(N, in, out, alg);
207  if (pl_ifft == NULL)
208  {
209  err_handler("Unable to create real data IFFTW plan");
210  }
211  }
212  fftw_execute_dft_c2r(pl_ifft, in, out);
213  x /= N;
214  }
215 
221  arma::vec ifft(arma::cx_vec &Pxx)
222  {
223  arma::vec x(N);
224  ifft(Pxx, x);
225  return x;
226  }
227 
233  void fft2(arma::mat &x, arma::cx_mat &Pxx)
234  {
235  arma::cx_mat Ptmp(R / 2 + 1, C, arma::fill::ones);
236  double *in = x.memptr();
237  fftw_complex *out = reinterpret_cast<fftw_complex *>(Ptmp.memptr());
238  if (pl_fft == NULL)
239  {
240  pl_fft = fftw_plan_dft_r2c_2d(
241  C, R, in, out,
242  alg); // Column to row-major order trick: switch C and R
243  if (pl_fft == NULL)
244  {
245  err_handler("Unable to create real data FFTW plan");
246  }
247  }
248 
249  fftw_execute_dft_r2c(pl_fft, in, out);
250 
251  // Reshape Pxx - upper half
252  std::complex<double> *ptr = reinterpret_cast<std::complex<double> *>(out);
253  const unsigned int Roff = R / 2 + 1;
254  for (unsigned int r = 0; r < Roff; r++)
255  {
256  for (unsigned int c = 0; c < C; c++)
257  {
258  Pxx(r, c) = ptr[r + c * Roff];
259  }
260  }
261  // Reshape Pxx - conj symmetry
262  for (unsigned int r = Roff; r < R; r++)
263  {
264  Pxx(r, 0) = conj(ptr[R - r]);
265  for (unsigned int c = 1; c < C; c++)
266  {
267  Pxx(r, c) = conj(ptr[R - r + (C - c) * Roff]);
268  }
269  }
270  }
271 
276  arma::cx_mat fft2(arma::mat &x)
277  {
278  arma::cx_mat Pxx(R, C, arma::fill::ones);
279  fft2(x, Pxx);
280  return Pxx;
281  }
282 
288  void ifft2(arma::cx_mat &Pxx, arma::mat &x)
289  {
290  // Reshape to row-major format
291  unsigned int Roff = R / 2 + 1;
292  arma::cx_mat Ptmp(Roff, C);
293  std::complex<double> *ptr =
294  reinterpret_cast<std::complex<double> *>(Ptmp.memptr());
295  for (unsigned int r = 0; r < Roff; r++)
296  {
297  for (unsigned int c = 0; c < C; c++)
298  {
299  ptr[r + c * Roff] = Pxx(r, c);
300  }
301  }
302 
303  fftw_complex *in = reinterpret_cast<fftw_complex *>(Ptmp.memptr());
304  double *out = x.memptr();
305  if (pl_ifft == NULL)
306  {
307  pl_ifft = fftw_plan_dft_c2r_2d(C, R, in, out, alg);
308  if (pl_ifft == NULL)
309  {
310  err_handler("Unable to create real data IFFTW plan");
311  }
312  }
313  fftw_execute_dft_c2r(pl_ifft, in, out);
314  x /= (R * C);
315  }
316 
322  arma::mat ifft2(arma::cx_mat &Pxx)
323  {
324  arma::mat x(R, C);
325  ifft2(Pxx, x);
326  return x;
327  }
328 
333  void import_wisdom_string(const std::string wisd)
334  {
335  int res = fftw_import_wisdom_from_string(wisd.c_str());
336  if (res == 0)
337  err_handler("Unable to import wisdom from string!");
338  }
339 
344  void import_wisdom_file(const std::string fname)
345  {
346  int res = fftw_import_wisdom_from_filename(fname.c_str());
347  if (res == 0)
348  err_handler("Unable to import wisdom from file!");
349  }
350 
355  void export_wisdom_fft(const std::string fname)
356  {
357  fftw_plan pl_w = NULL;
358  double *x_r;
359  fftw_complex *x_cx1;
360 
361  if (R == 0 || C == 0) // 1D
362  {
363  x_r = fftw_alloc_real(N);
364  x_cx1 = fftw_alloc_complex(N);
365 
366  // Replan using wisdom
367  pl_w = fftw_plan_dft_r2c_1d(N, x_r, x_cx1, export_alg);
368  }
369  else // 2D
370  {
371  x_r = fftw_alloc_real(R * C);
372  x_cx1 = fftw_alloc_complex(R * C);
373 
374  // Replan using wisdom
375  pl_w = fftw_plan_dft_r2c_2d(C, R, x_r, x_cx1, export_alg);
376  }
377  if (pl_w == NULL)
378  {
379  err_handler("Unable to create real data FFTW plan");
380  }
381 
382  // Export
383  if (fftw_export_wisdom_to_filename(fname.c_str()) == 0)
384  {
385  err_handler("Could not export wisdom to file!");
386  }
387 
388  fftw_destroy_plan(pl_w);
389  fftw_free(x_r);
390  fftw_free(x_cx1);
391  }
392 
397  void export_wisdom_ifft(const std::string fname)
398  {
399  fftw_plan pl_w = NULL;
400  double *x_r;
401  fftw_complex *x_cx1;
402 
403  if (R == 0 || C == 0) // 1D
404  {
405  x_r = fftw_alloc_real(N);
406  x_cx1 = fftw_alloc_complex(N);
407 
408  // Replan using wisdom
409  pl_w = fftw_plan_dft_c2r_1d(N, x_cx1, x_r, export_alg);
410  }
411  else // 2D
412  {
413  x_r = fftw_alloc_real(R * C);
414  x_cx1 = fftw_alloc_complex(R * C);
415 
416  // Replan using wisdom
417  pl_w = fftw_plan_dft_c2r_2d(C, R, x_cx1, x_r, export_alg);
418  }
419 
420  if (pl_w == NULL)
421  {
422  err_handler("Unable to create real data FFTW plan");
423  }
424 
425  // Export
426  if (fftw_export_wisdom_to_filename(fname.c_str()) == 0)
427  {
428  err_handler("Could not export wisdom to file!");
429  }
430 
431  fftw_destroy_plan(pl_w);
432  fftw_free(x_r);
433  fftw_free(x_cx1);
434  }
435 
440  void export_wisdom_fft_cx(const std::string fname)
441  {
442  fftw_plan pl_w = NULL;
443  fftw_complex *x_cx1, *x_cx2;
444 
445  if (R == 0 || C == 0) // 1D
446  {
447  x_cx1 = fftw_alloc_complex(N);
448  x_cx2 = fftw_alloc_complex(N);
449 
450  // Replan using wisdom
451  pl_w = fftw_plan_dft_1d(N, x_cx1, x_cx2, FFTW_FORWARD, export_alg);
452  }
453  else
454  {
455  x_cx1 = fftw_alloc_complex(R * C);
456  x_cx2 = fftw_alloc_complex(R * C);
457 
458  // Replan using wisdom
459  pl_w = fftw_plan_dft_2d(C, R, x_cx1, x_cx2, FFTW_FORWARD, export_alg);
460  }
461 
462  if (pl_w == NULL)
463  {
464  err_handler("Unable to create complex data FFTW plan");
465  }
466 
467  // Export
468  if (fftw_export_wisdom_to_filename(fname.c_str()) == 0)
469  {
470  err_handler("Could not export wisdom to file!");
471  }
472 
473  fftw_destroy_plan(pl_w);
474  fftw_free(x_cx1);
475  fftw_free(x_cx2);
476  }
481  void export_wisdom_ifft_cx(const std::string fname)
482  {
483  fftw_plan pl_w = NULL;
484  fftw_complex *x_cx1, *x_cx2;
485 
486  if (R == 0 || C == 0) // 1D
487  {
488  x_cx1 = fftw_alloc_complex(N);
489  x_cx2 = fftw_alloc_complex(N);
490 
491  // Replan using wisdom
492  pl_w = fftw_plan_dft_1d(N, x_cx2, x_cx1, FFTW_BACKWARD, export_alg);
493  }
494  else
495  {
496  x_cx1 = fftw_alloc_complex(R * C);
497  x_cx2 = fftw_alloc_complex(R * C);
498 
499  // Replan using wisdom
500  pl_w = fftw_plan_dft_2d(C, R, x_cx2, x_cx1, FFTW_BACKWARD, export_alg);
501  }
502  if (pl_w == NULL)
503  {
504  err_handler("Unable to create complex data IFFTW plan");
505  }
506 
507  // Export
508  if (fftw_export_wisdom_to_filename(fname.c_str()) == 0)
509  {
510  err_handler("Could not export wisdom to file!");
511  }
512 
513  fftw_destroy_plan(pl_w);
514  fftw_free(x_cx1);
515  fftw_free(x_cx2);
516  }
517 };
519 
520 } // namespace sp
521 #endif
void import_wisdom_file(const std::string fname)
Import wisdom from file.
Definition: fftw.h:344
unsigned int C
FFT 2D dims.
Definition: fftw.h:33
arma::mat ifft2(arma::cx_mat &Pxx)
Inverse FFT.
Definition: fftw.h:322
arma::cx_vec fft(arma::vec &x)
FFT of real input.
Definition: fftw.h:188
#define err_handler(msg)
Definition: base.h:213
void ifft2(arma::cx_mat &Pxx, arma::mat &x)
Inverse 2D FFT.
Definition: fftw.h:288
void export_wisdom_fft_cx(const std::string fname)
Export complex FFT wisdom to file.
Definition: fftw.h:440
Definition: base.h:7
void import_wisdom_string(const std::string wisd)
Import wisdom from string.
Definition: fftw.h:333
fftw_plan pl_ifft_cx
Complex IFFTW plan.
Definition: fftw.h:31
fftw_plan pl_fft
Real FFTW plan.
Definition: fftw.h:28
~FFTW()
Destructor.
Definition: fftw.h:79
FFTW(unsigned int _R, unsigned int _C, int _alg)
Constructor.
Definition: fftw.h:63
arma::cx_vec ifft_cx(arma::cx_vec &Pxx)
Inverse FFT.
Definition: fftw.h:150
unsigned int R
Definition: fftw.h:33
FFTW(unsigned int _N, int _alg=FFTW_ESTIMATE)
Constructor.
Definition: fftw.h:44
unsigned int N
FFT length.
Definition: fftw.h:32
arma::cx_mat fft2(arma::mat &x)
FFT of real 2D input.
Definition: fftw.h:276
void ifft(arma::cx_vec &Pxx, arma::vec &x)
Inverse FFT.
Definition: fftw.h:200
void export_wisdom_fft(const std::string fname)
Export real FFT wisdom to file.
Definition: fftw.h:355
arma::vec ifft(arma::cx_vec &Pxx)
Inverse FFT.
Definition: fftw.h:221
void fft2(arma::mat &x, arma::cx_mat &Pxx)
FFT of real 2D input.
Definition: fftw.h:233
fftw_plan pl_ifft
Real IFFTW plan.
Definition: fftw.h:29
arma::cx_vec fft_cx(arma::cx_vec &x)
FFT of complex input.
Definition: fftw.h:117
int export_alg
Alg used for exporting wisdom.
Definition: fftw.h:37
FFTW class.
Definition: fftw.h:25
void export_wisdom_ifft(const std::string fname)
Export real IFFT wisdom to file.
Definition: fftw.h:397
void ifft_cx(arma::cx_vec &Pxx, arma::cx_vec &x)
Inverse FFT.
Definition: fftw.h:129
void export_wisdom_ifft_cx(const std::string fname)
Export complex IFFT wisdom to file.
Definition: fftw.h:481
int alg
plans](http://fftw.org/fftw3_doc/Planner-Flags.html#Planner-Flags)
Definition: fftw.h:34
void fft_cx(arma::cx_vec &x, arma::cx_vec &Pxx)
FFT of complex input.
Definition: fftw.h:97
void fft(arma::vec &x, arma::cx_vec &Pxx)
FFT of real input.
Definition: fftw.h:162
fftw_plan pl_fft_cx
Complex FFTW plan.
Definition: fftw.h:30