#ifndef _Stat_H
#define _Stat_H

//! $Id: Stat.H,v 1.3 2003/04/04 14:18:33 nbecker Exp $

#include <complex>

namespace stat {

template<typename T>
struct Stat_t {
  typedef T value_t;
};

template<typename T>
struct Stat_t< std::complex<T> > {
  typedef typename std::complex<T>::value_type value_t;
};

template<typename T>
T Norm (T x) { return x * x ; }

template<typename T>
T Norm (std::complex<T> x) { return real (x) * real (x) + imag (x) * imag (x); }

template<typename coef_t, typename in_t = coef_t, typename out_t = in_t>
class one_pole_filter {
public:
  one_pole_filter (coef_t timeConstant) :
    inputGain (1./timeConstant),
    feedbackGain (1. - 1./timeConstant),
    sum (0)
  {}

  out_t Compute (in_t in) {
    sum = sum * feedbackGain + in * inputGain;
    return sum;
  }

  out_t Get() const { return sum; }

  void Reset() { sum = 0; }

private:
  const coef_t inputGain;
  const coef_t feedbackGain;
  out_t sum;
};

template<typename T, typename FILT1=one_pole_filter<typename Stat_t<T>::value_t>, typename FILT2=one_pole_filter<typename Stat_t<T>::value_t, T, T> >
class Stat {
  typedef typename Stat_t<T>::value_t value_t;

  FILT1 xSqrFilt;
  FILT2 xFilt;
public:
  Stat (value_t timeconstant=1000.) :
    xSqrFilt (timeconstant),
    xFilt (timeconstant)
    {}

  void Reset() {
    xFilt.Reset();
    xSqrFilt.Reset();
  }

  template<typename in_t>
  void Compute (in_t in, in_t inend) {
    for (; in != inend; in++)
      Compute (*in);
  }

  void Compute (T x) {
    operator += (x);
  }

  void operator += (T x) {
    xFilt.Compute (x);
    xSqrFilt.Compute (Norm (x));
  }

  void operator() (T x) {
    operator += (x);
  }

  T Mean() const {
    return (xFilt.Get());
  }

  value_t Var() const {
    return xSqrFilt.Get() - Norm (xFilt.Get());
  }

  value_t StdDev() const {
    return sqrt (Var());
  }

  value_t RMS() const {
    return StdDev();
  }

  value_t XsqrBar() const {
    return xSqrFilt.Get();
  }
  
};

}

#endif