IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
filter.h
Go to the documentation of this file.
1 
29 #ifndef FILTER_H
30 #define FILTER_H
31 
32 #include <itpp/base/vec.h>
33 
34 
35 namespace itpp
36 {
37 
53 template <class T1, class T2, class T3>
54 class Filter
55 {
56 public:
58  Filter() {}
60  virtual T3 operator()(const T1 Sample) { return filter(Sample); }
62  virtual Vec<T3> operator()(const Vec<T1> &v);
64  virtual ~Filter() {}
65 protected:
70  virtual T3 filter(const T1 Sample) = 0;
71 };
72 
97 template <class T1, class T2, class T3>
98 class MA_Filter : public Filter<T1, T2, T3>
99 {
100 public:
102  explicit MA_Filter();
104  explicit MA_Filter(const Vec<T2> &b);
106  virtual ~MA_Filter() { }
108  Vec<T2> get_coeffs() const { return coeffs; }
110  void set_coeffs(const Vec<T2> &b);
112  void clear() { mem.clear(); }
114  Vec<T3> get_state() const;
116  void set_state(const Vec<T3> &state);
117 
118 private:
119  virtual T3 filter(const T1 Sample);
120 
121  Vec<T3> mem;
122  Vec<T2> coeffs;
123  int inptr;
124  bool init;
125 };
126 
151 template <class T1, class T2, class T3>
152 class AR_Filter : public Filter<T1, T2, T3>
153 {
154 public:
156  explicit AR_Filter();
158  explicit AR_Filter(const Vec<T2> &a);
160  virtual ~AR_Filter() { }
162  Vec<T2> get_coeffs() const { return coeffs; }
164  void set_coeffs(const Vec<T2> &a);
166  void clear() { mem.clear(); }
168  Vec<T3> get_state() const;
170  void set_state(const Vec<T3> &state);
171 
172 private:
173  virtual T3 filter(const T1 Sample);
174 
175  Vec<T3> mem;
176  Vec<T2> coeffs;
177  T2 a0;
178  int inptr;
179  bool init;
180 };
181 
182 
209 template <class T1, class T2, class T3>
210 class ARMA_Filter : public Filter<T1, T2, T3>
211 {
212 public:
214  explicit ARMA_Filter();
216  explicit ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a);
218  virtual ~ARMA_Filter() { }
220  Vec<T2> get_coeffs_a() const { return acoeffs; }
222  Vec<T2> get_coeffs_b() const { return bcoeffs; }
224  void get_coeffs(Vec<T2> &b, Vec<T2> &a) const { b = bcoeffs; a = acoeffs; }
226  void set_coeffs(const Vec<T2> &b, const Vec<T2> &a);
228  void clear() { mem.clear(); }
230  Vec<T3> get_state() const;
232  void set_state(const Vec<T3> &state);
233 
234 private:
235  virtual T3 filter(const T1 Sample);
236 
237  Vec<T3> mem;
238  Vec<T2> acoeffs, bcoeffs;
239  int inptr;
240  bool init;
241 };
242 
243 
244 
268 vec filter(const vec &b, const vec &a, const vec &input);
269 cvec filter(const vec &b, const vec &a, const cvec &input);
270 cvec filter(const cvec &b, const cvec &a, const cvec &input);
271 cvec filter(const cvec &b, const cvec &a, const vec &input);
272 
273 vec filter(const vec &b, const int one, const vec &input);
274 cvec filter(const vec &b, const int one, const cvec &input);
275 cvec filter(const cvec &b, const int one, const cvec &input);
276 cvec filter(const cvec &b, const int one, const vec &input);
277 
278 vec filter(const int one, const vec &a, const vec &input);
279 cvec filter(const int one, const vec &a, const cvec &input);
280 cvec filter(const int one, const cvec &a, const cvec &input);
281 cvec filter(const int one, const cvec &a, const vec &input);
282 
283 
284 vec filter(const vec &b, const vec &a, const vec &input, const vec &state_in, vec &state_out);
285 cvec filter(const vec &b, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out);
286 cvec filter(const cvec &b, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out);
287 cvec filter(const cvec &b, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out);
288 
289 vec filter(const vec &b, const int one, const vec &input, const vec &state_in, vec &state_out);
290 cvec filter(const vec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out);
291 cvec filter(const cvec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out);
292 cvec filter(const cvec &b, const int one, const vec &input, const cvec &state_in, cvec &state_out);
293 
294 vec filter(const int one, const vec &a, const vec &input, const vec &state_in, vec &state_out);
295 cvec filter(const int one, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out);
296 cvec filter(const int one, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out);
297 cvec filter(const int one, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out);
306 vec fir1(int N, double cutoff);
307 
308 //----------------------------------------------------------------------------
309 // Implementation of templated functions starts here
310 //----------------------------------------------------------------------------
311 
312 //---------------------- class Filter ----------------------------
313 
314 template <class T1, class T2, class T3>
316 {
317  Vec<T3> y(x.length());
318 
319  for (int i = 0; i < x.length(); i++) {
320  y[i] = filter(x[i]);
321  }
322 
323  return y;
324 }
325 
326 //-------------------------- class MA_Filter ---------------------------------
327 
328 template <class T1, class T2, class T3>
330 {
331  inptr = 0;
332  init = false;
333 }
334 
335 template <class T1, class T2, class T3>
337 {
338  set_coeffs(b);
339 }
340 
341 
342 template <class T1, class T2, class T3>
344 {
345  it_assert(b.size() > 0, "MA_Filter: size of filter is 0!");
346 
347  coeffs = b;
348  mem.set_size(coeffs.size(), false);
349  mem.clear();
350  inptr = 0;
351  init = true;
352 }
353 
354 template <class T1, class T2, class T3>
356 {
357  it_assert(init == true, "MA_Filter: filter coefficients are not set!");
358 
359  int offset = inptr;
360  Vec<T3> state(mem.size());
361 
362  for (int n = 0; n < mem.size(); n++) {
363  state(n) = mem(offset);
364  offset = (offset + 1) % mem.size();
365  }
366 
367  return state;
368 }
369 
370 template <class T1, class T2, class T3>
372 {
373  it_assert(init == true, "MA_Filter: filter coefficients are not set!");
374  it_assert(state.size() == mem.size(), "MA_Filter: Invalid state vector!");
375 
376  mem = state;
377  inptr = 0;
378 }
379 
380 template <class T1, class T2, class T3>
381 T3 MA_Filter<T1, T2, T3>::filter(const T1 Sample)
382 {
383  it_assert(init == true, "MA_Filter: Filter coefficients are not set!");
384  T3 s = 0;
385 
386  mem(inptr) = Sample;
387  int L = mem.length() - inptr;
388 
389  for (int i = 0; i < L; i++) {
390  s += coeffs(i) * mem(inptr + i);
391  }
392  for (int i = 0; i < inptr; i++) {
393  s += coeffs(L + i) * mem(i);
394  }
395 
396  inptr--;
397  if (inptr < 0)
398  inptr += mem.length();
399 
400  return s;
401 }
402 
403 //---------------------- class AR_Filter ----------------------------------
404 
405 template <class T1, class T2, class T3>
407 {
408  inptr = 0;
409  init = false;
410 }
411 
412 template <class T1, class T2, class T3>
414 {
415  set_coeffs(a);
416 }
417 
418 template <class T1, class T2, class T3>
420 {
421  it_assert(a.size() > 0, "AR_Filter: size of filter is 0!");
422  it_assert(a(0) != T2(0), "AR_Filter: a(0) cannot be 0!");
423 
424  a0 = a(0);
425  coeffs = a / a0;
426 
427  mem.set_size(coeffs.size() - 1, false);
428  mem.clear();
429  inptr = 0;
430  init = true;
431 }
432 
433 
434 template <class T1, class T2, class T3>
436 {
437  it_assert(init == true, "AR_Filter: filter coefficients are not set!");
438 
439  int offset = inptr;
440  Vec<T3> state(mem.size());
441 
442  for (int n = 0; n < mem.size(); n++) {
443  state(n) = mem(offset);
444  offset = (offset + 1) % mem.size();
445  }
446 
447  return state;
448 }
449 
450 template <class T1, class T2, class T3>
452 {
453  it_assert(init == true, "AR_Filter: filter coefficients are not set!");
454  it_assert(state.size() == mem.size(), "AR_Filter: Invalid state vector!");
455 
456  mem = state;
457  inptr = 0;
458 }
459 
460 template <class T1, class T2, class T3>
461 T3 AR_Filter<T1, T2, T3>::filter(const T1 Sample)
462 {
463  it_assert(init == true, "AR_Filter: Filter coefficients are not set!");
464  T3 s = Sample;
465 
466  if (mem.size() == 0)
467  return (s / a0);
468 
469  int L = mem.size() - inptr;
470  for (int i = 0; i < L; i++) {
471  s -= mem(i + inptr) * coeffs(i + 1); // All coeffs except a(0)
472  }
473  for (int i = 0; i < inptr; i++) {
474  s -= mem(i) * coeffs(L + i + 1); // All coeffs except a(0)
475  }
476 
477  inptr--;
478  if (inptr < 0)
479  inptr += mem.size();
480  mem(inptr) = s;
481 
482  return (s / a0);
483 }
484 
485 
486 //---------------------- class ARMA_Filter ----------------------------------
487 template <class T1, class T2, class T3>
489 {
490  inptr = 0;
491  init = false;
492 }
493 
494 template <class T1, class T2, class T3>
495 ARMA_Filter<T1, T2, T3>::ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a) : Filter<T1, T2, T3>()
496 {
497  set_coeffs(b, a);
498 }
499 
500 template <class T1, class T2, class T3>
502 {
503  it_assert(a.size() > 0 && b.size() > 0, "ARMA_Filter: size of filter is 0!");
504  it_assert(a(0) != T2(0), "ARMA_Filter: a(0) cannot be 0!");
505 
506  acoeffs = a / a(0);
507  bcoeffs = b / a(0);
508 
509  mem.set_size(std::max(a.size(), b.size()) - 1, false);
510  mem.clear();
511  inptr = 0;
512  init = true;
513 }
514 
515 template <class T1, class T2, class T3>
517 {
518  it_assert(init == true, "ARMA_Filter: filter coefficients are not set!");
519 
520  int offset = inptr;
521  Vec<T3> state(mem.size());
522 
523  for (int n = 0; n < mem.size(); n++) {
524  state(n) = mem(offset);
525  offset = (offset + 1) % mem.size();
526  }
527 
528  return state;
529 }
530 
531 template <class T1, class T2, class T3>
533 {
534  it_assert(init == true, "ARMA_Filter: filter coefficients are not set!");
535  it_assert(state.size() == mem.size(), "ARMA_Filter: Invalid state vector!");
536 
537  mem = state;
538  inptr = 0;
539 }
540 
541 template <class T1, class T2, class T3>
542 T3 ARMA_Filter<T1, T2, T3>::filter(const T1 Sample)
543 {
544  it_assert(init == true, "ARMA_Filter: Filter coefficients are not set!");
545  T3 z = Sample;
546  T3 s;
547 
548  for (int i = 0; i < acoeffs.size() - 1; i++) { // All AR-coeff except a(0).
549  z -= mem((i + inptr) % mem.size()) * acoeffs(i + 1);
550  }
551  s = z * bcoeffs(0);
552 
553  for (int i = 0; i < bcoeffs.size() - 1; i++) { // All MA-coeff except b(0).
554  s += mem((i + inptr) % mem.size()) * bcoeffs(i + 1);
555  }
556 
557  inptr--;
558  if (inptr < 0)
559  inptr += mem.size();
560  mem(inptr) = z;
561 
562  mem(inptr) = z; // Store in the internal state.
563 
564  return s;
565 }
566 
568 
569 // ----------------------------------------------------------------------
570 // Instantiations
571 // ----------------------------------------------------------------------
572 
573 #ifndef _MSC_VER
574 
575 extern template class MA_Filter<double, double, double>;
576 extern template class MA_Filter < double, std::complex<double>,
577  std::complex<double> >;
578 extern template class MA_Filter < std::complex<double>, double,
579  std::complex<double> >;
580 extern template class MA_Filter < std::complex<double>, std::complex<double>,
581  std::complex<double> >;
582 
583 extern template class AR_Filter<double, double, double>;
584 extern template class AR_Filter < double, std::complex<double>,
585  std::complex<double> >;
586 extern template class AR_Filter < std::complex<double>,
587  double, std::complex<double> >;
588 extern template class AR_Filter < std::complex<double>, std::complex<double>,
589  std::complex<double> >;
590 
591 extern template class ARMA_Filter<double, double, double>;
592 extern template class ARMA_Filter < double, std::complex<double>,
593  std::complex<double> >;
594 extern template class ARMA_Filter < std::complex<double>,
595  double, std::complex<double> >;
596 extern template class ARMA_Filter < std::complex<double>, std::complex<double>,
597  std::complex<double> >;
598 
599 #endif // _MSC_VER
600 
602 
603 } // namespace itpp
604 
605 #endif // #ifndef FILTER_H
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:14 for IT++ by Doxygen 1.8.1.2