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
Generated on Wed Apr 18 11:45:34 2007 for IT++ by Doxygen 1.5.2