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