IT++ Logo

transforms.cpp

Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Tue Nov 23 08:45:11 2010 for IT++ by Doxygen 1.6.1