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