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