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