IT++ Logo Newcom Logo

filter.h

Go to the documentation of this file.
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
SourceForge Logo

Generated on Sat Aug 25 23:37:27 2007 for IT++ by Doxygen 1.5.2