00001 00030 #ifndef FILTER_H 00031 #define FILTER_H 00032 00033 #ifndef _MSC_VER 00034 # include <itpp/config.h> 00035 #else 00036 # include <itpp/config_msvc.h> 00037 #endif 00038 00039 #include <itpp/base/vec.h> 00040 00041 00042 namespace itpp { 00043 00059 template <class T1, class T2, class T3> 00060 class Filter { 00061 public: 00063 Filter() {} 00065 virtual T3 operator()(const T1 Sample) { return filter(Sample); } 00067 virtual Vec<T3> operator()(const Vec<T1> &v); 00069 virtual ~Filter() {} 00070 protected: 00075 virtual T3 filter(const T1 Sample)=0; 00076 }; 00077 00102 template <class T1, class T2, class T3> 00103 class MA_Filter : public Filter<T1,T2,T3> { 00104 public: 00106 explicit MA_Filter(); 00108 explicit MA_Filter(const Vec<T2> &b); 00110 virtual ~MA_Filter() { } 00112 Vec<T2> get_coeffs() const { return coeffs; } 00114 void set_coeffs(const Vec<T2> &b); 00116 void clear() { mem.clear(); } 00118 Vec<T3> get_state() const; 00120 void set_state(const Vec<T3> &state); 00121 00122 private: 00123 virtual T3 filter(const T1 Sample); 00124 00125 Vec<T3> mem; 00126 Vec<T2> coeffs; 00127 int inptr; 00128 bool init; 00129 }; 00130 00155 template <class T1, class T2, class T3> 00156 class AR_Filter : public Filter<T1,T2,T3> { 00157 public: 00159 explicit AR_Filter(); 00161 explicit AR_Filter(const Vec<T2> &a); 00163 virtual ~AR_Filter() { } 00165 Vec<T2> get_coeffs() const { return coeffs; } 00167 void set_coeffs(const Vec<T2> &a); 00169 void clear() { mem.clear(); } 00171 Vec<T3> get_state() const; 00173 void set_state(const Vec<T3> &state); 00174 00175 private: 00176 virtual T3 filter(const T1 Sample); 00177 00178 Vec<T3> mem; 00179 Vec<T2> coeffs; 00180 T2 a0; 00181 int inptr; 00182 bool init; 00183 }; 00184 00185 00212 template <class T1, class T2, class T3> 00213 class ARMA_Filter : public Filter<T1,T2,T3> { 00214 public: 00216 explicit ARMA_Filter(); 00218 explicit ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a); 00220 virtual ~ARMA_Filter() { } 00222 Vec<T2> get_coeffs_a() const { return acoeffs; } 00224 Vec<T2> get_coeffs_b() const { return bcoeffs; } 00226 void get_coeffs(Vec<T2> &b, Vec<T2> &a) const { b = bcoeffs; a = acoeffs; } 00228 void set_coeffs(const Vec<T2> &b, const Vec<T2> &a); 00230 void clear() { mem.clear(); } 00232 Vec<T3> get_state() const; 00234 void set_state(const Vec<T3> &state); 00235 00236 private: 00237 virtual T3 filter(const T1 Sample); 00238 00239 Vec<T3> mem; 00240 Vec<T2> acoeffs, bcoeffs; 00241 int inptr; 00242 bool init; 00243 }; 00244 00245 00246 00270 vec filter(const vec &b, const vec &a, const vec &input); 00271 cvec filter(const vec &b, const vec &a, const cvec &input); 00272 cvec filter(const cvec &b, const cvec &a, const cvec &input); 00273 cvec filter(const cvec &b, const cvec &a, const vec &input); 00274 00275 vec filter(const vec &b, const int one, const vec &input); 00276 cvec filter(const vec &b, const int one, const cvec &input); 00277 cvec filter(const cvec &b, const int one, const cvec &input); 00278 cvec filter(const cvec &b, const int one, const vec &input); 00279 00280 vec filter(const int one, const vec &a, const vec &input); 00281 cvec filter(const int one, const vec &a, const cvec &input); 00282 cvec filter(const int one, const cvec &a, const cvec &input); 00283 cvec filter(const int one, const cvec &a, const vec &input); 00284 00285 00286 vec filter(const vec &b, const vec &a, const vec &input, const vec &state_in, vec &state_out); 00287 cvec filter(const vec &b, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00288 cvec filter(const cvec &b, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00289 cvec filter(const cvec &b, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out); 00290 00291 vec filter(const vec &b, const int one, const vec &input, const vec &state_in, vec &state_out); 00292 cvec filter(const vec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out); 00293 cvec filter(const cvec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out); 00294 cvec filter(const cvec &b, const int one, const vec &input, const cvec &state_in, cvec &state_out); 00295 00296 vec filter(const int one, const vec &a, const vec &input, const vec &state_in, vec &state_out); 00297 cvec filter(const int one, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00298 cvec filter(const int one, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00299 cvec filter(const int one, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out); 00308 vec fir1(int N, double cutoff); 00309 00310 //---------------------------------------------------------------------------- 00311 // Implementation of templated functions starts here 00312 //---------------------------------------------------------------------------- 00313 00314 //---------------------- class Filter ---------------------------- 00315 00316 template <class T1, class T2, class T3> 00317 Vec<T3> Filter<T1,T2,T3>::operator()(const Vec<T1> &x) 00318 { 00319 Vec<T3> y(x.length()); 00320 00321 for (int i = 0; i < x.length(); i++) { 00322 y[i] = filter(x[i]); 00323 } 00324 00325 return y; 00326 } 00327 00328 //-------------------------- class MA_Filter --------------------------------- 00329 00330 template <class T1, class T2,class T3> 00331 MA_Filter<T1,T2,T3>::MA_Filter() : Filter<T1,T2,T3>() 00332 { 00333 inptr = 0; 00334 init = false; 00335 } 00336 00337 template <class T1, class T2,class T3> 00338 MA_Filter<T1,T2,T3>::MA_Filter(const Vec<T2> &b) : Filter<T1,T2,T3>() 00339 { 00340 set_coeffs(b); 00341 } 00342 00343 00344 template <class T1, class T2, class T3> 00345 void MA_Filter<T1,T2,T3>::set_coeffs(const Vec<T2> &b) 00346 { 00347 it_assert(b.size() > 0, "MA_Filter: size of filter is 0!"); 00348 00349 coeffs = b; 00350 mem.set_size(coeffs.size(), false); 00351 mem.clear(); 00352 inptr = 0; 00353 init = true; 00354 } 00355 00356 template <class T1, class T2, class T3> 00357 Vec<T3> MA_Filter<T1,T2,T3>::get_state() const 00358 { 00359 it_assert(init == true, "MA_Filter: filter coefficients are not set!"); 00360 00361 int offset = inptr; 00362 Vec<T3> state(mem.size()); 00363 00364 for(int n = 0; n < mem.size(); n++) { 00365 state(n) = mem(offset); 00366 offset = (offset + 1) % mem.size(); 00367 } 00368 00369 return state; 00370 } 00371 00372 template <class T1, class T2, class T3> 00373 void MA_Filter<T1,T2,T3>::set_state(const Vec<T3> &state) 00374 { 00375 it_assert(init == true, "MA_Filter: filter coefficients are not set!"); 00376 it_assert(state.size() == mem.size(), "MA_Filter: Invalid state vector!"); 00377 00378 mem = state; 00379 inptr = 0; 00380 } 00381 00382 template <class T1, class T2, class T3> 00383 T3 MA_Filter<T1,T2,T3>::filter(const T1 Sample) 00384 { 00385 it_assert(init == true, "MA_Filter: Filter coefficients are not set!"); 00386 T3 s = 0; 00387 00388 mem(inptr) = Sample; 00389 int L = mem.length() - inptr; 00390 00391 for (int i = 0; i < L; i++) { 00392 s += coeffs(i) * mem(inptr + i); 00393 } 00394 for (int i = 0; i < inptr; i++) { 00395 s += coeffs(L + i) * mem(i); 00396 } 00397 00398 inptr--; 00399 if (inptr < 0) 00400 inptr += mem.length(); 00401 00402 return s; 00403 } 00404 00405 //---------------------- class AR_Filter ---------------------------------- 00406 00407 template <class T1, class T2, class T3> 00408 AR_Filter<T1,T2,T3>::AR_Filter() : Filter<T1,T2,T3>() 00409 { 00410 inptr = 0; 00411 init = false; 00412 } 00413 00414 template <class T1, class T2, class T3> 00415 AR_Filter<T1,T2,T3>::AR_Filter(const Vec<T2> &a) : Filter<T1,T2,T3>() 00416 { 00417 set_coeffs(a); 00418 } 00419 00420 template <class T1, class T2, class T3> 00421 void AR_Filter<T1,T2,T3>::set_coeffs(const Vec<T2> &a) 00422 { 00423 it_assert(a.size() > 0, "AR_Filter: size of filter is 0!"); 00424 it_assert(a(0) != T2(0), "AR_Filter: a(0) cannot be 0!"); 00425 00426 a0 = a(0); 00427 coeffs = a / a0; 00428 00429 mem.set_size(coeffs.size() - 1, false); 00430 mem.clear(); 00431 inptr = 0; 00432 init = true; 00433 } 00434 00435 00436 template <class T1, class T2, class T3> 00437 Vec<T3> AR_Filter<T1,T2,T3>::get_state() const 00438 { 00439 it_assert(init == true, "AR_Filter: filter coefficients are not set!"); 00440 00441 int offset = inptr; 00442 Vec<T3> state(mem.size()); 00443 00444 for(int n = 0; n < mem.size(); n++) { 00445 state(n) = mem(offset); 00446 offset = (offset + 1) % mem.size(); 00447 } 00448 00449 return state; 00450 } 00451 00452 template <class T1, class T2, class T3> 00453 void AR_Filter<T1,T2,T3>::set_state(const Vec<T3> &state) 00454 { 00455 it_assert(init == true, "AR_Filter: filter coefficients are not set!"); 00456 it_assert(state.size() == mem.size(), "AR_Filter: Invalid state vector!"); 00457 00458 mem = state; 00459 inptr = 0; 00460 } 00461 00462 template <class T1, class T2, class T3> 00463 T3 AR_Filter<T1,T2,T3>::filter(const T1 Sample) 00464 { 00465 it_assert(init == true, "AR_Filter: Filter coefficients are not set!"); 00466 T3 s = Sample; 00467 00468 if (mem.size() == 0) 00469 return (s / a0); 00470 00471 int L = mem.size() - inptr; 00472 for (int i = 0; i < L; i++) { 00473 s -= mem(i + inptr) * coeffs(i + 1); // All coeffs except a(0) 00474 } 00475 for (int i = 0; i < inptr; i++) { 00476 s -= mem(i) * coeffs(L + i + 1); // All coeffs except a(0) 00477 } 00478 00479 inptr--; 00480 if (inptr < 0) 00481 inptr += mem.size(); 00482 mem(inptr) = s; 00483 00484 return (s / a0); 00485 } 00486 00487 00488 //---------------------- class ARMA_Filter ---------------------------------- 00489 template <class T1, class T2, class T3> 00490 ARMA_Filter<T1,T2,T3>::ARMA_Filter() : Filter<T1,T2,T3>() 00491 { 00492 inptr = 0; 00493 init = false; 00494 } 00495 00496 template <class T1, class T2, class T3> 00497 ARMA_Filter<T1,T2,T3>::ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a) : Filter<T1,T2,T3>() 00498 { 00499 set_coeffs(b, a); 00500 } 00501 00502 template <class T1, class T2, class T3> 00503 void ARMA_Filter<T1,T2,T3>::set_coeffs(const Vec<T2> &b, const Vec<T2> &a) 00504 { 00505 it_assert(a.size() > 0 && b.size() > 0, "ARMA_Filter: size of filter is 0!"); 00506 it_assert(a(0) != T2(0), "ARMA_Filter: a(0) cannot be 0!"); 00507 00508 acoeffs = a / a(0); 00509 bcoeffs = b / a(0); 00510 00511 mem.set_size(std::max(a.size(), b.size()) - 1, false); 00512 mem.clear(); 00513 inptr = 0; 00514 init = true; 00515 } 00516 00517 template <class T1, class T2, class T3> 00518 Vec<T3> ARMA_Filter<T1,T2,T3>::get_state() const 00519 { 00520 it_assert(init == true, "ARMA_Filter: filter coefficients are not set!"); 00521 00522 int offset = inptr; 00523 Vec<T3> state(mem.size()); 00524 00525 for(int n = 0; n < mem.size(); n++) { 00526 state(n) = mem(offset); 00527 offset = (offset + 1) % mem.size(); 00528 } 00529 00530 return state; 00531 } 00532 00533 template <class T1, class T2, class T3> 00534 void ARMA_Filter<T1,T2,T3>::set_state(const Vec<T3> &state) 00535 { 00536 it_assert(init == true, "ARMA_Filter: filter coefficients are not set!"); 00537 it_assert(state.size() == mem.size(), "ARMA_Filter: Invalid state vector!"); 00538 00539 mem = state; 00540 inptr = 0; 00541 } 00542 00543 template <class T1, class T2, class T3> 00544 T3 ARMA_Filter<T1,T2,T3>::filter(const T1 Sample) 00545 { 00546 it_assert(init == true, "ARMA_Filter: Filter coefficients are not set!"); 00547 T3 z = Sample; 00548 T3 s; 00549 00550 for(int i = 0; i < acoeffs.size() - 1; i++) { // All AR-coeff except a(0). 00551 z -= mem((i + inptr) % mem.size()) * acoeffs(i + 1); 00552 } 00553 s = z * bcoeffs(0); 00554 00555 for(int i = 0; i < bcoeffs.size() - 1; i++) { // All MA-coeff except b(0). 00556 s += mem((i + inptr) % mem.size()) * bcoeffs(i + 1); 00557 } 00558 00559 inptr--; 00560 if (inptr < 0) 00561 inptr += mem.size(); 00562 mem(inptr) = z; 00563 00564 mem(inptr) = z; // Store in the internal state. 00565 00566 return s; 00567 } 00568 00570 00571 // ---------------------------------------------------------------------- 00572 // Instantiations 00573 // ---------------------------------------------------------------------- 00574 00575 #ifdef HAVE_EXTERN_TEMPLATE 00576 00577 extern template class MA_Filter<double, double, double>; 00578 extern template class MA_Filter<double, std::complex<double>, 00579 std::complex<double> >; 00580 extern template class MA_Filter<std::complex<double>, double, 00581 std::complex<double> >; 00582 extern template class MA_Filter<std::complex<double>, std::complex<double>, 00583 std::complex<double> >; 00584 00585 extern template class AR_Filter<double, double, double>; 00586 extern template class AR_Filter<double, std::complex<double>, 00587 std::complex<double> >; 00588 extern template class AR_Filter<std::complex<double>, 00589 double,std::complex<double> >; 00590 extern template class AR_Filter<std::complex<double>, std::complex<double>, 00591 std::complex<double> >; 00592 00593 extern template class ARMA_Filter<double, double, double>; 00594 extern template class ARMA_Filter<double, std::complex<double>, 00595 std::complex<double> >; 00596 extern template class ARMA_Filter<std::complex<double>, 00597 double,std::complex<double> >; 00598 extern template class ARMA_Filter<std::complex<double>, std::complex<double>, 00599 std::complex<double> >; 00600 00601 #endif // HAVE_EXTERN_TEMPLATE 00602 00604 00605 } // namespace itpp 00606 00607 #endif // #ifndef FILTER_H
Generated on Sun Dec 9 17:30:27 2007 for IT++ by Doxygen 1.5.4