00001 00033 #ifndef FILTER_H 00034 #define FILTER_H 00035 00036 #include <itpp/base/vec.h> 00037 00038 00039 namespace itpp { 00040 00056 template <class T1, class T2, class T3> 00057 class Filter { 00058 public: 00060 Filter() {} 00062 virtual T3 operator()(const T1 Sample) { return filter(Sample); } 00064 virtual Vec<T3> operator()(const Vec<T1> &v); 00066 virtual ~Filter() {} 00067 protected: 00072 virtual T3 filter(const T1 Sample)=0; 00073 }; 00074 00099 template <class T1, class T2, class T3> 00100 class MA_Filter : public Filter<T1,T2,T3> { 00101 public: 00103 explicit MA_Filter(); 00105 explicit MA_Filter(const Vec<T2> &b); 00107 virtual ~MA_Filter() { } 00109 Vec<T2> get_coeffs() const { return coeffs; } 00111 void set_coeffs(const Vec<T2> &b); 00113 void clear() { mem.clear(); } 00115 Vec<T3> get_state() const; 00117 void set_state(const Vec<T3> &state); 00118 00119 private: 00120 virtual T3 filter(const T1 Sample); 00121 00122 Vec<T3> mem; 00123 Vec<T2> coeffs; 00124 int inptr; 00125 bool init; 00126 }; 00127 00152 template <class T1, class T2, class T3> 00153 class AR_Filter : public Filter<T1,T2,T3> { 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 public: 00213 explicit ARMA_Filter(); 00215 explicit ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a); 00217 virtual ~ARMA_Filter() { } 00219 Vec<T2> get_coeffs_a() const { return acoeffs; } 00221 Vec<T2> get_coeffs_b() const { return bcoeffs; } 00223 void get_coeffs(Vec<T2> &b, Vec<T2> &a) const { b = bcoeffs; a = acoeffs; } 00225 void set_coeffs(const Vec<T2> &b, const Vec<T2> &a); 00227 void clear() { mem.clear(); } 00229 Vec<T3> get_state() const; 00231 void set_state(const Vec<T3> &state); 00232 00233 private: 00234 virtual T3 filter(const T1 Sample); 00235 00236 Vec<T3> mem; 00237 Vec<T2> acoeffs, bcoeffs; 00238 int inptr; 00239 bool init; 00240 }; 00241 00242 00243 00266 vec filter(const vec &b, const vec &a, const vec &input); 00267 cvec filter(const vec &b, const vec &a, const cvec &input); 00268 cvec filter(const cvec &b, const cvec &a, const cvec &input); 00269 cvec filter(const cvec &b, const cvec &a, const vec &input); 00270 00271 vec filter(const vec &b, const int one, const vec &input); 00272 cvec filter(const vec &b, const int one, const cvec &input); 00273 cvec filter(const cvec &b, const int one, const cvec &input); 00274 cvec filter(const cvec &b, const int one, const vec &input); 00275 00276 vec filter(const int one, const vec &a, const vec &input); 00277 cvec filter(const int one, const vec &a, const cvec &input); 00278 cvec filter(const int one, const cvec &a, const cvec &input); 00279 cvec filter(const int one, const cvec &a, const vec &input); 00280 00281 00282 vec filter(const vec &b, const vec &a, const vec &input, const vec &state_in, vec &state_out); 00283 cvec filter(const vec &b, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00284 cvec filter(const cvec &b, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00285 cvec filter(const cvec &b, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out); 00286 00287 vec filter(const vec &b, const int one, const vec &input, const vec &state_in, vec &state_out); 00288 cvec filter(const vec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out); 00289 cvec filter(const cvec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out); 00290 cvec filter(const cvec &b, const int one, const vec &input, const cvec &state_in, cvec &state_out); 00291 00292 vec filter(const int one, const vec &a, const vec &input, const vec &state_in, vec &state_out); 00293 cvec filter(const int one, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00294 cvec filter(const int one, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out); 00295 cvec filter(const int one, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out); 00297 00298 00304 vec fir1(int N, double cutoff); 00305 00306 //---------------------------------------------------------------------------- 00307 // Implementation of templated functions starts here 00308 //---------------------------------------------------------------------------- 00309 00310 //---------------------- class Filter ---------------------------- 00311 00312 template <class T1, class T2, class T3> 00313 Vec<T3> Filter<T1,T2,T3>::operator()(const Vec<T1> &x) 00314 { 00315 Vec<T3> y(x.length()); 00316 00317 for (int i = 0; i < x.length(); i++) { 00318 y[i] = filter(x[i]); 00319 } 00320 00321 return y; 00322 } 00323 00324 //-------------------------- class MA_Filter --------------------------------- 00325 00326 template <class T1, class T2,class T3> 00327 MA_Filter<T1,T2,T3>::MA_Filter() : Filter<T1,T2,T3>() 00328 { 00329 inptr = 0; 00330 init = false; 00331 } 00332 00333 template <class T1, class T2,class T3> 00334 MA_Filter<T1,T2,T3>::MA_Filter(const Vec<T2> &b) : Filter<T1,T2,T3>() 00335 { 00336 set_coeffs(b); 00337 } 00338 00339 00340 template <class T1, class T2, class T3> 00341 void MA_Filter<T1,T2,T3>::set_coeffs(const Vec<T2> &b) 00342 { 00343 it_assert(b.size() > 0, "MA_Filter: size of filter is 0!"); 00344 00345 coeffs = b; 00346 mem.set_size(coeffs.size(), false); 00347 mem.clear(); 00348 inptr = 0; 00349 init = true; 00350 } 00351 00352 template <class T1, class T2, class T3> 00353 Vec<T3> MA_Filter<T1,T2,T3>::get_state() const 00354 { 00355 it_assert(init == true, "MA_Filter: filter coefficients are not set!"); 00356 00357 int offset = inptr; 00358 Vec<T3> state(mem.size()); 00359 00360 for(int n = 0; n < mem.size(); n++) { 00361 state(n) = mem(offset); 00362 offset = (offset + 1) % mem.size(); 00363 } 00364 00365 return state; 00366 } 00367 00368 template <class T1, class T2, class T3> 00369 void MA_Filter<T1,T2,T3>::set_state(const Vec<T3> &state) 00370 { 00371 it_assert(init == true, "MA_Filter: filter coefficients are not set!"); 00372 it_assert(state.size() == mem.size(), "MA_Filter: Invalid state vector!"); 00373 00374 mem = state; 00375 inptr = 0; 00376 } 00377 00378 template <class T1, class T2, class T3> 00379 T3 MA_Filter<T1,T2,T3>::filter(const T1 Sample) 00380 { 00381 it_assert(init == true, "MA_Filter: Filter coefficients are not set!"); 00382 T3 s = 0; 00383 00384 mem(inptr) = Sample; 00385 int L = mem.length() - inptr; 00386 00387 for (int i = 0; i < L; i++) { 00388 s += coeffs(i) * mem(inptr + i); 00389 } 00390 for (int i = 0; i < inptr; i++) { 00391 s += coeffs(L + i) * mem(i); 00392 } 00393 00394 inptr--; 00395 if (inptr < 0) 00396 inptr += mem.length(); 00397 00398 return s; 00399 } 00400 00401 //---------------------- class AR_Filter ---------------------------------- 00402 00403 template <class T1, class T2, class T3> 00404 AR_Filter<T1,T2,T3>::AR_Filter() : Filter<T1,T2,T3>() 00405 { 00406 inptr = 0; 00407 init = false; 00408 } 00409 00410 template <class T1, class T2, class T3> 00411 AR_Filter<T1,T2,T3>::AR_Filter(const Vec<T2> &a) : Filter<T1,T2,T3>() 00412 { 00413 set_coeffs(a); 00414 } 00415 00416 template <class T1, class T2, class T3> 00417 void AR_Filter<T1,T2,T3>::set_coeffs(const Vec<T2> &a) 00418 { 00419 it_assert(a.size() > 0, "AR_Filter: size of filter is 0!"); 00420 it_assert(a(0) != T2(0), "AR_Filter: a(0) cannot be 0!"); 00421 00422 a0 = a(0); 00423 coeffs = a / a0; 00424 00425 mem.set_size(coeffs.size() - 1, false); 00426 mem.clear(); 00427 inptr = 0; 00428 init = true; 00429 } 00430 00431 00432 template <class T1, class T2, class T3> 00433 Vec<T3> AR_Filter<T1,T2,T3>::get_state() const 00434 { 00435 it_assert(init == true, "AR_Filter: filter coefficients are not set!"); 00436 00437 int offset = inptr; 00438 Vec<T3> state(mem.size()); 00439 00440 for(int n = 0; n < mem.size(); n++) { 00441 state(n) = mem(offset); 00442 offset = (offset + 1) % mem.size(); 00443 } 00444 00445 return state; 00446 } 00447 00448 template <class T1, class T2, class T3> 00449 void AR_Filter<T1,T2,T3>::set_state(const Vec<T3> &state) 00450 { 00451 it_assert(init == true, "AR_Filter: filter coefficients are not set!"); 00452 it_assert(state.size() == mem.size(), "AR_Filter: Invalid state vector!"); 00453 00454 mem = state; 00455 inptr = 0; 00456 } 00457 00458 template <class T1, class T2, class T3> 00459 T3 AR_Filter<T1,T2,T3>::filter(const T1 Sample) 00460 { 00461 it_assert(init == true, "AR_Filter: Filter coefficients are not set!"); 00462 T3 s = Sample; 00463 00464 if (mem.size() == 0) 00465 return (s / a0); 00466 00467 int L = mem.size() - inptr; 00468 for (int i = 0; i < L; i++) { 00469 s -= mem(i + inptr) * coeffs(i + 1); // All coeffs except a(0) 00470 } 00471 for (int i = 0; i < inptr; i++) { 00472 s -= mem(i) * coeffs(L + i + 1); // All coeffs except a(0) 00473 } 00474 00475 inptr--; 00476 if (inptr < 0) 00477 inptr += mem.size(); 00478 mem(inptr) = s; 00479 00480 return (s / a0); 00481 } 00482 00483 00484 //---------------------- class ARMA_Filter ---------------------------------- 00485 template <class T1, class T2, class T3> 00486 ARMA_Filter<T1,T2,T3>::ARMA_Filter() : Filter<T1,T2,T3>() 00487 { 00488 inptr = 0; 00489 init = false; 00490 } 00491 00492 template <class T1, class T2, class T3> 00493 ARMA_Filter<T1,T2,T3>::ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a) : Filter<T1,T2,T3>() 00494 { 00495 set_coeffs(b, a); 00496 } 00497 00498 template <class T1, class T2, class T3> 00499 void ARMA_Filter<T1,T2,T3>::set_coeffs(const Vec<T2> &b, const Vec<T2> &a) 00500 { 00501 it_assert(a.size() > 0 && b.size() > 0, "ARMA_Filter: size of filter is 0!"); 00502 it_assert(a(0) != T2(0), "ARMA_Filter: a(0) cannot be 0!"); 00503 00504 acoeffs = a / a(0); 00505 bcoeffs = b / a(0); 00506 00507 mem.set_size(std::max(a.size(), b.size()) - 1, false); 00508 mem.clear(); 00509 inptr = 0; 00510 init = true; 00511 } 00512 00513 template <class T1, class T2, class T3> 00514 Vec<T3> ARMA_Filter<T1,T2,T3>::get_state() const 00515 { 00516 it_assert(init == true, "ARMA_Filter: filter coefficients are not set!"); 00517 00518 int offset = inptr; 00519 Vec<T3> state(mem.size()); 00520 00521 for(int n = 0; n < mem.size(); n++) { 00522 state(n) = mem(offset); 00523 offset = (offset + 1) % mem.size(); 00524 } 00525 00526 return state; 00527 } 00528 00529 template <class T1, class T2, class T3> 00530 void ARMA_Filter<T1,T2,T3>::set_state(const Vec<T3> &state) 00531 { 00532 it_assert(init == true, "ARMA_Filter: filter coefficients are not set!"); 00533 it_assert(state.size() == mem.size(), "ARMA_Filter: Invalid state vector!"); 00534 00535 mem = state; 00536 inptr = 0; 00537 } 00538 00539 template <class T1, class T2, class T3> 00540 T3 ARMA_Filter<T1,T2,T3>::filter(const T1 Sample) 00541 { 00542 it_assert(init == true, "ARMA_Filter: Filter coefficients are not set!"); 00543 T3 z = Sample; 00544 T3 s; 00545 00546 for(int i = 0; i < acoeffs.size() - 1; i++) { // All AR-coeff except a(0). 00547 z -= mem((i + inptr) % mem.size()) * acoeffs(i + 1); 00548 } 00549 s = z * bcoeffs(0); 00550 00551 for(int i = 0; i < bcoeffs.size() - 1; i++) { // All MA-coeff except b(0). 00552 s += mem((i + inptr) % mem.size()) * bcoeffs(i + 1); 00553 } 00554 00555 inptr--; 00556 if (inptr < 0) 00557 inptr += mem.size(); 00558 mem(inptr) = z; 00559 00560 mem(inptr) = z; // Store in the internal state. 00561 00562 return s; 00563 } 00564 00565 00566 #ifndef _MSC_VER 00567 00568 //----------------------------------------------------------------------- 00569 // class MA_Filter 00570 //----------------------------------------------------------------------- 00571 00573 extern template class MA_Filter<double,double,double>; 00575 extern template class MA_Filter<double,std::complex<double>,std::complex<double> >; 00577 extern template class MA_Filter<std::complex<double>,double,std::complex<double> >; 00579 extern template class MA_Filter<std::complex<double>,std::complex<double>,std::complex<double> >; 00580 00581 //----------------------------------------------------------------------- 00582 // class AR_Filter 00583 //----------------------------------------------------------------------- 00584 00586 extern template class AR_Filter<double,double,double>; 00588 extern template class AR_Filter<double,std::complex<double>,std::complex<double> >; 00590 extern template class AR_Filter<std::complex<double>,double,std::complex<double> >; 00592 extern template class AR_Filter<std::complex<double>,std::complex<double>,std::complex<double> >; 00593 00594 //----------------------------------------------------------------------- 00595 // class ARMA_Filter 00596 //----------------------------------------------------------------------- 00597 00599 extern template class ARMA_Filter<double,double,double>; 00601 extern template class ARMA_Filter<double,std::complex<double>,std::complex<double> >; 00603 extern template class ARMA_Filter<std::complex<double>,double,std::complex<double> >; 00605 extern template class ARMA_Filter<std::complex<double>,std::complex<double>,std::complex<double> >; 00606 00607 #endif // MSC_VER 00608 00609 } // namespace itpp 00610 00611 #endif // #ifndef FILTER_H
Generated on Fri Jan 11 08:51:36 2008 for IT++ by Doxygen 1.3.9.1