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

Generated on Wed Jul 27 2011 16:27:05 for IT++ by Doxygen 1.7.4