31 # include <itpp/config.h>
33 # include <itpp/config_msvc.h>
36 #if defined(HAVE_FFT_MKL)
39 # include <mkl_dfti.h>
40 # undef DftiCreateDescriptor
42 #elif defined(HAVE_FFT_ACML)
47 #elif defined(HAVE_FFTW3)
58 #if defined(HAVE_FFT_MKL)
64 void fft(
const cvec &in, cvec &out)
66 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
69 out.set_size(in.size(),
false);
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);
77 mkl::DftiComputeForward(fft_handle, (
void *)in._data(), out._data());
80 void ifft(
const cvec &in, cvec &out)
82 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
85 out.set_size(in.size(),
false);
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);
94 mkl::DftiComputeBackward(fft_handle, (
void *)in._data(), out._data());
97 void fft_real(
const vec &in, cvec &out)
99 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
102 out.set_size(in.size(),
false);
103 if (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);
110 mkl::DftiComputeForward(fft_handle, (
void *)in._data(), out._data());
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))));
122 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL;
125 out.set_size(in.size(),
false);
126 if (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);
134 mkl::DftiComputeBackward(fft_handle, (
void *)in._data(), out._data());
137 #endif // #ifdef HAVE_FFT_MKL
140 #if defined(HAVE_FFT_ACML)
146 void fft(
const cvec &in, cvec &out)
149 static cvec *comm_ptr = NULL;
151 out.set_size(in.size(),
false);
153 if (N != in.size()) {
155 if (comm_ptr != NULL)
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);
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);
167 void ifft(
const cvec &in, cvec &out)
170 static cvec *comm_ptr = NULL;
172 out.set_size(in.size(),
false);
174 if (N != in.size()) {
176 if (comm_ptr != NULL)
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);
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);
188 void fft_real(
const vec &in, cvec &out)
191 static double factor = 0;
192 static vec *comm_ptr = NULL;
196 if (N != in.size()) {
198 factor =
std::sqrt(static_cast<double>(N));
199 if (comm_ptr != NULL)
201 comm_ptr =
new vec(5 * in.size() + 100);
202 acml::dzfft(0, N, out_re._data(), comm_ptr->_data(), &info);
204 acml::dzfft(1, N, out_re._data(), comm_ptr->_data(), &info);
210 vec out_im(in.size());
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)));
222 static double factor = 0;
223 static vec *comm_ptr = NULL;
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)));
231 if (N != in.size()) {
233 factor = 1.0 /
std::sqrt(static_cast<double>(N));
234 if (comm_ptr != NULL)
236 comm_ptr =
new vec(5 * in.size() + 100);
237 acml::zdfft(0, N, out._data(), comm_ptr->_data(), &info);
239 acml::zdfft(1, N, out._data(), comm_ptr->_data(), &info);
240 out.set_subvector(1,
reverse(out(1, N - 1)));
246 #endif // defined(HAVE_FFT_ACML)
249 #if defined(HAVE_FFTW3)
255 void fft(
const cvec &in, cvec &out)
258 static fftw_plan p = NULL;
259 out.set_size(in.size(),
false);
261 if (N != in.size()) {
264 fftw_destroy_plan(p);
266 p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(),
267 (fftw_complex *)out._data(),
268 FFTW_FORWARD, FFTW_ESTIMATE);
272 fftw_execute_dft(p, (fftw_complex *)in._data(),
273 (fftw_complex *)out._data());
276 void ifft(
const cvec &in, cvec &out)
280 static fftw_plan p = NULL;
281 out.set_size(in.size(),
false);
283 if (N != in.size()) {
287 fftw_destroy_plan(p);
289 p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(),
290 (fftw_complex *)out._data(),
291 FFTW_BACKWARD, FFTW_ESTIMATE);
295 fftw_execute_dft(p, (fftw_complex *)in._data(),
296 (fftw_complex *)out._data());
302 void fft_real(
const vec &in, cvec &out)
305 static fftw_plan p = NULL;
306 out.set_size(in.size(),
false);
308 if (N != in.size()) {
311 fftw_destroy_plan(p);
314 p = fftw_plan_dft_r2c_1d(N, (
double *)in._data(),
315 (fftw_complex *)out._data(),
320 fftw_execute_dft_r2c(p, (
double *)in._data(),
321 (fftw_complex *)out._data());
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));
333 void ifft_real(
const cvec &in, vec & out)
337 static fftw_plan p = NULL;
338 out.set_size(in.size(),
false);
340 if (N != in.size()) {
344 fftw_destroy_plan(p);
347 p = fftw_plan_dft_c2r_1d(N, (fftw_complex *)in._data(),
348 (
double *)out._data(),
349 FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
353 fftw_execute_dft_c2r(p, (fftw_complex *)in._data(),
354 (
double *)out._data());
363 void dct(
const vec &in, vec &out)
366 static fftw_plan p = NULL;
367 out.set_size(in.size(),
false);
369 if (N != in.size()) {
372 fftw_destroy_plan(p);
375 p = fftw_plan_r2r_1d(N, (
double *)in._data(),
376 (
double *)out._data(),
377 FFTW_REDFT10, FFTW_ESTIMATE);
381 fftw_execute_r2r(p, (
double *)in._data(), (
double *)out._data());
389 void idct(
const vec &in, vec &out)
392 static fftw_plan p = NULL;
399 if (N != in.size()) {
402 fftw_destroy_plan(p);
405 p = fftw_plan_r2r_1d(N, (
double *)out._data(),
406 (
double *)out._data(),
407 FFTW_REDFT01, FFTW_ESTIMATE);
411 fftw_execute_r2r(p, (
double *)out._data(), (
double *)out._data());
414 #endif // defined(HAVE_FFTW3)
417 #if defined(HAVE_FFT_MKL) || defined(HAVE_FFT_ACML)
423 void dct(
const vec &in, vec &out)
431 for (
int i = 0; i < N; i++) {
440 void idct(
const vec &in, vec &out)
447 c.set_size(2*N,
true);
449 for (
int i = 0; i < N; i++) {
453 for (
int i = N - 1; i >= 1; i--) {
454 c(c.size() - i) = c(i) * std::complex<double>(
std::cos(
pi * i / N),
458 out.set_size(N,
true);
462 #endif // defined(HAVE_FFT_MKL) || defined(HAVE_FFT_ACML)
465 #if !defined(HAVE_FFT)
467 void fft(
const cvec &in, cvec &out)
469 it_error(
"FFT library is needed to use fft() function");
472 void ifft(
const cvec &in, cvec &out)
474 it_error(
"FFT library is needed to use ifft() function");
477 void fft_real(
const vec &in, cvec &out)
479 it_error(
"FFT library is needed to use fft_real() function");
482 void ifft_real(
const cvec &in, vec & out)
484 it_error(
"FFT library is needed to use ifft_real() function");
487 void dct(
const vec &in, vec &out)
489 it_error(
"FFT library is needed to use dct() function");
492 void idct(
const vec &in, vec &out)
494 it_error(
"FFT library is needed to use idct() function");
497 #endif // !defined(HAVE_FFT)
499 cvec
fft(
const cvec &in)
506 cvec
fft(
const cvec &in,
const int N)
510 in2.set_size(N,
true);
515 cvec
ifft(
const cvec &in)
522 cvec
ifft(
const cvec &in,
const int N)
526 in2.set_size(N,
true);
538 cvec
fft_real(
const vec& in,
const int N)
542 in2.set_size(N,
true);
554 vec
ifft_real(
const cvec &in,
const int N)
557 in2.set_size(N,
true);
563 vec
dct(
const vec &in)
570 vec
idct(
const vec &in)
582 template vec
dht(
const vec &v);
583 template cvec
dht(
const cvec &v);
585 template void dht(
const vec &vin, vec &vout);
586 template void dht(
const cvec &vin, cvec &vout);
591 template vec
dwht(
const vec &v);
592 template cvec
dwht(
const cvec &v);
594 template void dwht(
const vec &vin, vec &vout);
595 template void dwht(
const cvec &vin, cvec &vout);
600 template mat
dht2(
const mat &m);
601 template cmat
dht2(
const cmat &m);
603 template mat
dwht2(
const mat &m);
604 template cmat
dwht2(
const cmat &m);