filter.h

Go to the documentation of this file.

#ifndef FILTER_H
#define FILTER_H

#include <itpp/base/vec.h>
#include <itpp/itexports.h>

namespace itpp
{

template <class T1, class T2, class T3>
class Filter
{
public:
  Filter() {}
  virtual T3 operator()(const T1 Sample) { return filter(Sample); }
  virtual Vec<T3> operator()(const Vec<T1> &v);
  virtual ~Filter() {}
protected:
  virtual T3 filter(const T1 Sample) = 0;
};

template <class T1, class T2, class T3>
class MA_Filter : public Filter<T1, T2, T3>
{
public:
  explicit MA_Filter();
  explicit MA_Filter(const Vec<T2> &b);
  virtual ~MA_Filter() { }
  Vec<T2> get_coeffs() const { return coeffs; }
  void set_coeffs(const Vec<T2> &b);
  void clear() { mem.clear(); }
  Vec<T3> get_state() const;
  void set_state(const Vec<T3> &state);

private:
  virtual T3 filter(const T1 Sample);

  Vec<T3> mem;
  Vec<T2> coeffs;
  int inptr;
  bool init;
};

template <class T1, class T2, class T3>
class AR_Filter : public Filter<T1, T2, T3>
{
public:
  explicit AR_Filter();
  explicit AR_Filter(const Vec<T2> &a);
  virtual ~AR_Filter() { }
  Vec<T2> get_coeffs() const { return coeffs; }
  void set_coeffs(const Vec<T2> &a);
  void clear() { mem.clear(); }
  Vec<T3> get_state() const;
  void set_state(const Vec<T3> &state);

private:
  virtual T3 filter(const T1 Sample);

  Vec<T3> mem;
  Vec<T2> coeffs;
  Vec<T2> a0;
  int inptr;
  bool init;
};


template <class T1, class T2, class T3>
class ARMA_Filter : public Filter<T1, T2, T3>
{
public:
  explicit ARMA_Filter();
  explicit ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a);
  virtual ~ARMA_Filter() { }
  Vec<T2> get_coeffs_a() const { return acoeffs; }
  Vec<T2> get_coeffs_b() const { return bcoeffs; }
  void get_coeffs(Vec<T2> &b, Vec<T2> &a) const { b = bcoeffs; a = acoeffs; }
  void set_coeffs(const Vec<T2> &b, const Vec<T2> &a);
  void clear() { mem.clear(); }
  Vec<T3> get_state() const;
  void set_state(const Vec<T3> &state);

private:
  virtual T3 filter(const T1 Sample);

  Vec<T3> mem;
  Vec<T2> acoeffs, bcoeffs;
  int inptr;
  bool init;
};



ITPP_EXPORT vec filter(const vec &b, const vec &a, const vec &input);
ITPP_EXPORT cvec filter(const vec &b, const vec &a, const cvec &input);
ITPP_EXPORT cvec filter(const cvec &b, const cvec &a, const cvec &input);
ITPP_EXPORT cvec filter(const cvec &b, const cvec &a, const vec &input);

ITPP_EXPORT vec filter(const vec &b, const int one, const vec &input);
ITPP_EXPORT cvec filter(const vec &b, const int one, const cvec &input);
ITPP_EXPORT cvec filter(const cvec &b, const int one, const cvec &input);
ITPP_EXPORT cvec filter(const cvec &b, const int one, const vec &input);

ITPP_EXPORT vec filter(const int one, const vec &a, const vec &input);
ITPP_EXPORT cvec filter(const int one, const vec &a, const cvec &input);
ITPP_EXPORT cvec filter(const int one, const cvec &a, const cvec &input);
ITPP_EXPORT cvec filter(const int one, const cvec &a, const vec &input);


ITPP_EXPORT vec filter(const vec &b, const vec &a, const vec &input, const vec &state_in, vec &state_out);
ITPP_EXPORT cvec filter(const vec &b, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out);
ITPP_EXPORT cvec filter(const cvec &b, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out);
ITPP_EXPORT cvec filter(const cvec &b, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out);

ITPP_EXPORT vec filter(const vec &b, const int one, const vec &input, const vec &state_in, vec &state_out);
ITPP_EXPORT cvec filter(const vec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out);
ITPP_EXPORT cvec filter(const cvec &b, const int one, const cvec &input, const cvec &state_in, cvec &state_out);
ITPP_EXPORT cvec filter(const cvec &b, const int one, const vec &input, const cvec &state_in, cvec &state_out);

ITPP_EXPORT vec filter(const int one, const vec &a, const vec &input, const vec &state_in, vec &state_out);
ITPP_EXPORT cvec filter(const int one, const vec &a, const cvec &input, const cvec &state_in, cvec &state_out);
ITPP_EXPORT cvec filter(const int one, const cvec &a, const cvec &input, const cvec &state_in, cvec &state_out);
ITPP_EXPORT cvec filter(const int one, const cvec &a, const vec &input, const cvec &state_in, cvec &state_out);
ITPP_EXPORT vec fir1(int N, double cutoff);

//----------------------------------------------------------------------------
// Implementation of templated functions starts here
//----------------------------------------------------------------------------

//---------------------- class Filter ----------------------------

template <class T1, class T2, class T3>
Vec<T3> Filter<T1, T2, T3>::operator()(const Vec<T1> &x)
{
  Vec<T3> y(x.length());

  for (int i = 0; i < x.length(); i++) {
    y[i] = filter(x[i]);
  }

  return y;
}

//-------------------------- class MA_Filter ---------------------------------

template <class T1, class T2, class T3>
MA_Filter<T1, T2, T3>::MA_Filter() : Filter<T1, T2, T3>()
{
  inptr = 0;
  init = false;
}

template <class T1, class T2, class T3>
MA_Filter<T1, T2, T3>::MA_Filter(const Vec<T2> &b) : Filter<T1, T2, T3>()
{
  set_coeffs(b);
}


template <class T1, class T2, class T3>
void MA_Filter<T1, T2, T3>::set_coeffs(const Vec<T2> &b)
{
  it_assert(b.size() > 0, "MA_Filter: size of filter is 0!");

  coeffs = b;
  mem.set_size(coeffs.size(), false);
  mem.clear();
  inptr = 0;
  init = true;
}

template <class T1, class T2, class T3>
Vec<T3> MA_Filter<T1, T2, T3>::get_state() const
{
  it_assert(init == true, "MA_Filter: filter coefficients are not set!");

  int offset = inptr;
  Vec<T3> state(mem.size());

  for (int n = 0; n < mem.size(); n++) {
    state(n) = mem(offset);
    offset = (offset + 1) % mem.size();
  }

  return state;
}

template <class T1, class T2, class T3>
void MA_Filter<T1, T2, T3>::set_state(const Vec<T3> &state)
{
  it_assert(init == true, "MA_Filter: filter coefficients are not set!");
  it_assert(state.size() == mem.size(), "MA_Filter: Invalid state vector!");

  mem = state;
  inptr = 0;
}

template <class T1, class T2, class T3>
T3 MA_Filter<T1, T2, T3>::filter(const T1 Sample)
{
  it_assert(init == true, "MA_Filter: Filter coefficients are not set!");
  T3 s = 0;

  mem(inptr) = Sample;
  int L = mem.length() - inptr;

  for (int i = 0; i < L; i++) {
    s += coeffs(i) * mem(inptr + i);
  }
  for (int i = 0; i < inptr; i++) {
    s += coeffs(L + i) * mem(i);
  }

  inptr--;
  if (inptr < 0)
    inptr += mem.length();

  return s;
}

//---------------------- class AR_Filter ----------------------------------

template <class T1, class T2, class T3>
AR_Filter<T1, T2, T3>::AR_Filter() : Filter<T1, T2, T3>()
{
  inptr = 0;
  init = false;
}

template <class T1, class T2, class T3>
AR_Filter<T1, T2, T3>::AR_Filter(const Vec<T2> &a) : Filter<T1, T2, T3>()
{
  set_coeffs(a);
}

template <class T1, class T2, class T3>
void AR_Filter<T1, T2, T3>::set_coeffs(const Vec<T2> &a)
{
  it_assert(a.size() > 0, "AR_Filter: size of filter is 0!");
  it_assert(a(0) != T2(0), "AR_Filter: a(0) cannot be 0!");

  a0.set_size(1);//needed to keep the first coefficient for future reuse
  a0(0) = a(0);
  coeffs = a / a0(0);

  mem.set_size(coeffs.size() - 1, false);
  mem.clear();
  inptr = 0;
  init = true;
}


template <class T1, class T2, class T3>
Vec<T3> AR_Filter<T1, T2, T3>::get_state() const
{
  it_assert(init == true, "AR_Filter: filter coefficients are not set!");

  int offset = inptr;
  Vec<T3> state(mem.size());

  for (int n = 0; n < mem.size(); n++) {
    state(n) = mem(offset);
    offset = (offset + 1) % mem.size();
  }

  return state;
}

template <class T1, class T2, class T3>
void AR_Filter<T1, T2, T3>::set_state(const Vec<T3> &state)
{
  it_assert(init == true, "AR_Filter: filter coefficients are not set!");
  it_assert(state.size() == mem.size(), "AR_Filter: Invalid state vector!");

  mem = state;
  inptr = 0;
}

template <class T1, class T2, class T3>
T3 AR_Filter<T1, T2, T3>::filter(const T1 Sample)
{
  it_assert(init == true, "AR_Filter: Filter coefficients are not set!");
  T3 s = Sample;

  if (mem.size() == 0)
    return (s / a0(0));

  int L = mem.size() - inptr;
  for (int i = 0; i < L; i++) {
    s -= mem(i + inptr) * coeffs(i + 1); // All coeffs except a(0)
  }
  for (int i = 0; i < inptr; i++) {
    s -= mem(i) * coeffs(L + i + 1); // All coeffs except a(0)
  }

  inptr--;
  if (inptr < 0)
    inptr += mem.size();
  mem(inptr) = s;

  return (s / a0(0));
}


//---------------------- class ARMA_Filter ----------------------------------
template <class T1, class T2, class T3>
ARMA_Filter<T1, T2, T3>::ARMA_Filter() : Filter<T1, T2, T3>()
{
  inptr = 0;
  init = false;
}

template <class T1, class T2, class T3>
ARMA_Filter<T1, T2, T3>::ARMA_Filter(const Vec<T2> &b, const Vec<T2> &a) : Filter<T1, T2, T3>()
{
  set_coeffs(b, a);
}

template <class T1, class T2, class T3>
void ARMA_Filter<T1, T2, T3>::set_coeffs(const Vec<T2> &b, const Vec<T2> &a)
{
  it_assert(a.size() > 0 && b.size() > 0, "ARMA_Filter: size of filter is 0!");
  it_assert(a(0) != T2(0), "ARMA_Filter: a(0) cannot be 0!");

  acoeffs = a / a(0);
  bcoeffs = b / a(0);

  mem.set_size(std::max(a.size(), b.size()) - 1, false);
  mem.clear();
  inptr = 0;
  init = true;
}

template <class T1, class T2, class T3>
Vec<T3> ARMA_Filter<T1, T2, T3>::get_state() const
{
  it_assert(init == true, "ARMA_Filter: filter coefficients are not set!");

  int offset = inptr;
  Vec<T3> state(mem.size());

  for (int n = 0; n < mem.size(); n++) {
    state(n) = mem(offset);
    offset = (offset + 1) % mem.size();
  }

  return state;
}

template <class T1, class T2, class T3>
void ARMA_Filter<T1, T2, T3>::set_state(const Vec<T3> &state)
{
  it_assert(init == true, "ARMA_Filter: filter coefficients are not set!");
  it_assert(state.size() == mem.size(), "ARMA_Filter: Invalid state vector!");

  mem = state;
  inptr = 0;
}

template <class T1, class T2, class T3>
T3 ARMA_Filter<T1, T2, T3>::filter(const T1 Sample)
{
  it_assert(init == true, "ARMA_Filter: Filter coefficients are not set!");
  T3 z = Sample;
  T3 s;

  for (int i = 0; i < acoeffs.size() - 1; i++) { // All AR-coeff except a(0).
    z -= mem((i + inptr) % mem.size()) * acoeffs(i + 1);
  }
  s = z * bcoeffs(0);

  for (int i = 0; i < bcoeffs.size() - 1; i++) { // All MA-coeff except b(0).
    s += mem((i + inptr) % mem.size()) * bcoeffs(i + 1);
  }

  inptr--;
  if (inptr < 0)
    inptr += mem.size();
  mem(inptr) = z;

  mem(inptr) = z; // Store in the internal state.

  return s;
}


// ----------------------------------------------------------------------
// Instantiations
// ----------------------------------------------------------------------
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT MA_Filter<double, double, double>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT MA_Filter< double, std::complex<double>,
  std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT MA_Filter< std::complex<double>, double,
  std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT MA_Filter< std::complex<double>, std::complex<double>,
  std::complex<double> >;

ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT AR_Filter<double, double, double>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT AR_Filter< double, std::complex<double>,
  std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT AR_Filter< std::complex<double>,
  double, std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT AR_Filter< std::complex<double>, std::complex<double>,
  std::complex<double> >;

ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT ARMA_Filter<double, double, double>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT ARMA_Filter< double, std::complex<double>,
  std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT ARMA_Filter< std::complex<double>,
  double, std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT ARMA_Filter< std::complex<double>, std::complex<double>,
  std::complex<double> >;


} // namespace itpp

#endif // #ifndef FILTER_H