IT++ Logo
filter.h
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4