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