IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
transforms.cpp
Go to the documentation of this file.
1 
30 #ifndef _MSC_VER
31 # include <itpp/config.h>
32 #else
33 # include <itpp/config_msvc.h>
34 #endif
35 
36 #if defined(HAVE_FFT_MKL)
37 namespace mkl
38 {
39 # include <mkl_dfti.h>
40 # undef DftiCreateDescriptor
41 }
42 #elif defined(HAVE_FFT_ACML)
43 namespace acml
44 {
45 # include <acml.h>
46 }
47 #elif defined(HAVE_FFTW3)
48 # include <fftw3.h>
49 #endif
50 
51 #include <itpp/signal/transforms.h>
52 
54 
55 namespace itpp
56 {
57 
58 #if defined(HAVE_FFT_MKL)
59 
60 //---------------------------------------------------------------------------
61 // FFT/IFFT based on MKL
62 //---------------------------------------------------------------------------
63 
64 void fft(const cvec &in, cvec &out)
65 {
66  static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
67  static int N;
68 
69  out.set_size(in.size(), false);
70  if (N != in.size()) {
71  N = in.size();
72  if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle);
73  mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_COMPLEX, 1, N);
74  mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE);
75  mkl::DftiCommitDescriptor(fft_handle);
76  }
77  mkl::DftiComputeForward(fft_handle, (void *)in._data(), out._data());
78 }
79 
80 void ifft(const cvec &in, cvec &out)
81 {
82  static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
83  static int N;
84 
85  out.set_size(in.size(), false);
86  if (N != in.size()) {
87  N = in.size();
88  if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle);
89  mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_COMPLEX, 1, N);
90  mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE);
91  mkl::DftiSetValue(fft_handle, mkl::DFTI_BACKWARD_SCALE, 1.0 / N);
92  mkl::DftiCommitDescriptor(fft_handle);
93  }
94  mkl::DftiComputeBackward(fft_handle, (void *)in._data(), out._data());
95 }
96 
97 void fft_real(const vec &in, cvec &out)
98 {
99  static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
100  static int N;
101 
102  out.set_size(in.size(), false);
103  if (N != in.size()) {
104  N = in.size();
105  if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle);
106  mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_REAL, 1, N);
107  mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE);
108  mkl::DftiCommitDescriptor(fft_handle);
109  }
110  mkl::DftiComputeForward(fft_handle, (void *)in._data(), out._data());
111 
112  // Real FFT does not compute the 2nd half of the FFT points because it
113  // is redundant to the 1st half. However, we want all of the data so we
114  // fill it in. This is consistent with Matlab's functionality
115  int istart = ceil_i(in.size() / 2.0);
116  int idelta = in.size() - istart;
117  out.set_subvector(istart, reverse(conj(out(1, idelta))));
118 }
119 
120 void ifft_real(const cvec &in, vec &out)
121 {
122  static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
123  static int N;
124 
125  out.set_size(in.size(), false);
126  if (N != in.size()) {
127  N = in.size();
128  if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle);
129  mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_REAL, 1, N);
130  mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE);
131  mkl::DftiSetValue(fft_handle, mkl::DFTI_BACKWARD_SCALE, 1.0 / N);
132  mkl::DftiCommitDescriptor(fft_handle);
133  }
134  mkl::DftiComputeBackward(fft_handle, (void *)in._data(), out._data());
135 }
136 
137 #endif // #ifdef HAVE_FFT_MKL
138 
139 
140 #if defined(HAVE_FFT_ACML)
141 
142 //---------------------------------------------------------------------------
143 // FFT/IFFT based on ACML
144 //---------------------------------------------------------------------------
145 
146 void fft(const cvec &in, cvec &out)
147 {
148  static int N = 0;
149  static cvec *comm_ptr = NULL;
150  int info;
151  out.set_size(in.size(), false);
152 
153  if (N != in.size()) {
154  N = in.size();
155  if (comm_ptr != NULL)
156  delete comm_ptr;
157  comm_ptr = new cvec(5 * in.size() + 100);
158  acml::zfft1dx(0, 1.0, false, N, (acml::doublecomplex *)in._data(), 1,
159  (acml::doublecomplex *)out._data(), 1,
160  (acml::doublecomplex *)comm_ptr->_data(), &info);
161  }
162  acml::zfft1dx(-1, 1.0, false, N, (acml::doublecomplex *)in._data(), 1,
163  (acml::doublecomplex *)out._data(), 1,
164  (acml::doublecomplex *)comm_ptr->_data(), &info);
165 }
166 
167 void ifft(const cvec &in, cvec &out)
168 {
169  static int N = 0;
170  static cvec *comm_ptr = NULL;
171  int info;
172  out.set_size(in.size(), false);
173 
174  if (N != in.size()) {
175  N = in.size();
176  if (comm_ptr != NULL)
177  delete comm_ptr;
178  comm_ptr = new cvec(5 * in.size() + 100);
179  acml::zfft1dx(0, 1.0 / N, false, N, (acml::doublecomplex *)in._data(), 1,
180  (acml::doublecomplex *)out._data(), 1,
181  (acml::doublecomplex *)comm_ptr->_data(), &info);
182  }
183  acml::zfft1dx(1, 1.0 / N, false, N, (acml::doublecomplex *)in._data(), 1,
184  (acml::doublecomplex *)out._data(), 1,
185  (acml::doublecomplex *)comm_ptr->_data(), &info);
186 }
187 
188 void fft_real(const vec &in, cvec &out)
189 {
190  static int N = 0;
191  static double factor = 0;
192  static vec *comm_ptr = NULL;
193  int info;
194  vec out_re = in;
195 
196  if (N != in.size()) {
197  N = in.size();
198  factor = std::sqrt(static_cast<double>(N));
199  if (comm_ptr != NULL)
200  delete comm_ptr;
201  comm_ptr = new vec(5 * in.size() + 100);
202  acml::dzfft(0, N, out_re._data(), comm_ptr->_data(), &info);
203  }
204  acml::dzfft(1, N, out_re._data(), comm_ptr->_data(), &info);
205 
206  // Normalise output data
207  out_re *= factor;
208 
209  // Convert the real Hermitian DZFFT's output to the Matlab's complex form
210  vec out_im(in.size());
211  out_im.zeros();
212  out.set_size(in.size(), false);
213  out_im.set_subvector(1, reverse(out_re(N / 2 + 1, N - 1)));
214  out_im.set_subvector(N / 2 + 1, -out_re(N / 2 + 1, N - 1));
215  out_re.set_subvector(N / 2 + 1, reverse(out_re(1, (N - 1) / 2)));
216  out = to_cvec(out_re, out_im);
217 }
218 
219 void ifft_real(const cvec &in, vec &out)
220 {
221  static int N = 0;
222  static double factor = 0;
223  static vec *comm_ptr = NULL;
224  int info;
225 
226  // Convert Matlab's complex input to the real Hermitian form
227  out.set_size(in.size());
228  out.set_subvector(0, real(in(0, in.size() / 2)));
229  out.set_subvector(in.size() / 2 + 1, -imag(in(in.size() / 2 + 1, in.size() - 1)));
230 
231  if (N != in.size()) {
232  N = in.size();
233  factor = 1.0 / std::sqrt(static_cast<double>(N));
234  if (comm_ptr != NULL)
235  delete comm_ptr;
236  comm_ptr = new vec(5 * in.size() + 100);
237  acml::zdfft(0, N, out._data(), comm_ptr->_data(), &info);
238  }
239  acml::zdfft(1, N, out._data(), comm_ptr->_data(), &info);
240  out.set_subvector(1, reverse(out(1, N - 1)));
241 
242  // Normalise output data
243  out *= factor;
244 }
245 
246 #endif // defined(HAVE_FFT_ACML)
247 
248 
249 #if defined(HAVE_FFTW3)
250 
251 //---------------------------------------------------------------------------
252 // FFT/IFFT based on FFTW
253 //---------------------------------------------------------------------------
254 
255 void fft(const cvec &in, cvec &out)
256 {
257  static int N = 0;
258  static fftw_plan p = NULL;
259  out.set_size(in.size(), false);
260 
261  if (N != in.size()) {
262  N = in.size();
263  if (p != NULL)
264  fftw_destroy_plan(p); // destroy the previous plan
265  // create a new plan
266  p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(),
267  (fftw_complex *)out._data(),
268  FFTW_FORWARD, FFTW_ESTIMATE);
269  }
270 
271  // compute FFT using the GURU FFTW interface
272  fftw_execute_dft(p, (fftw_complex *)in._data(),
273  (fftw_complex *)out._data());
274 }
275 
276 void ifft(const cvec &in, cvec &out)
277 {
278  static int N = 0;
279  static double inv_N;
280  static fftw_plan p = NULL;
281  out.set_size(in.size(), false);
282 
283  if (N != in.size()) {
284  N = in.size();
285  inv_N = 1.0 / N;
286  if (p != NULL)
287  fftw_destroy_plan(p); // destroy the previous plan
288  // create a new plan
289  p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(),
290  (fftw_complex *)out._data(),
291  FFTW_BACKWARD, FFTW_ESTIMATE);
292  }
293 
294  // compute IFFT using the GURU FFTW interface
295  fftw_execute_dft(p, (fftw_complex *)in._data(),
296  (fftw_complex *)out._data());
297 
298  // scale output
299  out *= inv_N;
300 }
301 
302 void fft_real(const vec &in, cvec &out)
303 {
304  static int N = 0;
305  static fftw_plan p = NULL;
306  out.set_size(in.size(), false);
307 
308  if (N != in.size()) {
309  N = in.size();
310  if (p != NULL)
311  fftw_destroy_plan(p); //destroy the previous plan
312 
313  // create a new plan
314  p = fftw_plan_dft_r2c_1d(N, (double *)in._data(),
315  (fftw_complex *)out._data(),
316  FFTW_ESTIMATE);
317  }
318 
319  // compute FFT using the GURU FFTW interface
320  fftw_execute_dft_r2c(p, (double *)in._data(),
321  (fftw_complex *)out._data());
322 
323  // Real FFT does not compute the 2nd half of the FFT points because it
324  // is redundant to the 1st half. However, we want all of the data so we
325  // fill it in. This is consistent with Matlab's functionality
326  int offset = ceil_i(in.size() / 2.0);
327  int n_elem = in.size() - offset;
328  for (int i = 0; i < n_elem; ++i) {
329  out(offset + i) = std::conj(out(n_elem - i));
330  }
331 }
332 
333 void ifft_real(const cvec &in, vec & out)
334 {
335  static int N = 0;
336  static double inv_N;
337  static fftw_plan p = NULL;
338  out.set_size(in.size(), false);
339 
340  if (N != in.size()) {
341  N = in.size();
342  inv_N = 1.0 / N;
343  if (p != NULL)
344  fftw_destroy_plan(p); // destroy the previous plan
345 
346  // create a new plan
347  p = fftw_plan_dft_c2r_1d(N, (fftw_complex *)in._data(),
348  (double *)out._data(),
349  FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
350  }
351 
352  // compute IFFT using the GURU FFTW interface
353  fftw_execute_dft_c2r(p, (fftw_complex *)in._data(),
354  (double *)out._data());
355 
356  out *= inv_N;
357 }
358 
359 //---------------------------------------------------------------------------
360 // DCT/IDCT based on FFTW
361 //---------------------------------------------------------------------------
362 
363 void dct(const vec &in, vec &out)
364 {
365  static int N;
366  static fftw_plan p = NULL;
367  out.set_size(in.size(), false);
368 
369  if (N != in.size()) {
370  N = in.size();
371  if (p != NULL)
372  fftw_destroy_plan(p); // destroy the previous plan
373 
374  // create a new plan
375  p = fftw_plan_r2r_1d(N, (double *)in._data(),
376  (double *)out._data(),
377  FFTW_REDFT10, FFTW_ESTIMATE);
378  }
379 
380  // compute FFT using the GURU FFTW interface
381  fftw_execute_r2r(p, (double *)in._data(), (double *)out._data());
382 
383  // Scale to matlab definition format
384  out /= std::sqrt(2.0 * N);
385  out(0) /= std::sqrt(2.0);
386 }
387 
388 // IDCT
389 void idct(const vec &in, vec &out)
390 {
391  static int N;
392  static fftw_plan p = NULL;
393  out = in;
394 
395  // Rescale to FFTW format
396  out(0) *= std::sqrt(2.0);
397  out /= std::sqrt(2.0 * in.size());
398 
399  if (N != in.size()) {
400  N = in.size();
401  if (p != NULL)
402  fftw_destroy_plan(p); // destroy the previous plan
403 
404  // create a new plan
405  p = fftw_plan_r2r_1d(N, (double *)out._data(),
406  (double *)out._data(),
407  FFTW_REDFT01, FFTW_ESTIMATE);
408  }
409 
410  // compute FFT using the GURU FFTW interface
411  fftw_execute_r2r(p, (double *)out._data(), (double *)out._data());
412 }
413 
414 #endif // defined(HAVE_FFTW3)
415 
416 
417 #if defined(HAVE_FFT_MKL) || defined(HAVE_FFT_ACML)
418 
419 //---------------------------------------------------------------------------
420 // DCT/IDCT based on MKL or ACML
421 //---------------------------------------------------------------------------
422 
423 void dct(const vec &in, vec &out)
424 {
425  int N = in.size();
426  if (N == 1)
427  out = in;
428  else {
429  cvec c = fft_real(concat(in, reverse(in)));
430  c.set_size(N, true);
431  for (int i = 0; i < N; i++) {
432  c(i) *= std::complex<double>(std::cos(pi * i / N / 2), std::sin(-pi * i / N / 2))
433  / std::sqrt(2.0 * N);
434  }
435  out = real(c);
436  out(0) /= std::sqrt(2.0);
437  }
438 }
439 
440 void idct(const vec &in, vec &out)
441 {
442  int N = in.size();
443  if (N == 1)
444  out = in;
445  else {
446  cvec c = to_cvec(in);
447  c.set_size(2*N, true);
448  c(0) *= std::sqrt(2.0);
449  for (int i = 0; i < N; i++) {
450  c(i) *= std::complex<double>(std::cos(pi * i / N / 2), std::sin(pi * i / N / 2))
451  * std::sqrt(2.0 * N);
452  }
453  for (int i = N - 1; i >= 1; i--) {
454  c(c.size() - i) = c(i) * std::complex<double>(std::cos(pi * i / N),
455  std::sin(-pi * i / N));
456  }
457  out = ifft_real(c);
458  out.set_size(N, true);
459  }
460 }
461 
462 #endif // defined(HAVE_FFT_MKL) || defined(HAVE_FFT_ACML)
463 
464 
465 #if !defined(HAVE_FFT)
466 
467 void fft(const cvec &in, cvec &out)
468 {
469  it_error("FFT library is needed to use fft() function");
470 }
471 
472 void ifft(const cvec &in, cvec &out)
473 {
474  it_error("FFT library is needed to use ifft() function");
475 }
476 
477 void fft_real(const vec &in, cvec &out)
478 {
479  it_error("FFT library is needed to use fft_real() function");
480 }
481 
482 void ifft_real(const cvec &in, vec & out)
483 {
484  it_error("FFT library is needed to use ifft_real() function");
485 }
486 
487 void dct(const vec &in, vec &out)
488 {
489  it_error("FFT library is needed to use dct() function");
490 }
491 
492 void idct(const vec &in, vec &out)
493 {
494  it_error("FFT library is needed to use idct() function");
495 }
496 
497 #endif // !defined(HAVE_FFT)
498 
499 cvec fft(const cvec &in)
500 {
501  cvec out;
502  fft(in, out);
503  return out;
504 }
505 
506 cvec fft(const cvec &in, const int N)
507 {
508  cvec in2 = in;
509  cvec out;
510  in2.set_size(N, true);
511  fft(in2, out);
512  return out;
513 }
514 
515 cvec ifft(const cvec &in)
516 {
517  cvec out;
518  ifft(in, out);
519  return out;
520 }
521 
522 cvec ifft(const cvec &in, const int N)
523 {
524  cvec in2 = in;
525  cvec out;
526  in2.set_size(N, true);
527  ifft(in2, out);
528  return out;
529 }
530 
531 cvec fft_real(const vec& in)
532 {
533  cvec out;
534  fft_real(in, out);
535  return out;
536 }
537 
538 cvec fft_real(const vec& in, const int N)
539 {
540  vec in2 = in;
541  cvec out;
542  in2.set_size(N, true);
543  fft_real(in2, out);
544  return out;
545 }
546 
547 vec ifft_real(const cvec &in)
548 {
549  vec out;
550  ifft_real(in, out);
551  return out;
552 }
553 
554 vec ifft_real(const cvec &in, const int N)
555 {
556  cvec in2 = in;
557  in2.set_size(N, true);
558  vec out;
559  ifft_real(in2, out);
560  return out;
561 }
562 
563 vec dct(const vec &in)
564 {
565  vec out;
566  dct(in, out);
567  return out;
568 }
569 
570 vec idct(const vec &in)
571 {
572  vec out;
573  idct(in, out);
574  return out;
575 }
576 
577 
578 // ----------------------------------------------------------------------
579 // Instantiation
580 // ----------------------------------------------------------------------
581 
582 template vec dht(const vec &v);
583 template cvec dht(const cvec &v);
584 
585 template void dht(const vec &vin, vec &vout);
586 template void dht(const cvec &vin, cvec &vout);
587 
588 template void self_dht(vec &v);
589 template void self_dht(cvec &v);
590 
591 template vec dwht(const vec &v);
592 template cvec dwht(const cvec &v);
593 
594 template void dwht(const vec &vin, vec &vout);
595 template void dwht(const cvec &vin, cvec &vout);
596 
597 template void self_dwht(vec &v);
598 template void self_dwht(cvec &v);
599 
600 template mat dht2(const mat &m);
601 template cmat dht2(const cmat &m);
602 
603 template mat dwht2(const mat &m);
604 template cmat dwht2(const cmat &m);
605 
606 } // namespace itpp
607 
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:14 for IT++ by Doxygen 1.8.1.2