IT++ Logo

filter.h

Go to the documentation of this file.
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 {
00044 
00060 template <class T1, class T2, class T3>
00061 class Filter
00062 {
00063 public:
00065   Filter() {}
00067   virtual T3 operator()(const T1 Sample) { return filter(Sample); }
00069   virtual Vec<T3> operator()(const Vec<T1> &v);
00071   virtual ~Filter() {}
00072 protected:
00077   virtual T3 filter(const T1 Sample) = 0;
00078 };
00079 
00104 template <class T1, class T2, class T3>
00105 class MA_Filter : public Filter<T1, T2, T3>
00106 {
00107 public:
00109   explicit MA_Filter();
00111   explicit MA_Filter(const Vec<T2> &b);
00113   virtual ~MA_Filter() { }
00115   Vec<T2> get_coeffs() const { return coeffs; }
00117   void set_coeffs(const Vec<T2> &b);
00119   void clear() { mem.clear(); }
00121   Vec<T3> get_state() const;
00123   void set_state(const Vec<T3> &state);
00124 
00125 private:
00126   virtual T3 filter(const T1 Sample);
00127 
00128   Vec<T3> mem;
00129   Vec<T2> coeffs;
00130   int inptr;
00131   bool init;
00132 };
00133 
00158 template <class T1, class T2, class T3>
00159 class AR_Filter : public Filter<T1, T2, T3>
00160 {
00161 public:
00163   explicit AR_Filter();
00165   explicit AR_Filter(const Vec<T2> &a);
00167   virtual ~AR_Filter() { }
00169   Vec<T2> get_coeffs() const { return coeffs; }
00171   void set_coeffs(const Vec<T2> &a);
00173   void clear() { mem.clear(); }
00175   Vec<T3> get_state() const;
00177   void set_state(const Vec<T3> &state);
00178 
00179 private:
00180   virtual T3 filter(const T1 Sample);
00181 
00182   Vec<T3> mem;
00183   Vec<T2> coeffs;
00184   T2 a0;
00185   int inptr;
00186   bool init;
00187 };
00188 
00189 
00216 template <class T1, class T2, class T3>
00217 class ARMA_Filter : public Filter<T1, T2, T3>
00218 {
00219 public:
00221   explicit ARMA_Filter();
00223   explicit ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a);
00225   virtual ~ARMA_Filter() { }
00227   Vec<T2> get_coeffs_a() const { return acoeffs; }
00229   Vec<T2> get_coeffs_b() const { return bcoeffs; }
00231   void get_coeffs(Vec<T2> &b, Vec<T2> &a) const { b = bcoeffs; a = acoeffs; }
00233   void set_coeffs(const Vec<T2> &b, const Vec<T2> &a);
00235   void clear() { mem.clear(); }
00237   Vec<T3> get_state() const;
00239   void set_state(const Vec<T3> &state);
00240 
00241 private:
00242   virtual T3 filter(const T1 Sample);
00243 
00244   Vec<T3> mem;
00245   Vec<T2> acoeffs, bcoeffs;
00246   int inptr;
00247   bool init;
00248 };
00249 
00250 
00251 
00275 vec filter(const vec &b, const vec &a, const vec &input);
00276 cvec filter(const vec &b, const vec &a, const cvec &input);
00277 cvec filter(const cvec &b, const cvec &a, const cvec &input);
00278 cvec filter(const cvec &b, const cvec &a, const vec &input);
00279 
00280 vec filter(const vec &b, const int one, const vec &input);
00281 cvec filter(const vec &b, const int one, const cvec &input);
00282 cvec filter(const cvec &b, const int one, const cvec &input);
00283 cvec filter(const cvec &b, const int one, const vec &input);
00284 
00285 vec filter(const int one, const vec &a, const vec &input);
00286 cvec filter(const int one, const vec &a, const cvec &input);
00287 cvec filter(const int one, const cvec &a, const cvec &input);
00288 cvec filter(const int one, const cvec &a, const vec &input);
00289 
00290 
00291 vec filter(const vec &b, const vec &a, const vec &input, const vec &state_in, vec &state_out);
00292 cvec filter(const vec &b, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out);
00293 cvec filter(const cvec &b, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out);
00294 cvec filter(const cvec &b, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out);
00295 
00296 vec filter(const vec &b, const int one, const vec &input, const vec &state_in, vec &state_out);
00297 cvec filter(const vec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out);
00298 cvec filter(const cvec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out);
00299 cvec filter(const cvec &b, const int one, const vec &input, const cvec &state_in, cvec &state_out);
00300 
00301 vec filter(const int one, const vec &a, const vec &input, const vec &state_in, vec &state_out);
00302 cvec filter(const int one, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out);
00303 cvec filter(const int one, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out);
00304 cvec filter(const int one, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out);
00313 vec fir1(int N, double cutoff);
00314 
00315 //----------------------------------------------------------------------------
00316 // Implementation of templated functions starts here
00317 //----------------------------------------------------------------------------
00318 
00319 //---------------------- class Filter ----------------------------
00320 
00321 template <class T1, class T2, class T3>
00322 Vec<T3> Filter<T1, T2, T3>::operator()(const Vec<T1> &x)
00323 {
00324   Vec<T3> y(x.length());
00325 
00326   for (int i = 0; i < x.length(); i++) {
00327     y[i] = filter(x[i]);
00328   }
00329 
00330   return y;
00331 }
00332 
00333 //-------------------------- class MA_Filter ---------------------------------
00334 
00335 template <class T1, class T2, class T3>
00336 MA_Filter<T1, T2, T3>::MA_Filter() : Filter<T1, T2, T3>()
00337 {
00338   inptr = 0;
00339   init = false;
00340 }
00341 
00342 template <class T1, class T2, class T3>
00343 MA_Filter<T1, T2, T3>::MA_Filter(const Vec<T2> &b) : Filter<T1, T2, T3>()
00344 {
00345   set_coeffs(b);
00346 }
00347 
00348 
00349 template <class T1, class T2, class T3>
00350 void MA_Filter<T1, T2, T3>::set_coeffs(const Vec<T2> &b)
00351 {
00352   it_assert(b.size() > 0, "MA_Filter: size of filter is 0!");
00353 
00354   coeffs = b;
00355   mem.set_size(coeffs.size(), false);
00356   mem.clear();
00357   inptr = 0;
00358   init = true;
00359 }
00360 
00361 template <class T1, class T2, class T3>
00362 Vec<T3> MA_Filter<T1, T2, T3>::get_state() const
00363 {
00364   it_assert(init == true, "MA_Filter: filter coefficients are not set!");
00365 
00366   int offset = inptr;
00367   Vec<T3> state(mem.size());
00368 
00369   for (int n = 0; n < mem.size(); n++) {
00370     state(n) = mem(offset);
00371     offset = (offset + 1) % mem.size();
00372   }
00373 
00374   return state;
00375 }
00376 
00377 template <class T1, class T2, class T3>
00378 void MA_Filter<T1, T2, T3>::set_state(const Vec<T3> &state)
00379 {
00380   it_assert(init == true, "MA_Filter: filter coefficients are not set!");
00381   it_assert(state.size() == mem.size(), "MA_Filter: Invalid state vector!");
00382 
00383   mem = state;
00384   inptr = 0;
00385 }
00386 
00387 template <class T1, class T2, class T3>
00388 T3 MA_Filter<T1, T2, T3>::filter(const T1 Sample)
00389 {
00390   it_assert(init == true, "MA_Filter: Filter coefficients are not set!");
00391   T3 s = 0;
00392 
00393   mem(inptr) = Sample;
00394   int L = mem.length() - inptr;
00395 
00396   for (int i = 0; i < L; i++) {
00397     s += coeffs(i) * mem(inptr + i);
00398   }
00399   for (int i = 0; i < inptr; i++) {
00400     s += coeffs(L + i) * mem(i);
00401   }
00402 
00403   inptr--;
00404   if (inptr < 0)
00405     inptr += mem.length();
00406 
00407   return s;
00408 }
00409 
00410 //---------------------- class AR_Filter ----------------------------------
00411 
00412 template <class T1, class T2, class T3>
00413 AR_Filter<T1, T2, T3>::AR_Filter() : Filter<T1, T2, T3>()
00414 {
00415   inptr = 0;
00416   init = false;
00417 }
00418 
00419 template <class T1, class T2, class T3>
00420 AR_Filter<T1, T2, T3>::AR_Filter(const Vec<T2> &a) : Filter<T1, T2, T3>()
00421 {
00422   set_coeffs(a);
00423 }
00424 
00425 template <class T1, class T2, class T3>
00426 void AR_Filter<T1, T2, T3>::set_coeffs(const Vec<T2> &a)
00427 {
00428   it_assert(a.size() > 0, "AR_Filter: size of filter is 0!");
00429   it_assert(a(0) != T2(0), "AR_Filter: a(0) cannot be 0!");
00430 
00431   a0 = a(0);
00432   coeffs = a / a0;
00433 
00434   mem.set_size(coeffs.size() - 1, false);
00435   mem.clear();
00436   inptr = 0;
00437   init = true;
00438 }
00439 
00440 
00441 template <class T1, class T2, class T3>
00442 Vec<T3> AR_Filter<T1, T2, T3>::get_state() const
00443 {
00444   it_assert(init == true, "AR_Filter: filter coefficients are not set!");
00445 
00446   int offset = inptr;
00447   Vec<T3> state(mem.size());
00448 
00449   for (int n = 0; n < mem.size(); n++) {
00450     state(n) = mem(offset);
00451     offset = (offset + 1) % mem.size();
00452   }
00453 
00454   return state;
00455 }
00456 
00457 template <class T1, class T2, class T3>
00458 void AR_Filter<T1, T2, T3>::set_state(const Vec<T3> &state)
00459 {
00460   it_assert(init == true, "AR_Filter: filter coefficients are not set!");
00461   it_assert(state.size() == mem.size(), "AR_Filter: Invalid state vector!");
00462 
00463   mem = state;
00464   inptr = 0;
00465 }
00466 
00467 template <class T1, class T2, class T3>
00468 T3 AR_Filter<T1, T2, T3>::filter(const T1 Sample)
00469 {
00470   it_assert(init == true, "AR_Filter: Filter coefficients are not set!");
00471   T3 s = Sample;
00472 
00473   if (mem.size() == 0)
00474     return (s / a0);
00475 
00476   int L = mem.size() - inptr;
00477   for (int i = 0; i < L; i++) {
00478     s -= mem(i + inptr) * coeffs(i + 1); // All coeffs except a(0)
00479   }
00480   for (int i = 0; i < inptr; i++) {
00481     s -= mem(i) * coeffs(L + i + 1); // All coeffs except a(0)
00482   }
00483 
00484   inptr--;
00485   if (inptr < 0)
00486     inptr += mem.size();
00487   mem(inptr) = s;
00488 
00489   return (s / a0);
00490 }
00491 
00492 
00493 //---------------------- class ARMA_Filter ----------------------------------
00494 template <class T1, class T2, class T3>
00495 ARMA_Filter<T1, T2, T3>::ARMA_Filter() : Filter<T1, T2, T3>()
00496 {
00497   inptr = 0;
00498   init = false;
00499 }
00500 
00501 template <class T1, class T2, class T3>
00502 ARMA_Filter<T1, T2, T3>::ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a) : Filter<T1, T2, T3>()
00503 {
00504   set_coeffs(b, a);
00505 }
00506 
00507 template <class T1, class T2, class T3>
00508 void ARMA_Filter<T1, T2, T3>::set_coeffs(const Vec<T2> &b, const Vec<T2> &a)
00509 {
00510   it_assert(a.size() > 0 && b.size() > 0, "ARMA_Filter: size of filter is 0!");
00511   it_assert(a(0) != T2(0), "ARMA_Filter: a(0) cannot be 0!");
00512 
00513   acoeffs = a / a(0);
00514   bcoeffs = b / a(0);
00515 
00516   mem.set_size(std::max(a.size(), b.size()) - 1, false);
00517   mem.clear();
00518   inptr = 0;
00519   init = true;
00520 }
00521 
00522 template <class T1, class T2, class T3>
00523 Vec<T3> ARMA_Filter<T1, T2, T3>::get_state() const
00524 {
00525   it_assert(init == true, "ARMA_Filter: filter coefficients are not set!");
00526 
00527   int offset = inptr;
00528   Vec<T3> state(mem.size());
00529 
00530   for (int n = 0; n < mem.size(); n++) {
00531     state(n) = mem(offset);
00532     offset = (offset + 1) % mem.size();
00533   }
00534 
00535   return state;
00536 }
00537 
00538 template <class T1, class T2, class T3>
00539 void ARMA_Filter<T1, T2, T3>::set_state(const Vec<T3> &state)
00540 {
00541   it_assert(init == true, "ARMA_Filter: filter coefficients are not set!");
00542   it_assert(state.size() == mem.size(), "ARMA_Filter: Invalid state vector!");
00543 
00544   mem = state;
00545   inptr = 0;
00546 }
00547 
00548 template <class T1, class T2, class T3>
00549 T3 ARMA_Filter<T1, T2, T3>::filter(const T1 Sample)
00550 {
00551   it_assert(init == true, "ARMA_Filter: Filter coefficients are not set!");
00552   T3 z = Sample;
00553   T3 s;
00554 
00555   for (int i = 0; i < acoeffs.size() - 1; i++) { // All AR-coeff except a(0).
00556     z -= mem((i + inptr) % mem.size()) * acoeffs(i + 1);
00557   }
00558   s = z * bcoeffs(0);
00559 
00560   for (int i = 0; i < bcoeffs.size() - 1; i++) { // All MA-coeff except b(0).
00561     s += mem((i + inptr) % mem.size()) * bcoeffs(i + 1);
00562   }
00563 
00564   inptr--;
00565   if (inptr < 0)
00566     inptr += mem.size();
00567   mem(inptr) = z;
00568 
00569   mem(inptr) = z; // Store in the internal state.
00570 
00571   return s;
00572 }
00573 
00575 
00576 // ----------------------------------------------------------------------
00577 // Instantiations
00578 // ----------------------------------------------------------------------
00579 
00580 #ifdef HAVE_EXTERN_TEMPLATE
00581 
00582 extern template class MA_Filter<double, double, double>;
00583 extern template class MA_Filter < double, std::complex<double>,
00584   std::complex<double> >;
00585 extern template class MA_Filter < std::complex<double>, double,
00586   std::complex<double> >;
00587 extern template class MA_Filter < std::complex<double>, std::complex<double>,
00588   std::complex<double> >;
00589 
00590 extern template class AR_Filter<double, double, double>;
00591 extern template class AR_Filter < double, std::complex<double>,
00592   std::complex<double> >;
00593 extern template class AR_Filter < std::complex<double>,
00594   double, std::complex<double> >;
00595 extern template class AR_Filter < std::complex<double>, std::complex<double>,
00596   std::complex<double> >;
00597 
00598 extern template class ARMA_Filter<double, double, double>;
00599 extern template class ARMA_Filter < double, std::complex<double>,
00600   std::complex<double> >;
00601 extern template class ARMA_Filter < std::complex<double>,
00602   double, std::complex<double> >;
00603 extern template class ARMA_Filter < std::complex<double>, std::complex<double>,
00604   std::complex<double> >;
00605 
00606 #endif // HAVE_EXTERN_TEMPLATE
00607 
00609 
00610 } // namespace itpp
00611 
00612 #endif // #ifndef FILTER_H
SourceForge Logo

Generated on Sun Jul 26 08:54:30 2009 for IT++ by Doxygen 1.5.9