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