Butterworth.h

#ifndef _MSPASS_BUTTERWORTH_H_
#define _MSPASS_BUTTERWORTH_H_
#include "mspass/algorithms/deconvolution/ComplexArray.h"
#include "mspass/seismic/Seismogram.h"
#include "mspass/seismic/TimeSeries.h"
namespace mspass::algorithms {
class Butterworth {
public:
  Butterworth();
  Butterworth(const bool zerophase, const bool enable_lo, const bool enable_hi,
              const double fstoplo, const double astoplo, const double fpasslo,
              const double apasslo, const double fpasshi, const double apasshi,
              const double fstophi, const double astophi,
              const double sample_interval);
  Butterworth(const mspass::utility::Metadata &md);
  Butterworth(const bool zerophase, const bool enable_lo, const bool enable_hi,
              const int npolelo, const double f3dblo, const int npolehi,
              const double f3dbhi, const double sample_interval);
  Butterworth(const Butterworth &parent);
  Butterworth &operator=(const Butterworth &parent);
  mspass::seismic::CoreTimeSeries impulse_response(const int n);
  void apply(mspass::seismic::CoreTimeSeries &d);
  void apply(mspass::seismic::TimeSeries &d);
  void apply(std::vector<double> &d);
  void apply(mspass::seismic::CoreSeismogram &d);
  void apply(mspass::seismic::Seismogram &d);
  mspass::algorithms::deconvolution::ComplexArray
  transfer_function(const int n);
  void change_dt(const double dtnew) {
    this->f3db_lo *= (dtnew / (this->dt));
    this->f3db_hi *= (dtnew / (this->dt));
    this->dt = dtnew;
  };
  double low_corner() const { return f3db_lo / dt; };
  double high_corner() const { return f3db_hi / dt; };
  int npoles_low() const { return npoles_lo; };
  int npoles_high() const { return npoles_hi; };
  double current_dt() const { return dt; };
  std::string filter_type() const {
    if (use_lo) {
      if (use_hi)
        return std::string("bandpass");
      else
        return std::string("highpass");
    } else {
      if (use_hi)
        return std::string("lowpass");
      else
        return std::string("Undefined");
    }
  };
  bool is_zerophase() const { return zerophase; };
 
private:
  bool use_lo, use_hi;
  bool zerophase;
 
  /* bfdesign sets npoles and 3db points based on specified band properties.
  These stored these when the constructor calls bfdesign for high and
  low pass elements. Names should make clear which is which. These are
  always stored in nondimensional form (f_nd=f*dt)*/
  double f3db_lo, f3db_hi;
  int npoles_lo, npoles_hi;
  double dt;
  /* These three functions are nearly exact copies of seismic unix
  functions with a couple of differences:
  1.  All floats in the su code are made double here to mesh with modern arch
  2.  bfhighpass is changed to bflowpass and bflowpass is changed to bfhighcut
  to mesh with the names used here more cleanly.  that is, pass and cut define
  the opposite sense of a filter and it gets very confusing when the
  terms get mixed up in the same set of code.  This way low always means the
  low end of the pass band and high is the high end of the pass band.
  */
 
  void bfdesign(double fpass, double apass, double fstop, double astop,
                int *npoles, double *f3db);
  void bfhighcut(int npoles, double f3db, int n, double p[], double q[]);
  void bflowcut(int npoles, double f3db, int n, double p[], double q[]);
  /* These internal methods use internal dt value to call bfdesign with
  nondimensional frequencies. Hence they should always be called after
  a change in dt.  They are also handy to contain common code for constructors.
*/
  void set_lo(const double fstop, const double fpass, const double astop,
              const double apass);
  void set_hi(const double fstop, const double fpass, const double astop,
              const double apass);
};
} // namespace mspass::algorithms
#endif