00001 00030 #ifndef _MSC_VER 00031 # include <itpp/config.h> 00032 #else 00033 # include <itpp/config_msvc.h> 00034 #endif 00035 00036 #if defined(HAVE_FFT_MKL) 00037 namespace mkl 00038 { 00039 # include <mkl_dfti.h> 00040 # undef DftiCreateDescriptor 00041 } 00042 #elif defined(HAVE_FFT_ACML) 00043 namespace acml 00044 { 00045 # include <acml.h> 00046 } 00047 #elif defined(HAVE_FFTW3) 00048 # include <fftw3.h> 00049 #endif 00050 00051 #include <itpp/signal/transforms.h> 00052 00054 00055 namespace itpp 00056 { 00057 00058 #if defined(HAVE_FFT_MKL) 00059 00060 //--------------------------------------------------------------------------- 00061 // FFT/IFFT based on MKL 00062 //--------------------------------------------------------------------------- 00063 00064 void fft(const cvec &in, cvec &out) 00065 { 00066 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL; 00067 static int N; 00068 00069 out.set_size(in.size(), false); 00070 if (N != in.size()) { 00071 N = in.size(); 00072 if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle); 00073 mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_COMPLEX, 1, N); 00074 mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE); 00075 mkl::DftiCommitDescriptor(fft_handle); 00076 } 00077 mkl::DftiComputeForward(fft_handle, (void *)in._data(), out._data()); 00078 } 00079 00080 void ifft(const cvec &in, cvec &out) 00081 { 00082 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL; 00083 static int N; 00084 00085 out.set_size(in.size(), false); 00086 if (N != in.size()) { 00087 N = in.size(); 00088 if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle); 00089 mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_COMPLEX, 1, N); 00090 mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE); 00091 mkl::DftiSetValue(fft_handle, mkl::DFTI_BACKWARD_SCALE, 1.0 / N); 00092 mkl::DftiCommitDescriptor(fft_handle); 00093 } 00094 mkl::DftiComputeBackward(fft_handle, (void *)in._data(), out._data()); 00095 } 00096 00097 void fft_real(const vec &in, cvec &out) 00098 { 00099 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL; 00100 static int N; 00101 00102 out.set_size(in.size(), false); 00103 if (N != in.size()) { 00104 N = in.size(); 00105 if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle); 00106 mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_REAL, 1, N); 00107 mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE); 00108 mkl::DftiCommitDescriptor(fft_handle); 00109 } 00110 mkl::DftiComputeForward(fft_handle, (void *)in._data(), out._data()); 00111 00112 // Real FFT does not compute the 2nd half of the FFT points because it 00113 // is redundant to the 1st half. However, we want all of the data so we 00114 // fill it in. This is consistent with Matlab's functionality 00115 int istart = ceil_i(in.size() / 2.0); 00116 int idelta = in.size() - istart; 00117 out.set_subvector(istart, reverse(conj(out(1, idelta)))); 00118 } 00119 00120 void ifft_real(const cvec &in, vec &out) 00121 { 00122 static mkl::DFTI_DESCRIPTOR* fft_handle = NULL; 00123 static int N; 00124 00125 out.set_size(in.size(), false); 00126 if (N != in.size()) { 00127 N = in.size(); 00128 if (fft_handle != NULL) mkl::DftiFreeDescriptor(&fft_handle); 00129 mkl::DftiCreateDescriptor(&fft_handle, mkl::DFTI_DOUBLE, mkl::DFTI_REAL, 1, N); 00130 mkl::DftiSetValue(fft_handle, mkl::DFTI_PLACEMENT, mkl::DFTI_NOT_INPLACE); 00131 mkl::DftiSetValue(fft_handle, mkl::DFTI_BACKWARD_SCALE, 1.0 / N); 00132 mkl::DftiCommitDescriptor(fft_handle); 00133 } 00134 mkl::DftiComputeBackward(fft_handle, (void *)in._data(), out._data()); 00135 } 00136 00137 #endif // #ifdef HAVE_FFT_MKL 00138 00139 00140 #if defined(HAVE_FFT_ACML) 00141 00142 //--------------------------------------------------------------------------- 00143 // FFT/IFFT based on ACML 00144 //--------------------------------------------------------------------------- 00145 00146 void fft(const cvec &in, cvec &out) 00147 { 00148 static int N = 0; 00149 static cvec *comm_ptr = NULL; 00150 int info; 00151 out.set_size(in.size(), false); 00152 00153 if (N != in.size()) { 00154 N = in.size(); 00155 if (comm_ptr != NULL) 00156 delete comm_ptr; 00157 comm_ptr = new cvec(5 * in.size() + 100); 00158 acml::zfft1dx(0, 1.0, false, N, (acml::doublecomplex *)in._data(), 1, 00159 (acml::doublecomplex *)out._data(), 1, 00160 (acml::doublecomplex *)comm_ptr->_data(), &info); 00161 } 00162 acml::zfft1dx(-1, 1.0, false, N, (acml::doublecomplex *)in._data(), 1, 00163 (acml::doublecomplex *)out._data(), 1, 00164 (acml::doublecomplex *)comm_ptr->_data(), &info); 00165 } 00166 00167 void ifft(const cvec &in, cvec &out) 00168 { 00169 static int N = 0; 00170 static cvec *comm_ptr = NULL; 00171 int info; 00172 out.set_size(in.size(), false); 00173 00174 if (N != in.size()) { 00175 N = in.size(); 00176 if (comm_ptr != NULL) 00177 delete comm_ptr; 00178 comm_ptr = new cvec(5 * in.size() + 100); 00179 acml::zfft1dx(0, 1.0 / N, false, N, (acml::doublecomplex *)in._data(), 1, 00180 (acml::doublecomplex *)out._data(), 1, 00181 (acml::doublecomplex *)comm_ptr->_data(), &info); 00182 } 00183 acml::zfft1dx(1, 1.0 / N, false, N, (acml::doublecomplex *)in._data(), 1, 00184 (acml::doublecomplex *)out._data(), 1, 00185 (acml::doublecomplex *)comm_ptr->_data(), &info); 00186 } 00187 00188 void fft_real(const vec &in, cvec &out) 00189 { 00190 static int N = 0; 00191 static double factor = 0; 00192 static vec *comm_ptr = NULL; 00193 int info; 00194 vec out_re = in; 00195 00196 if (N != in.size()) { 00197 N = in.size(); 00198 factor = std::sqrt(static_cast<double>(N)); 00199 if (comm_ptr != NULL) 00200 delete comm_ptr; 00201 comm_ptr = new vec(5 * in.size() + 100); 00202 acml::dzfft(0, N, out_re._data(), comm_ptr->_data(), &info); 00203 } 00204 acml::dzfft(1, N, out_re._data(), comm_ptr->_data(), &info); 00205 00206 // Normalise output data 00207 out_re *= factor; 00208 00209 // Convert the real Hermitian DZFFT's output to the Matlab's complex form 00210 vec out_im(in.size()); 00211 out_im.zeros(); 00212 out.set_size(in.size(), false); 00213 out_im.set_subvector(1, reverse(out_re(N / 2 + 1, N - 1))); 00214 out_im.set_subvector(N / 2 + 1, -out_re(N / 2 + 1, N - 1)); 00215 out_re.set_subvector(N / 2 + 1, reverse(out_re(1, (N - 1) / 2))); 00216 out = to_cvec(out_re, out_im); 00217 } 00218 00219 void ifft_real(const cvec &in, vec &out) 00220 { 00221 static int N = 0; 00222 static double factor = 0; 00223 static vec *comm_ptr = NULL; 00224 int info; 00225 00226 // Convert Matlab's complex input to the real Hermitian form 00227 out.set_size(in.size()); 00228 out.set_subvector(0, real(in(0, in.size() / 2))); 00229 out.set_subvector(in.size() / 2 + 1, -imag(in(in.size() / 2 + 1, in.size() - 1))); 00230 00231 if (N != in.size()) { 00232 N = in.size(); 00233 factor = 1.0 / std::sqrt(static_cast<double>(N)); 00234 if (comm_ptr != NULL) 00235 delete comm_ptr; 00236 comm_ptr = new vec(5 * in.size() + 100); 00237 acml::zdfft(0, N, out._data(), comm_ptr->_data(), &info); 00238 } 00239 acml::zdfft(1, N, out._data(), comm_ptr->_data(), &info); 00240 out.set_subvector(1, reverse(out(1, N - 1))); 00241 00242 // Normalise output data 00243 out *= factor; 00244 } 00245 00246 #endif // defined(HAVE_FFT_ACML) 00247 00248 00249 #if defined(HAVE_FFTW3) 00250 00251 //--------------------------------------------------------------------------- 00252 // FFT/IFFT based on FFTW 00253 //--------------------------------------------------------------------------- 00254 00255 void fft(const cvec &in, cvec &out) 00256 { 00257 static int N = 0; 00258 static fftw_plan p = NULL; 00259 out.set_size(in.size(), false); 00260 00261 if (N != in.size()) { 00262 N = in.size(); 00263 if (p != NULL) 00264 fftw_destroy_plan(p); // destroy the previous plan 00265 // create a new plan 00266 p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(), 00267 (fftw_complex *)out._data(), 00268 FFTW_FORWARD, FFTW_ESTIMATE); 00269 } 00270 00271 // compute FFT using the GURU FFTW interface 00272 fftw_execute_dft(p, (fftw_complex *)in._data(), 00273 (fftw_complex *)out._data()); 00274 } 00275 00276 void ifft(const cvec &in, cvec &out) 00277 { 00278 static int N = 0; 00279 static double inv_N; 00280 static fftw_plan p = NULL; 00281 out.set_size(in.size(), false); 00282 00283 if (N != in.size()) { 00284 N = in.size(); 00285 inv_N = 1.0 / N; 00286 if (p != NULL) 00287 fftw_destroy_plan(p); // destroy the previous plan 00288 // create a new plan 00289 p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(), 00290 (fftw_complex *)out._data(), 00291 FFTW_BACKWARD, FFTW_ESTIMATE); 00292 } 00293 00294 // compute IFFT using the GURU FFTW interface 00295 fftw_execute_dft(p, (fftw_complex *)in._data(), 00296 (fftw_complex *)out._data()); 00297 00298 // scale output 00299 out *= inv_N; 00300 } 00301 00302 void fft_real(const vec &in, cvec &out) 00303 { 00304 static int N = 0; 00305 static fftw_plan p = NULL; 00306 out.set_size(in.size(), false); 00307 00308 if (N != in.size()) { 00309 N = in.size(); 00310 if (p != NULL) 00311 fftw_destroy_plan(p); //destroy the previous plan 00312 00313 // create a new plan 00314 p = fftw_plan_dft_r2c_1d(N, (double *)in._data(), 00315 (fftw_complex *)out._data(), 00316 FFTW_ESTIMATE); 00317 } 00318 00319 // compute FFT using the GURU FFTW interface 00320 fftw_execute_dft_r2c(p, (double *)in._data(), 00321 (fftw_complex *)out._data()); 00322 00323 // Real FFT does not compute the 2nd half of the FFT points because it 00324 // is redundant to the 1st half. However, we want all of the data so we 00325 // fill it in. This is consistent with Matlab's functionality 00326 int offset = ceil_i(in.size() / 2.0); 00327 int n_elem = in.size() - offset; 00328 for (int i = 0; i < n_elem; ++i) { 00329 out(offset + i) = std::conj(out(n_elem - i)); 00330 } 00331 } 00332 00333 void ifft_real(const cvec &in, vec & out) 00334 { 00335 static int N = 0; 00336 static double inv_N; 00337 static fftw_plan p = NULL; 00338 out.set_size(in.size(), false); 00339 00340 if (N != in.size()) { 00341 N = in.size(); 00342 inv_N = 1.0 / N; 00343 if (p != NULL) 00344 fftw_destroy_plan(p); // destroy the previous plan 00345 00346 // create a new plan 00347 p = fftw_plan_dft_c2r_1d(N, (fftw_complex *)in._data(), 00348 (double *)out._data(), 00349 FFTW_ESTIMATE | FFTW_PRESERVE_INPUT); 00350 } 00351 00352 // compute IFFT using the GURU FFTW interface 00353 fftw_execute_dft_c2r(p, (fftw_complex *)in._data(), 00354 (double *)out._data()); 00355 00356 out *= inv_N; 00357 } 00358 00359 //--------------------------------------------------------------------------- 00360 // DCT/IDCT based on FFTW 00361 //--------------------------------------------------------------------------- 00362 00363 void dct(const vec &in, vec &out) 00364 { 00365 static int N; 00366 static fftw_plan p = NULL; 00367 out.set_size(in.size(), false); 00368 00369 if (N != in.size()) { 00370 N = in.size(); 00371 if (p != NULL) 00372 fftw_destroy_plan(p); // destroy the previous plan 00373 00374 // create a new plan 00375 p = fftw_plan_r2r_1d(N, (double *)in._data(), 00376 (double *)out._data(), 00377 FFTW_REDFT10, FFTW_ESTIMATE); 00378 } 00379 00380 // compute FFT using the GURU FFTW interface 00381 fftw_execute_r2r(p, (double *)in._data(), (double *)out._data()); 00382 00383 // Scale to matlab definition format 00384 out /= std::sqrt(2.0 * N); 00385 out(0) /= std::sqrt(2.0); 00386 } 00387 00388 // IDCT 00389 void idct(const vec &in, vec &out) 00390 { 00391 static int N; 00392 static fftw_plan p = NULL; 00393 out = in; 00394 00395 // Rescale to FFTW format 00396 out(0) *= std::sqrt(2.0); 00397 out /= std::sqrt(2.0 * in.size()); 00398 00399 if (N != in.size()) { 00400 N = in.size(); 00401 if (p != NULL) 00402 fftw_destroy_plan(p); // destroy the previous plan 00403 00404 // create a new plan 00405 p = fftw_plan_r2r_1d(N, (double *)out._data(), 00406 (double *)out._data(), 00407 FFTW_REDFT01, FFTW_ESTIMATE); 00408 } 00409 00410 // compute FFT using the GURU FFTW interface 00411 fftw_execute_r2r(p, (double *)out._data(), (double *)out._data()); 00412 } 00413 00414 #endif // defined(HAVE_FFTW3) 00415 00416 00417 #if defined(HAVE_FFT_MKL) || defined(HAVE_FFT_ACML) 00418 00419 //--------------------------------------------------------------------------- 00420 // DCT/IDCT based on MKL or ACML 00421 //--------------------------------------------------------------------------- 00422 00423 void dct(const vec &in, vec &out) 00424 { 00425 int N = in.size(); 00426 if (N == 1) 00427 out = in; 00428 else { 00429 cvec c = fft_real(concat(in, reverse(in))); 00430 c.set_size(N, true); 00431 for (int i = 0; i < N; i++) { 00432 c(i) *= std::complex<double>(std::cos(pi * i / N / 2), std::sin(-pi * i / N / 2)) 00433 / std::sqrt(2.0 * N); 00434 } 00435 out = real(c); 00436 out(0) /= std::sqrt(2.0); 00437 } 00438 } 00439 00440 void idct(const vec &in, vec &out) 00441 { 00442 int N = in.size(); 00443 if (N == 1) 00444 out = in; 00445 else { 00446 cvec c = to_cvec(in); 00447 c.set_size(2*N, true); 00448 c(0) *= std::sqrt(2.0); 00449 for (int i = 0; i < N; i++) { 00450 c(i) *= std::complex<double>(std::cos(pi * i / N / 2), std::sin(pi * i / N / 2)) 00451 * std::sqrt(2.0 * N); 00452 } 00453 for (int i = N - 1; i >= 1; i--) { 00454 c(c.size() - i) = c(i) * std::complex<double>(std::cos(pi * i / N), 00455 std::sin(-pi * i / N)); 00456 } 00457 out = ifft_real(c); 00458 out.set_size(N, true); 00459 } 00460 } 00461 00462 #endif // defined(HAVE_FFT_MKL) || defined(HAVE_FFT_ACML) 00463 00464 00465 #if !defined(HAVE_FFT) 00466 00467 void fft(const cvec &in, cvec &out) 00468 { 00469 it_error("FFT library is needed to use fft() function"); 00470 } 00471 00472 void ifft(const cvec &in, cvec &out) 00473 { 00474 it_error("FFT library is needed to use ifft() function"); 00475 } 00476 00477 void fft_real(const vec &in, cvec &out) 00478 { 00479 it_error("FFT library is needed to use fft_real() function"); 00480 } 00481 00482 void ifft_real(const cvec &in, vec & out) 00483 { 00484 it_error("FFT library is needed to use ifft_real() function"); 00485 } 00486 00487 void dct(const vec &in, vec &out) 00488 { 00489 it_error("FFT library is needed to use dct() function"); 00490 } 00491 00492 void idct(const vec &in, vec &out) 00493 { 00494 it_error("FFT library is needed to use idct() function"); 00495 } 00496 00497 #endif // !defined(HAVE_FFT) 00498 00499 cvec fft(const cvec &in) 00500 { 00501 cvec out; 00502 fft(in, out); 00503 return out; 00504 } 00505 00506 cvec fft(const cvec &in, const int N) 00507 { 00508 cvec in2 = in; 00509 cvec out; 00510 in2.set_size(N, true); 00511 fft(in2, out); 00512 return out; 00513 } 00514 00515 cvec ifft(const cvec &in) 00516 { 00517 cvec out; 00518 ifft(in, out); 00519 return out; 00520 } 00521 00522 cvec ifft(const cvec &in, const int N) 00523 { 00524 cvec in2 = in; 00525 cvec out; 00526 in2.set_size(N, true); 00527 ifft(in2, out); 00528 return out; 00529 } 00530 00531 cvec fft_real(const vec& in) 00532 { 00533 cvec out; 00534 fft_real(in, out); 00535 return out; 00536 } 00537 00538 cvec fft_real(const vec& in, const int N) 00539 { 00540 vec in2 = in; 00541 cvec out; 00542 in2.set_size(N, true); 00543 fft_real(in2, out); 00544 return out; 00545 } 00546 00547 vec ifft_real(const cvec &in) 00548 { 00549 vec out; 00550 ifft_real(in, out); 00551 return out; 00552 } 00553 00554 vec ifft_real(const cvec &in, const int N) 00555 { 00556 cvec in2 = in; 00557 in2.set_size(N, true); 00558 vec out; 00559 ifft_real(in2, out); 00560 return out; 00561 } 00562 00563 vec dct(const vec &in) 00564 { 00565 vec out; 00566 dct(in, out); 00567 return out; 00568 } 00569 00570 vec idct(const vec &in) 00571 { 00572 vec out; 00573 idct(in, out); 00574 return out; 00575 } 00576 00577 00578 // ---------------------------------------------------------------------- 00579 // Instantiation 00580 // ---------------------------------------------------------------------- 00581 00582 template vec dht(const vec &v); 00583 template cvec dht(const cvec &v); 00584 00585 template void dht(const vec &vin, vec &vout); 00586 template void dht(const cvec &vin, cvec &vout); 00587 00588 template void self_dht(vec &v); 00589 template void self_dht(cvec &v); 00590 00591 template vec dwht(const vec &v); 00592 template cvec dwht(const cvec &v); 00593 00594 template void dwht(const vec &vin, vec &vout); 00595 template void dwht(const cvec &vin, cvec &vout); 00596 00597 template void self_dwht(vec &v); 00598 template void self_dwht(cvec &v); 00599 00600 template mat dht2(const mat &m); 00601 template cmat dht2(const cmat &m); 00602 00603 template mat dwht2(const mat &m); 00604 template cmat dwht2(const cmat &m); 00605 00606 } // namespace itpp 00607
Generated on Sat Jul 9 2011 15:21:33 for IT++ by Doxygen 1.7.4