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 #  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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jan 20 23:03:06 2010 for IT++ by Doxygen 1.6.2