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

Generated on Sun Dec 9 17:30:27 2007 for IT++ by Doxygen 1.5.4