IT++ Logo Newcom Logo

transforms.cpp

Go to the documentation of this file.
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_MKL8)
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_MKL8
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_MKL8 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
SourceForge Logo

Generated on Fri Jun 8 00:27:14 2007 for IT++ by Doxygen 1.5.2