IT++ Logo Newcom Logo

transforms.cpp

Go to the documentation of this file.
00001 
00037 #ifndef _MSC_VER
00038 #  include <itpp/config.h>
00039 #else
00040 #  include <itpp/config_msvc.h>
00041 #endif
00042 
00043 #if defined(HAVE_FFT_MKL8)
00044 #  include <mkl_dfti.h>
00045 #elif defined(HAVE_FFT_ACML)
00046 #  include <acml.h>
00047 #elif defined(HAVE_FFTW3)
00048 #  include <fftw3.h>
00049 #endif
00050 
00051 #include <itpp/base/matfunc.h>
00052 #include <itpp/base/transforms.h>
00053 #include <itpp/base/elmatfunc.h>
00054 
00055 
00056 namespace itpp { 
00057 
00058 #ifdef HAVE_FFT_MKL8
00059 
00060   //---------------------------------------------------------------------------
00061   // FFT/IFFT based on MKL
00062   //---------------------------------------------------------------------------
00063 
00064   void fft(const cvec &in, cvec &out) 
00065   {
00066     static 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) DftiFreeDescriptor(&fft_handle);
00073       DftiCreateDescriptor(&fft_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, N);
00074       DftiSetValue(fft_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00075       DftiCommitDescriptor(fft_handle);
00076     }
00077     DftiComputeForward(fft_handle, (void *)in._data(), out._data());
00078   }
00079 
00080   void ifft(const cvec &in, cvec &out)
00081   {
00082     static 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) DftiFreeDescriptor(&fft_handle);
00089       DftiCreateDescriptor(&fft_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, N);
00090       DftiSetValue(fft_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00091       DftiSetValue(fft_handle, DFTI_BACKWARD_SCALE, 1.0/N);
00092       DftiCommitDescriptor(fft_handle);
00093     }
00094     DftiComputeBackward(fft_handle, (void *)in._data(), out._data());
00095   }
00096 
00097   void fft_real(const vec &in, cvec &out)
00098   {
00099     static 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) DftiFreeDescriptor(&fft_handle);
00106       DftiCreateDescriptor(&fft_handle, DFTI_DOUBLE, DFTI_REAL, 1, N);
00107       DftiSetValue(fft_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00108       DftiCommitDescriptor(fft_handle);
00109     }
00110     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 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) DftiFreeDescriptor(&fft_handle);
00130       DftiCreateDescriptor( &fft_handle, DFTI_DOUBLE, DFTI_REAL, 1, N);
00131       DftiSetValue(fft_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
00132       DftiSetValue(fft_handle, DFTI_BACKWARD_SCALE, 1.0/N);
00133       DftiCommitDescriptor(fft_handle);
00134     }
00135     DftiComputeBackward(fft_handle, (void *)in._data(), out._data());
00136   }
00137 
00138   //---------------------------------------------------------------------------
00139   // DCT/IDCT based on MKL
00140   //---------------------------------------------------------------------------
00141 
00142   void dct(const vec &in, vec &out)
00143   {
00144     int N = in.size();
00145     if (N == 1)
00146       out = in;
00147     else {
00148       cvec c = fft_real(concat(in, reverse(in)));
00149       c.set_size(N, true);
00150       for (int i = 0; i < N; i++) {
00151         c(i) *= std::complex<double>(std::cos(pi*i/N/2), std::sin(-pi*i/N/2))
00152           / std::sqrt(2.0 * N);
00153       }
00154       out = real(c);
00155       out(0) /= std::sqrt(2.0);
00156     }
00157   }
00158 
00159   void idct(const vec &in, vec &out)
00160   {
00161     int N = in.size();
00162     if (N == 1)
00163       out = in;
00164     else {
00165       cvec c = to_cvec(in);
00166       c.set_size(2*N, true);
00167       c(0) *= std::sqrt(2.0);
00168       for (int i = 0; i < N; i++) {
00169         c(i) *= std::complex<double>(std::cos(pi*i/N/2), std::sin(pi*i/N/2))
00170           * std::sqrt(2.0 * N);
00171       }
00172       for (int i = N - 1; i >= 1; i--) {
00173         c(c.size() - i) = c(i) * std::complex<double>(std::cos(pi*i/N), 
00174                                                       std::sin(-pi*i/N));
00175       }
00176       out = ifft_real(c);
00177       out.set_size(N, true);
00178     }
00179   }
00180 
00181 #elif defined(HAVE_FFT_ACML)
00182 
00183   //---------------------------------------------------------------------------
00184   // FFT/IFFT based on ACML
00185   //---------------------------------------------------------------------------
00186 
00187   void fft(const cvec &in, cvec &out) 
00188   {
00189     static int N = 0;
00190     static cvec *comm_ptr = NULL;
00191     int info;
00192     out.set_size(in.size(), false);
00193 
00194     if (N != in.size()) {
00195       N = in.size();
00196       if (comm_ptr != NULL)
00197         delete comm_ptr;
00198       comm_ptr = new cvec(5 * in.size() + 100);
00199       zfft1dx(0, 1.0, false, N, (doublecomplex *)in._data(), 1, 
00200               (doublecomplex *)out._data(), 1, 
00201               (doublecomplex *)comm_ptr->_data(), &info);
00202     }
00203     zfft1dx(-1, 1.0, false, N, (doublecomplex *)in._data(), 1, 
00204             (doublecomplex *)out._data(), 1, 
00205             (doublecomplex *)comm_ptr->_data(), &info);
00206   }
00207 
00208   void ifft(const cvec &in, cvec &out) 
00209   {
00210     static int N = 0;
00211     static cvec *comm_ptr = NULL;
00212     int info;
00213     out.set_size(in.size(), false);
00214 
00215     if (N != in.size()) {
00216       N = in.size();
00217       if (comm_ptr != NULL)
00218         delete comm_ptr;
00219       comm_ptr = new cvec(5 * in.size() + 100);
00220       zfft1dx(0, 1.0/N, false, N, (doublecomplex *)in._data(), 1, 
00221               (doublecomplex *)out._data(), 1, 
00222               (doublecomplex *)comm_ptr->_data(), &info);
00223     }
00224     zfft1dx(1, 1.0/N, false, N, (doublecomplex *)in._data(), 1, 
00225             (doublecomplex *)out._data(), 1, 
00226             (doublecomplex *)comm_ptr->_data(), &info);
00227   }
00228 
00229   void fft_real(const vec &in, cvec &out)
00230   {
00231     static int N = 0;
00232     static double factor = 0;
00233     static vec *comm_ptr = NULL;
00234     int info;
00235     vec out_re = in;
00236 
00237     if (N != in.size()) {
00238       N = in.size();
00239       factor = std::sqrt(static_cast<double>(N));
00240       if (comm_ptr != NULL)
00241         delete comm_ptr;
00242       comm_ptr = new vec(5 * in.size() + 100);
00243       dzfft(0, N, out_re._data(), comm_ptr->_data(), &info);
00244     }
00245     dzfft(1, N, out_re._data(), comm_ptr->_data(), &info);
00246 
00247     // Normalise output data
00248     out_re *= factor;
00249 
00250     // Convert the real Hermitian DZFFT's output to the Matlab's complex form
00251     vec out_im(in.size()); out_im.zeros();
00252     out.set_size(in.size(), false);
00253     out_im.set_subvector(1, reverse(out_re(N/2 + 1, N-1)));
00254     out_im.set_subvector(N/2 + 1, -out_re(N/2 + 1, N-1));
00255     out_re.set_subvector(N/2 + 1, reverse(out_re(1, (N-1)/2)));
00256     out = to_cvec(out_re, out_im);
00257   }
00258 
00259   void ifft_real(const cvec &in, vec &out)
00260   {
00261     static int N = 0;
00262     static double factor = 0;
00263     static vec *comm_ptr = NULL;
00264     int info;
00265 
00266     // Convert Matlab's complex input to the real Hermitian form
00267     out.set_size(in.size());
00268     out.set_subvector(0, real(in(0, in.size()/2)));
00269     out.set_subvector(in.size()/2 + 1, -imag(in(in.size()/2 + 1, in.size()-1)));
00270 
00271     if (N != in.size()) {
00272       N = in.size();
00273       factor = 1.0 / std::sqrt(static_cast<double>(N));
00274       if (comm_ptr != NULL)
00275         delete comm_ptr;
00276       comm_ptr = new vec(5 * in.size() + 100);
00277       zdfft(0, N, out._data(), comm_ptr->_data(), &info);
00278     }
00279     zdfft(1, N, out._data(), comm_ptr->_data(), &info);
00280     out.set_subvector(1, reverse(out(1, N-1)));
00281 
00282     // Normalise output data
00283     out *= factor;
00284   }
00285 
00286   //---------------------------------------------------------------------------
00287   // DCT/IDCT based on ACML
00288   //---------------------------------------------------------------------------
00289 
00290   void dct(const vec &in, vec &out)
00291   {
00292     int N = in.size();
00293     if (N == 1)
00294       out = in;
00295     else {
00296       cvec c = fft_real(concat(in, reverse(in)));
00297       c.set_size(N, true);
00298       for (int i = 0; i < N; i++) {
00299         c(i) *= std::complex<double>(std::cos(pi*i/N/2), std::sin(-pi*i/N/2))
00300           / std::sqrt(2.0 * N);
00301       }
00302       out = real(c);
00303       out(0) /= std::sqrt(2.0);
00304     }
00305   }
00306 
00307   void idct(const vec &in, vec &out)
00308   {
00309     int N = in.size();
00310     if (N == 1)
00311       out = in;
00312     else {
00313       cvec c = to_cvec(in);
00314       c.set_size(2*N, true);
00315       c(0) *= std::sqrt(2.0);
00316       for (int i = 0; i < N; i++) {
00317         c(i) *= std::complex<double>(std::cos(pi*i/N/2), std::sin(pi*i/N/2))
00318           * std::sqrt(2.0 * N);
00319       }
00320       for (int i = N - 1; i >= 1; i--) {
00321         c(c.size() - i) = c(i) * std::complex<double>(std::cos(pi*i/N), 
00322                                                       std::sin(-pi*i/N));
00323       }
00324       out = ifft_real(c);
00325       out.set_size(N, true);
00326     }
00327   }
00328 
00329 #elif defined(HAVE_FFTW3)
00330 
00331   //---------------------------------------------------------------------------
00332   // FFT/IFFT based on FFTW
00333   //---------------------------------------------------------------------------
00334 
00335   /* Note: The FFTW_UNALIGNED flag is set when calling the FFTW
00336      library, as memory alignment may change between calls (this was
00337      done to fix bug 1418707).
00338 
00339      Setting the FFTW_UNALIGNED flag fixes bug 1418707 for the time
00340      being but it does not result in maximum possible performance in
00341      all cases. It appears that a solution that achieves maximum
00342      performance from FFTW in IT++ would require either to (i)
00343      redefine IT++ memory management functions to use the type of
00344      memory alignment required by FFTW, (ii) implement a memory
00345      re-allocation and data copying mechanism in the IT++/FFTW
00346      interface function to ensure FFTW is called with properly aligned
00347      data, or (iii) force the IT++/FFTW interface function to recreate
00348      the FFT plan whenever memory alignment has changed.  None of
00349      these solutions was found attractive.
00350      
00351   */
00352 
00353   void fft(const cvec &in, cvec &out)
00354   {
00355     static int N = 0;
00356     static fftw_plan p = NULL;
00357     out.set_size(in.size(), false);
00358 
00359     if (N != in.size()) {
00360       N = in.size();
00361       if (p != NULL)
00362         fftw_destroy_plan(p); // destroy the previous plan
00363       // create a new plan
00364       p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(), (fftw_complex *)out._data(),
00365                            FFTW_FORWARD, FFTW_ESTIMATE | FFTW_UNALIGNED);
00366     }
00367 
00368     // compute FFT using the GURU FFTW interface
00369     fftw_execute_dft(p, (fftw_complex *)in._data(), (fftw_complex *)out._data());
00370   }
00371 
00372   void ifft(const cvec &in, cvec &out)
00373   {
00374     static int N = 0;
00375     static double inv_N;
00376     static fftw_plan p = NULL;
00377     out.set_size(in.size(), false);
00378 
00379     if (N != in.size()) {
00380       N = in.size();
00381       inv_N = 1.0/N;
00382       if (p != NULL)
00383         fftw_destroy_plan(p); // destroy the previous plan
00384       // create a new plan
00385       p = fftw_plan_dft_1d(N, (fftw_complex *)in._data(), (fftw_complex *)out._data(),
00386                            FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_UNALIGNED);
00387     }
00388 
00389     // compute IFFT using the GURU FFTW interface
00390     fftw_execute_dft(p, (fftw_complex *)in._data(), (fftw_complex *)out._data());
00391 
00392     // scale output
00393     out *= inv_N;
00394   }
00395 
00396   void fft_real(const vec &in, cvec &out)
00397   {
00398     static int N = 0;
00399     static fftw_plan p = NULL;
00400     out.set_size(in.size(), false);
00401 
00402     if (N != in.size()) {
00403       N = in.size();
00404       if (p!= NULL)
00405         fftw_destroy_plan(p); //destroy the previous plan
00406 
00407       // create a new plan
00408       p = fftw_plan_dft_r2c_1d(N, (double *)in._data(), (fftw_complex *)out._data(),
00409                                FFTW_ESTIMATE | FFTW_UNALIGNED);
00410     }
00411 
00412     // compute FFT using the GURU FFTW interface
00413     fftw_execute_dft_r2c(p, (double *)in._data(), (fftw_complex *)out._data());
00414 
00415     // Real FFT does not compute the 2nd half of the FFT points because it
00416     // is redundant to the 1st half. However, we want all of the data so we
00417     // fill it in. This is consistent with Matlab's functionality
00418     int istart = ceil_i(in.size() / 2.0);
00419     int iend = in.size() - 1;
00420     int idelta = iend - istart + 1;
00421     out.set_subvector(istart, iend, reverse(conj(out(1, idelta))));
00422   }
00423 
00424   void ifft_real(const cvec &in, vec & out)
00425   {
00426     static int N = 0;
00427     static double inv_N;
00428     static fftw_plan p = NULL;
00429     out.set_size(in.size(), false);
00430 
00431     if (N != in.size()) {
00432       N = in.size();
00433       inv_N = 1.0/N;
00434       if (p != NULL) 
00435         fftw_destroy_plan(p); // destroy the previous plan
00436 
00437       // create a new plan
00438       p = fftw_plan_dft_c2r_1d(N, (fftw_complex *)in._data(), (double *)out._data(), 
00439                                FFTW_ESTIMATE | FFTW_UNALIGNED);
00440     }
00441 
00442     // compute IFFT using the GURU FFTW interface
00443     fftw_execute_dft_c2r(p, (fftw_complex *)in._data(), (double *)out._data());
00444 
00445     out *= inv_N;
00446   }
00447 
00448   //---------------------------------------------------------------------------
00449   // DCT/IDCT based on FFTW
00450   //---------------------------------------------------------------------------
00451 
00452   void dct(const vec &in, vec &out)
00453   {
00454     static int N;
00455     static fftw_plan p = NULL;
00456     out.set_size(in.size(), false);
00457 
00458     if (N != in.size()) {
00459       N = in.size();
00460       if (p!= NULL)
00461         fftw_destroy_plan(p); // destroy the previous plan
00462 
00463       // create a new plan
00464       p = fftw_plan_r2r_1d(N, (double *)in._data(), (double *)out._data(), 
00465                            FFTW_REDFT10, FFTW_ESTIMATE | FFTW_UNALIGNED);
00466     }
00467 
00468     // compute FFT using the GURU FFTW interface
00469     fftw_execute_r2r(p, (double *)in._data(), (double *)out._data());
00470 
00471     // Scale to matlab definition format
00472     out /= std::sqrt(2.0 * N);
00473     out(0) /= std::sqrt(2.0);
00474   }
00475 
00476   // IDCT
00477   void idct(const vec &in, vec &out)
00478   {
00479     static int N;
00480     static fftw_plan p = NULL;
00481     out = in;
00482 
00483     // Rescale to FFTW format
00484     out(0) *= std::sqrt(2.0);
00485     out /= std::sqrt(2.0 * in.size());
00486 
00487     if (N != in.size()) {
00488       N = in.size();
00489       if (p != NULL)
00490         fftw_destroy_plan(p); // destroy the previous plan
00491       
00492       // create a new plan
00493       p = fftw_plan_r2r_1d(N, (double *)out._data(), (double *)out._data(),
00494                            FFTW_REDFT01, FFTW_ESTIMATE | FFTW_UNALIGNED);
00495     }
00496 
00497     // compute FFT using the GURU FFTW interface
00498     fftw_execute_r2r(p, (double *)out._data(), (double *)out._data());
00499   }
00500 
00501 #else
00502 
00503   void fft(const cvec &in, cvec &out)
00504   {
00505     it_error("FFT library is needed to use fft() function");
00506   }
00507 
00508   void ifft(const cvec &in, cvec &out)
00509   {
00510     it_error("FFT library is needed to use ifft() function");
00511   }
00512 
00513   void fft_real(const vec &in, cvec &out)
00514   {
00515     it_error("FFT library is needed to use fft_real() function");
00516   }
00517 
00518   void ifft_real(const cvec &in, vec & out)
00519   {
00520     it_error("FFT library is needed to use ifft_real() function");
00521   }
00522 
00523   void dct(const vec &in, vec &out)
00524   {
00525     it_error("FFT library is needed to use dct() function");
00526   }
00527 
00528   void idct(const vec &in, vec &out)
00529   {
00530     it_error("FFT library is needed to use idct() function");
00531   }
00532 
00533 #endif // HAVE_FFT_MKL8 or HAVE_FFT_ACML or HAVE_FFTW3
00534 
00535   cvec fft(const cvec &in) 
00536   { 
00537     cvec out; 
00538     fft(in, out); 
00539     return out; 
00540   }
00541 
00542   cvec fft(const cvec &in, const int N)  
00543   { 
00544     cvec in2 = in;
00545     cvec out;
00546     in2.set_size(N, true); 
00547     fft(in2, out); 
00548     return out; 
00549   }
00550 
00551   cvec ifft(const cvec &in)  
00552   { 
00553     cvec out;
00554     ifft(in, out);
00555     return out; 
00556   }
00557 
00558   cvec ifft(const cvec &in, const int N)  
00559   { 
00560     cvec in2 = in;
00561     cvec out;
00562     in2.set_size(N, true); 
00563     ifft(in2, out); 
00564     return out; 
00565   }
00566 
00567   cvec fft_real(const vec& in)  
00568   { 
00569     cvec out;
00570     fft_real(in, out);
00571     return out; 
00572   }
00573 
00574   cvec fft_real(const vec& in, const int N)  
00575   { 
00576     vec in2 = in;
00577     cvec out;
00578     in2.set_size(N, true); 
00579     fft_real(in2, out);
00580     return out; 
00581   }
00582 
00583   vec ifft_real(const cvec &in) 
00584   {
00585     vec out;
00586     ifft_real(in, out);
00587     return out; 
00588   }
00589 
00590   vec ifft_real(const cvec &in, const int N) 
00591   {
00592     cvec in2 = in; 
00593     in2.set_size(N, true);
00594     vec out;
00595     ifft_real(in2, out);
00596     return out; 
00597   }
00598 
00599   vec dct(const vec &in)
00600   { 
00601     vec out;
00602     dct(in, out);
00603     return out; 
00604   }
00605 
00606   vec idct(const vec &in) 
00607   { 
00608     vec out;
00609     idct(in, out);
00610     return out; 
00611   }
00612 
00613 
00614   template <class T>
00615   Vec<T> dht(const Vec<T> &v)
00616   {
00617     Vec<T> ret(v.size());
00618     dht(v, ret);
00619     return ret;
00620   }
00621 
00622   template <class T>
00623   void bitrv(Vec<T> &out) {
00624     int i,j,N1,K,N=out.size();
00625     T TEMP;
00626 
00627     j=0;
00628     N1=N-1;
00629     for (i=0;i<N1;i++) {
00630       if (i<j) {
00631         TEMP=out[j];
00632         out[j]=out[i];
00633         out[i]=TEMP;
00634       }
00635       K=N/2;
00636       while (K<=j) {
00637         j=j-K;
00638         K=K/2;
00639       }
00640       j=j+K;
00641     }
00642   }
00643 
00644   template <class T>
00645   void dht(const Vec<T> &vin, Vec<T> &vout)
00646   {
00647     T t;
00648     int m,N,l,k,j,ib,i;
00649 
00650     N = vin.size();
00651     vout.set_size(N);
00652 
00653     m = levels2bits(N);
00654     it_assert1((1<<m)==N, "dht: The vector size must be a power of two!");
00655 
00656     // This step is separated because it copies vin to vout
00657     for (ib=0; ib<N; ib+=2) {
00658       vout(ib) = vin(ib) + vin(ib+1);
00659       vout(ib+1) = vin(ib) - vin(ib+1);
00660     }
00661     N /= 2;
00662 
00663     l = 2;    
00664     for (i=1; i<m; i++) {
00665       N /= 2;
00666       ib = 0;
00667       for (k=0; k<N; k++) {
00668         for (j=0; j<l; j++) {
00669           t = vout(ib+j);
00670           vout(ib+j) += vout(ib+j+l);
00671           vout(ib+j+l) = t - vout(ib+j+l);
00672         }
00673         ib += 2*l;
00674       }
00675       l *= 2;
00676     }
00677 
00678     vout /= T(std::sqrt(double(vin.size())));
00679   }
00680 
00681   template <class T>
00682   void self_dht(Vec<T> &v)
00683   {
00684     T t;
00685     int m,N,l=1,k,j,ib,i;
00686 
00687     N = v.size();
00688 
00689     m = levels2bits(N);
00690     it_assert1((1<<m)==N, "dht: The vector size must be a power of two!");
00691 
00692     for (i=0; i<m; i++) {
00693       N /= 2;
00694       ib = 0;
00695       for (k=0; k<N; k++) {
00696         for (j=0; j<l; j++) {
00697           t = v(ib+j);
00698           v(ib+j) += v(ib+j+l);
00699           v(ib+j+l) = t - v(ib+j+l);
00700         }
00701         ib += 2*l;
00702       }
00703       l *= 2;
00704     }
00705 
00706     v /= T(std::sqrt(double(v.size())));
00707   }
00708 
00709   template <class T>
00710   Vec<T> dwht(const Vec<T> &v)
00711   {
00712     Vec<T> ret(v.size());
00713 
00714     dwht(v, ret);
00715 
00716     return ret;
00717   }
00718 
00719   template <class T>
00720   void dwht(const Vec<T> &vin, Vec<T> &vout)
00721   {
00722     dht(vin, vout);
00723     bitrv(vout);
00724   }
00725 
00726 
00727   template <class T>
00728   void self_dwht(Vec<T> &v)
00729   {
00730     self_dht(v);
00731     bitrv(v);
00732   }
00733 
00734   template <class T>
00735   Mat<T> dht2(const Mat<T> &m)
00736   {
00737     Mat<T> ret(m.rows(), m.cols());
00738     Vec<T> v;
00739     int i;
00740 
00741     for (i=0; i<m.rows(); i++) {
00742       v = m.get_row(i);
00743       self_dht(v);
00744       ret.set_row(i, v);
00745     }
00746     for (i=0; i<m.cols(); i++) {
00747       v = ret.get_col(i);
00748       self_dht(v);
00749       ret.set_col(i, v);
00750     }
00751 
00752     return transpose(ret);
00753   }
00754 
00755   template <class T>
00756   Mat<T> dwht2(const Mat<T> &m)
00757   {
00758     Mat<T> ret(m.rows(), m.cols());
00759     Vec<T> v;
00760     int i;
00761 
00762     for (i=0; i<m.rows(); i++) {
00763       v = m.get_row(i);
00764       self_dwht(v);
00765       ret.set_row(i, v);
00766     }
00767     for (i=0; i<m.cols(); i++) {
00768       v = ret.get_col(i);
00769       self_dwht(v);
00770       ret.set_col(i, v);
00771     }
00772 
00773     return transpose(ret);
00774   }
00775 
00777   //
00778   //  template instantiation
00779   //
00781 
00782   template vec dht(const vec &v);
00783   template cvec dht(const cvec &v);
00784 
00785   template void dht(const vec &vin, vec &vout);
00786   template void dht(const cvec &vin, cvec &vout);
00787 
00788   template void self_dht(vec &v);
00789   template void self_dht(cvec &v);
00790 
00791   template vec dwht(const vec &v);
00792   template cvec dwht(const cvec &v);
00793 
00794   template void dwht(const vec &vin, vec &vout);
00795   template void dwht(const cvec &vin, cvec &vout);
00796 
00797   template void self_dwht(vec &v);
00798   template void self_dwht(cvec &v);
00799 
00800   template mat  dht2(const mat &m);
00801   template cmat dht2(const cmat &m);
00802 
00803   template mat  dwht2(const mat &m);
00804   template cmat dwht2(const cmat &m);
00805 
00806 } // namespace itpp
SourceForge Logo

Generated on Wed Apr 18 11:19:58 2007 for IT++ by Doxygen 1.5.2