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 #  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 
SourceForge Logo

Generated on Sat Apr 19 10:43:56 2008 for IT++ by Doxygen 1.5.5