mspass::algorithms::Butterworth Class Reference

MsPASS implementation of Butterworth filter as processing object. More...

#include <Butterworth.h>

Public Member Functions
	Butterworth ()
	Default constructor.
	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)
	Fully parameterized constructor with args similar to subfilt.
	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)
	Construct by defining corner frequencies and number of npoles.
	Butterworth (const Butterworth &parent)
Butterworth &	operator= (const Butterworth &parent)
mspass::seismic::CoreTimeSeries	impulse_response (const int n)
	Return the impulse response.
void	apply (mspass::seismic::CoreTimeSeries &d)
	pply the filter to a CoreTimeSeries object.
void	apply (mspass::seismic::TimeSeries &d)
	Apply the filter to a CoreTimeSeries object.
void	apply (std::vector< double > &d)
	Filter a raw vector of data.
void	apply (mspass::seismic::CoreSeismogram &d)
	Apply the filter to a CoreSeismogram object.
void	apply (mspass::seismic::Seismogram &d)
	Apply the filter to a CoreTimeSeries object.
mspass::algorithms::deconvolution::ComplexArray	transfer_function (const int n)
	Return the response of the filter in the frequency domain.
void	change_dt (const double dtnew)
	set the sample interval assumed for input data.
double	low_corner () const
double	high_corner () const
int	npoles_low () const
int	npoles_high () const
double	current_dt () const
std::string	filter_type () const
bool	is_zerophase () const

Detailed Description

MsPASS implementation of Butterworth filter as processing object.

MsPASS has an existing filter routine that can implement buterworth filters via obspy. This class was created to allow a clean interface to Butterworth filtering from C++ code that needs such an operator. The original use was an experimental deconvolution code, but there will likely be others because simple, efficient filters are a common internal need for potential applications.

This C++ class can be viewed as a wrapper for Seismic Unix functions that implement Butterworth filters. The parent functions were found in the cwp lib and were called bfdesign, bfhighpass, and bflowpass. Those three functions are the central tools used to implement this class. All the rest is really just a wrapper to provide an object oriented api to the su functions.

This class implements a processing object concept. That is, it is intended to be constructed and then used for processing multiple data objects with fixed parameters. For Butterwoth filtering the parameters are relatively simple (mainly two corner frequencies and number of poles defining the filter rolloff). A complexity, however, is that the class was designed to allow automatic handling of multiple sample rate data. That is handled internally by caching the sample interval of the data and automatically adjusting the coefficients when the sample interval changes. Note that feature only works for MsPASS data objects CoreTimeSeries and Seismogram where the sample interval is embedded in the object. The raw interface with a simple vector cannot know that. The method to change the expected sample interval has some sanity checks to reduce, but not eliminate the possibility of mistakes that will create unstable filters.

Constructor & Destructor Documentation

Butterworth() [1/5]

mspass::algorithms::Butterworth::Butterworth

(

)

Default constructor.

The default constructor does not define a null. The default generates an antialiasing filter identical to the default in the antialias function in seismic unix. That is it produces a low pass filter with a band edge (pass parameter) at 60% of Nyquist and a stop edge at Nyquist. It then calls the su bfdesign function to compute the number of poles and the 3db frequency of the corner to define this filter. These can be retrieved with getters (see below)

                         {
  /* This is a translation of a minimum phase antialias filter used in
  the antialias function in seismic unix.  We define it as a default
  because it does almost nothing - almost because it depends upon
  form of antialias applied to get to current state.  */
  double fnyq, fpass, apass, fstop, astop;
  fnyq = 0.5;
  fpass = 0.6 * fnyq;
  apass = 0.99;
  fstop = fnyq;
  astop = 0.01;
  this->bfdesign(fpass, apass, fstop, astop, &(this->npoles_hi),
                 &(this->f3db_hi));
  dt = 1.0;
  use_lo = false;
  use_hi = true;
  f3db_lo = 0.0;
  npoles_lo = 0;
  zerophase = false;
};

Butterworth() [2/5]

mspass::algorithms::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
	)

Fully parameterized constructor with args similar to subfilt.

A butterworth filter can be described two ways: (1) corner frequency and number of poles and (2) by band stop and band frequencies. This constructor is used to define the filter by stop and pass band parameters. Frequencies must satisfy fstoplo<fpasslo and fpasshi<fstophi as the four frequencis define a bandd pass between the fpasslo and fpasshi. The stop frequencies define where the response should near zero. Thus for band pass filters the apasslo and apasshi should be 1 and the stop Amplitudes a small number like 0.01. For a band reject filter set stop amplitudes to 1 and pass amplitudes to small numbers like 0.01. (For a reject filter the pass frequencies act like stop frequencies for a bandbpass filer - this is mostly like the subfilt seismic unix program). The booleans control which terms are enabld. When enable_lo is true the lo components are used an when enable_hi is true the high componens are used.

There is a confusing nomenclature related to "high" and "low". In this implmentation I always take low to mean the low side of the passband and high to be the high side of the passband as described above for the 4 frequency point defiitions. The issue is "lowcut" versus "lowpass". Seismic Unix really mixes this up as their implmenetation (which I used here) refernces bflowpass and bfhighpass but subfilt uses the inverse lowcut and higcut terminology. A geeky implementation detail is I actually changed the names of the functions to eliminate the confusion in the implementation. That matters only if you want to compare what we did here to the original seismic unix code.

Note the iir filter coefficiets are always derived from the poles and Frequencies so this constructor is just an alternate way to define the filter without the abstraction of number of poles.

Parameters

zerophase	when true use a zerophase filter. When false defines a one pass minimum phase filter.
enable_lo	is a boolean that when true enables the low band parameters for the filter (i.e. the highpass=low-cut components)
enable_hi	is a boolean taht when true enables the parameters defining the upper frequency band edge (i.e. lowpass=high-cut parameters)
fstoplo	- stop band frequency for lower band edge
astoplo	- amplitude at stop frequency (small number for band pass, 1 for band reject)
fpasslo	- pass band frequency for lower band edge
apasslo	- amplitude at fpasslo frequency (1 for bandpass, small number for band reject)
fstophi	- stop band frequency for upper band edge
astophi	- amplitude at stop frequency (small number for band pass, 1 for band reject)
fpasshi	- pass band frequency for upper band edge
apasshi	- amplitude at fpasshi frequency (1 for bandpass, small number for band reject)
sample_interval	is the expected data sample interval

                                                       {
  dt = sample_interval;
  zerophase = zp;
  use_lo = lcon;
  use_hi = hcon;
  this->bfdesign(fpasslo * dt, apasslo, fstoplo * dt, astoplo,
                 &(this->npoles_lo), &(this->f3db_lo));
  this->bfdesign(fpasshi * dt, apasshi, fstophi * dt, astophi,
                 &(this->npoles_hi), &(this->f3db_hi));
};

Butterworth() [3/5]

mspass::algorithms::Butterworth::Butterworth

(

const mspass::utility::Metadata &

md

)

Construct using tagged valus created from a Metadata container.

This behaves exactly like the fully parameterized contructor except it gets the parameters from metadata. Metadata keys in initial implementation are identical to the argument names defined above. The best guidance for using this constuctor is to look a the comments in the default parameter file.

Construct using tagged valus created from a Metadata container

                                           {
  string base_error("Butterworth pf constructor:  ");
  try {
    dt = md.get<double>("sample_interval");
    /* We define these here even if they aren't all used for all options.
    Better for tiny memory cost than declaration in each block */
    double fpass, apass, fstop, astop;
    zerophase = md.get_bool("zerophase");
    /* These options simplify setup but complicate this constructor.  */
    string ftype, fdmeth;
    ftype = md.get_string("filter_type");
    fdmeth = md.get_string("filter_definition_method");
    if (ftype == "bandpass") {
      use_lo = true;
      use_hi = true;
      if (fdmeth == "corner_pole") {
        /* Note we allow flow>figh to create strongly peaked filters although
        we don't scale these any any rational way */
        npoles_lo = md.get<int>("npoles_low");
        npoles_hi = md.get<int>("npoles_high");
        f3db_lo = md.get<double>("corner_low");
        f3db_hi = md.get<double>("corner_high");
        /* Stored as nondimensional internally*/
        f3db_lo *= dt;
        f3db_hi *= dt;
      } else {
        /* note the set routines will test validity of input and throw
        an error if the points do not define the right sense of pass versus
        cut*/
        fpass = md.get<double>("fpass_low");
        apass = md.get<double>("apass_low");
        fstop = md.get<double>("fstop_low");
        astop = md.get<double>("astop_low");
        this->set_lo(fstop, fpass, astop, apass);
        fpass = md.get<double>("fpass_high");
        apass = md.get<double>("apass_high");
        fstop = md.get<double>("fstop_high");
        astop = md.get<double>("astop_high");
        this->set_hi(fstop, fpass, astop, apass);
      }
    } else if (ftype == "lowpass") {
      /* Note the confusing naming hre - lo means lo corner not low pass*/
      use_lo = false;
      use_hi = true;
      /* Explicity initialization useful here */
      npoles_lo = 0;
      f3db_lo = 0.0;
      if (fdmeth == "corner_pole") {
        npoles_hi = md.get<int>("npoles_high");
        f3db_hi = md.get<double>("corner_high");
        f3db_hi *= dt;
      } else {
        fpass = md.get<double>("fpass_high");
        apass = md.get<double>("apass_high");
        fstop = md.get<double>("fstop_high");
        astop = md.get<double>("astop_high");
        /* note this handles inconsistencies and can throw an error */
        this->set_hi(fstop, fpass, astop, apass);
      }
    } else if (ftype == "highpass") {
      use_lo = true;
      use_hi = false;
      /* Explicity initialization useful here */
      npoles_hi = 0;
      f3db_hi = 0.0;
      if (fdmeth == "corner_pole") {
        npoles_lo = md.get<int>("npoles_low");
        f3db_lo = md.get<double>("corner_low");
        f3db_lo *= dt;
      } else {
        fpass = md.get<double>("fpass_low");
        apass = md.get<double>("apass_low");
        fstop = md.get<double>("fstop_low");
        astop = md.get<double>("astop_low");
        /* note this handles inconsistencies and can throw an error */
        this->set_lo(fstop, fpass, astop, apass);
      }
    } else if (ftype == "bandreject") {
      use_lo = true;
      use_hi = true;
      if (fdmeth == "corner_pole") {
        /* Note we allow flow>figh to create strongly peaked filters although
        we don't scale these any any rational way */
        npoles_lo = md.get<int>("npoles_low");
        npoles_hi = md.get<int>("npoles_high");
        f3db_lo = md.get<double>("corner_low");
        f3db_hi = md.get<double>("corner_high");
        /* Enforce this constraint or we don't have a band reject filtr */
        if (f3db_lo < f3db_hi) {
          throw MsPASSError(base_error +
                                "Illegal corner frequencies for band reject "
                                "filter definition" +
                                "For a band reject low corner must be larger "
                                "than high corner (pass becomes cut to reject)",
                            ErrorSeverity::Invalid);
        }
        f3db_lo *= dt;
        f3db_hi *= dt;
      } else {
        /* We should test these values for defining band reject but bypass
        these safeties for now as I don't expect this feature to get much
        use. */
        fpass = md.get<double>("fpass_low");
        apass = md.get<double>("apass_low");
        fstop = md.get<double>("fstop_low");
        astop = md.get<double>("astop_low");
        this->bfdesign(fpass * dt, apass, fstop * dt, astop, &(this->npoles_lo),
                       &(this->f3db_lo));
        fpass = md.get<double>("fpass_high");
        apass = md.get<double>("apass_high");
        fstop = md.get<double>("fstop_high");
        astop = md.get<double>("astop_high");
        this->bfdesign(fpass * dt, apass, fstop * dt, astop, &(this->npoles_hi),
                       &(this->f3db_hi));
      }
    } else {
      throw MsPASSError(base_error + "Unsupported filtr_type=" + ftype,
                        ErrorSeverity::Invalid);
    }
  } catch (...) {
    throw;
  };
}

References mspass::utility::Metadata::get(), mspass::utility::Metadata::get_bool(), and mspass::utility::Metadata::get_string().

Butterworth() [4/5]

mspass::algorithms::Butterworth::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
	)

Construct by defining corner frequencies and number of npoles.

Butterworth filters can also be defind by a corner frequency and number of poles. In fact, only the nondimensional form of these parameters are stored as private attributes to define the filter.

Parameters

zerophase	when true use a zerophase filter. When false defines a one pass minimum phase filter.
enable_lo	is a boolean that when true enables the low band parameters for the filter (i.e. the highpass=low-cut components)
enable_hi	is a boolean taht when true enables the parameters defining the upper frequency band edge (i.e. lowpass=high-cut parameters)
npolelo	is the number of poles for the low frequency corner (highpass)
f3dblo	is the corner frequency for the low frequency corner (highpass)
npolehi	is the number of poles for the high frequency corner (lowpass)
f3dbhi	is the corner frequency for the high frequency corner (lowpass)
sample_interval	is the expected data sample interval

                                                       {
  this->dt = sample_interval;
  this->zerophase = zero;
  this->use_lo = locut;
  this->use_hi = hicut;
  this->f3db_hi = f3dbhi;
  this->f3db_lo = f3dblo;
  this->npoles_lo = npolelo;
  this->npoles_hi = npolehi;
  /* Convert the frequencies to nondimensional form */
  this->f3db_lo *= dt;
  this->f3db_hi *= dt;
}

Butterworth() [5/5]

mspass::algorithms::Butterworth::Butterworth

(

const Butterworth &

parent

)

Standard copy conststructor.

                                                  {
  use_lo = parent.use_lo;
  use_hi = parent.use_hi;
  zerophase = parent.zerophase;
  f3db_lo = parent.f3db_lo;
  f3db_hi = parent.f3db_hi;
  npoles_lo = parent.npoles_lo;
  npoles_hi = parent.npoles_hi;
  dt = parent.dt;
}

Member Function Documentation

apply() [1/5]

void mspass::algorithms::Butterworth::apply

(

mspass::seismic::CoreSeismogram &

d

)

Apply the filter to a CoreSeismogram object.

This method alters the data vector inside d in place and changes no other parts of the data. Automatic switching of data sample rate is used on the operator. That is, if the sample rate of the data is different than the operator sample rate the internal operator coefficients will be adjusted to the new sample rate. The operator sample rate will also be changed to the sample rate of d whenever the sample rate changes from the previous call.

This method has a safety to prevent irrational sample rate changes. The IRR filter used to compute a Butterworth filter becomes unstable if the low pass filter component (high corner) approach Nyquist or worse exceed Nyquist. This method will throw a MsPASSError exception if the sample rate of d is too low for the filter high corner. (current 90% of Nyquist). When this error is throw the data will be unaltered and the internal sample rate will be left in the previous state.

Parameters

d	input data to be filtered - altered in place.

Exceptions

throws

a MsPASSError if the hi corner is inconsistent with the sample rate of d

                                         {
  /* we could just handle a MsPASSError to take care of exceptions
  thrown by apply if the sample rate is illegal, but this is more adaptable. */
  try {
    if (d.npts() == 0)
      return;
    double d_dt = d.dt();
    if ((this->dt != d_dt) && use_hi &&
        (this->f3db_hi * d_dt / this->dt > FHighFloor)) {
      stringstream ss;
      ss << "Butterworth::apply:  automatic dt change error" << endl
         << "Current operator dt=" << this->dt << " data dt=" << d_dt << endl
         << "Change would produce a corner too close to Nyquist"
         << " and create an unstable filter" << endl
         << "Use a different filter operator for data with this sample rate"
         << endl;
      throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
    }
    if (this->dt != d_dt)
      this->change_dt(d_dt);
    /* We copy each component to this buffer, filter, and then copy
    back.  Not as efficient as if we used a strike parameter in the filter
    functions, but I did not want to rewrite those functions. */
    vector<double> comp;
    int npts = d.npts();
    comp.reserve(npts);
    /* This initializes the buffer to allow us a simpler loop from 0 to 2
    using the blas function dcopy to handle the skip for a fortran style
    array used for storing the 3c data in  u*/
    for (auto i = 0; i < npts; ++i)
      comp.push_back(0.0);
    for (auto k = 0; k < 3; ++k) {
      dcopy(npts, d.u.get_address(k, 0), 3, &(comp[0]), 1);
      this->apply(comp);
      dcopy(npts, &(comp[0]), 1, d.u.get_address(k, 0), 3);
    }
  } catch (...) {
    throw;
  };
}

References apply(), change_dt(), mspass::seismic::BasicTimeSeries::dt(), mspass::utility::dmatrix::get_address(), mspass::seismic::BasicTimeSeries::npts(), and mspass::seismic::CoreSeismogram::u.

apply() [2/5]

void mspass::algorithms::Butterworth::apply

(

mspass::seismic::CoreTimeSeries &

d

)

pply the filter to a CoreTimeSeries object.

This method alters the data vector inside d in place and changes no other parts of the data. Automatic switching of data sample rate is used on the operator. That is, if the sample rate of the data is different than the operator sample rate the internal operator coefficients will be adjusted to the new sample rate. The operator sample rate will also be changed to the sample rate of d whenever the sample rate changes from the previous call.

This method has a safety to prevent irrational sample rate changes. The IRR filter used to compute a Butterworth filter becomes unstable if the low pass filter component (high corner) approach Nyquist or worse exceed Nyquist. This method will throw a MsPASSError exception if the sample rate of d is too low for the filter high corner. (current 90% of Nyquist). When this error is throw the data will be unaltered and the internal sample rate will be left in the previous state.

Parameters

d	input data to be filtered - altered in place.

Exceptions

throws

a MsPASSError if the hi corner is inconsistent with the sample rate of d

                                                      {
  if (d.npts() == 0)
    return;
  double d_dt = d.dt();
  if (this->dt != d_dt) {
    /* Here we throw an exception if a requested sample rate is
    illegal */
    double fhtest = this->f3db_hi * d_dt / this->dt;
    if (use_hi && (fhtest > FHighFloor)) {
      stringstream ss;
      ss << "Butterworth::apply:  automatic dt change error" << endl
         << "Current operator dt=" << this->dt << " data dt=" << d_dt << endl
         << "Change would produce a corner too close to Nyquist"
         << " and create an unstable filter" << endl
         << "Use a different filter operator for data with this sample rate"
         << endl;
      throw MsPASSError(ss.str(), ErrorSeverity::Invalid);
    }
    this->change_dt(d_dt);
  }
  this->apply(d.s);
}

References apply(), change_dt(), mspass::seismic::BasicTimeSeries::dt(), mspass::seismic::BasicTimeSeries::npts(), and mspass::seismic::CoreTimeSeries::s.

apply() [3/5]

void mspass::algorithms::Butterworth::apply

(

mspass::seismic::Seismogram &

d

)

Apply the filter to a CoreTimeSeries object.

This method alters the data vector inside d in place and changes no other parts of the data. Automatic switching of data sample rate is used on the operator. That is, if the sample rate of the data is different than the operator sample rate the internal operator coefficients will be adjusted to the new sample rate. The operator sample rate will also be changed to the sample rate of d whenever the sample rate changes from the previous call.

This method has a safety to prevent irrational sample rate changes. The IRR filter used to compute a Butterworth filter becomes unstable if the low pass filter component (high corner) approach Nyquist or worse exceed Nyquist. This method will automatically disable the high corner (lowpass) component of the filter if the corner approaches or exceed Nyquist. When that happens the internal sample rate is restored to the previous value and a complaint message is posted to elog of d.

Parameters

d	input data to be filtered - altered in place.

Exceptions

none,but

callers should consider checking for errors posted to elog

                                                  {
  if (d.npts() == 0)
    return;
  double d_dt = d.dt();
  if (this->dt != d_dt) {
    /* Here we throw an exception if a requested sample rate is
    illegal */
    double fhtest = this->f3db_hi * d_dt / this->dt;
    if (use_hi && (fhtest > FHighFloor)) {
      Butterworth adjusted(*this);
      adjusted.use_hi = false;
      adjusted.change_dt(d_dt);
      adjusted.apply(d);
      stringstream ss;
      ss << "Auto adjust for sample rate change error" << endl
         << "Upper corner of filter=" << this->high_corner()
         << " is near or above Nyquist frequency for requested sample "
         << "interval=" << d_dt << endl
         << "Disabling upper corner (lowpass) and applying filter anyway"
         << endl;
      d.elog.log_error(string("Butterworth::apply"), ss.str(),
                       ErrorSeverity::Complaint);
      return;
    } else {
      this->change_dt(d_dt);
    }
  }
  vector<double> comp;
  int npts = d.npts();
  comp.reserve(npts);
  /* This initializes the buffer to allow us a simpler loop from 0 to 2
  using the blas function dcopy to handle the skip for a fortran style
  array used for storing the 3c data in  u*/
  for (auto i = 0; i < npts; ++i)
    comp.push_back(0.0);
  for (auto k = 0; k < 3; ++k) {
    dcopy(npts, d.u.get_address(k, 0), 3, &(comp[0]), 1);
    this->apply(comp);
    dcopy(npts, &(comp[0]), 1, d.u.get_address(k, 0), 3);
  }
}

References apply(), change_dt(), mspass::seismic::BasicTimeSeries::dt(), mspass::utility::ProcessingHistory::elog, mspass::utility::dmatrix::get_address(), high_corner(), mspass::utility::ErrorLogger::log_error(), mspass::seismic::BasicTimeSeries::npts(), and mspass::seismic::CoreSeismogram::u.

apply() [4/5]

void mspass::algorithms::Butterworth::apply

(

mspass::seismic::TimeSeries &

d

)

Apply the filter to a CoreTimeSeries object.

This method alters the data vector inside d in place and changes no other parts of the data. Automatic switching of data sample rate is used on the operator. That is, if the sample rate of the data is different than the operator sample rate the internal operator coefficients will be adjusted to the new sample rate. The operator sample rate will also be changed to the sample rate of d whenever the sample rate changes from the previous call.

This method has a safety to prevent irrational sample rate changes. The IRR filter used to compute a Butterworth filter becomes unstable if the low pass filter component (high corner) approach Nyquist or worse exceed Nyquist. This method will automatically disable the high corner (lowpass) component of the filter if the corner approaches or exceed Nyquist. When that happens the internal sample rate is restored to the previous value and a complaint message is posted to elog of d.

Parameters

d	input data to be filtered - altered in place.

Exceptions

none,but

callers should consider checking for errors posted to elog

                                                  {
  if (d.npts() == 0)
    return;
  double d_dt = d.dt();
  if (this->dt != d_dt) {
    /* Here we throw an exception if a requested sample rate is
    illegal */
    double fhtest = this->f3db_hi * d_dt / this->dt;
    if (use_hi && (fhtest > FHighFloor)) {
      Butterworth adjusted(*this);
      adjusted.use_hi = false;
      adjusted.change_dt(d_dt);
      adjusted.apply(d.s);
      stringstream ss;
      ss << "Auto adjust for sample rate change error" << endl
         << "Upper corner of filter=" << this->high_corner()
         << " is near or above Nyquist frequency for requested sample "
         << "interval=" << d_dt << endl
         << "Disabling upper corner (lowpass) and applying filter anyway"
         << endl;
      d.elog.log_error(string("Butterworth::apply"), ss.str(),
                       ErrorSeverity::Complaint);
      return;
    } else {
      this->change_dt(d_dt);
    }
  }
  this->apply(d.s);
}

References apply(), change_dt(), mspass::seismic::BasicTimeSeries::dt(), mspass::utility::ProcessingHistory::elog, high_corner(), mspass::utility::ErrorLogger::log_error(), mspass::seismic::BasicTimeSeries::npts(), and mspass::seismic::CoreTimeSeries::s.

apply() [5/5]

void mspass::algorithms::Butterworth::apply

(

std::vector< double > &

d

)

Filter a raw vector of data.

Use this method to apply the filter to a raw vector of data. The C++ interface uses an std::vector container, but the python api in MsPASS allows this to be a double numpy array or any iterable version of a vector container (meaning storage as a contiguous block of memory). If this method is used it is assumed the sample interval defined for the operator is the same as the for the input data.

Parameters

d	is the data to be filtered (note the data are altered in place)

                                            {
  if (d.empty())
    return;
  if (use_lo) {
    this->bflowcut(npoles_lo, f3db_lo, d.size(), &(d[0]), &(d[0]));
    if (zerophase) {
      reverse_vector(d.size(), &(d[0]));
      this->bflowcut(npoles_lo, f3db_lo, d.size(), &(d[0]), &(d[0]));
      reverse_vector(d.size(), &(d[0]));
    }
  }
  if (use_hi) {
    this->bfhighcut(npoles_hi, f3db_hi, d.size(), &(d[0]), &(d[0]));
    if (zerophase) {
      reverse_vector(d.size(), &(d[0]));
      this->bfhighcut(npoles_hi, f3db_hi, d.size(), &(d[0]), &(d[0]));
      reverse_vector(d.size(), &(d[0]));
    }
  }
}

change_dt()

void mspass::algorithms::Butterworth::change_dt ( const double  dtnew )

inline

set the sample interval assumed for input data.

This function can be used when running with raw data vectors if the sample interval of the data series is different from that called on construction or set previously. This is a nontrivial change because the filter coefficients depend upon sample interval. In particular, for this implementation npoles and the 3db frequency points stored internally are altered when this function is called. If the frequency intervals change the expectation is the user will create a new instance of this object.

Warning: this routine does not implement the safeties built into TimeSeries and Seismogram apply methods. It will silently change the upper corner to an unstable position if called inappropriately.

Parameters

dtnew

is the new sample interval to set for the operator.

                                     {
    this->f3db_lo *= (dtnew / (this->dt));
    this->f3db_hi *= (dtnew / (this->dt));
    this->dt = dtnew;
  };

current_dt()

double mspass::algorithms::Butterworth::current_dt ( ) const

inline

Return the current operator sample interval.

{ return dt; };

filter_type()

std::string mspass::algorithms::Butterworth::filter_type ( ) const

inline

Return a string defining the type of operator this filter defines. Currently can be one of the following: bandpass, lowpass, or highpass. It is possible to construct a band reject filter with the right constructor, but the implementation of this method will not detect that condition. A band reject filter will be incorrectly tagged bandpass. The algorithm just looks to see which of the band edges are defined (the hi and lo concepts described above) and guesses the filter type. If both are off it returns "Undefined".

                                {
    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");
    }
  };

high_corner()

double mspass::algorithms::Butterworth::high_corner ( ) const

inline

Return the high frequency 3db corner (in Hz).

{ return f3db_hi / dt; };

impulse_response()

CoreTimeSeries mspass::algorithms::Butterworth::impulse_response

(

const int

n

)

Return the impulse response.

The response of a linear filter like the butterworth filter can always be described by either the time domain impulse response or its fourier transform commonly called the tranfer function. This function returns the impulse response centered in a time window with a specified number of samples using the current sample interval cached in the object. Note the return has dt and the impulse is at the center of the data window (n/2) with t0 set so the functions zero is correct if using the implict time scale (time method) of a time series object.

Parameters

n	is the number of samples to generate to characterize the impulse response. The function is always returned centered on the vector of length n and t0 of the TimeSeries is set to make that impulse point be time 0.

                                                        {
  CoreTimeSeries result(n);
  /* We use a feature that the above constructor initiallizes the buffer
  to zeros so just a spike at n/2. */
  result.s[n / 2] = 1.0;
  result.set_t0(((double)(-n / 2)) * (this->dt));
  result.set_dt(this->dt);
  result.set_tref(TimeReferenceType::Relative);
  result.set_live();
  this->apply(result.s);
  result.s = normalize<double>(result.s);
  return result;
}

References apply(), mspass::seismic::CoreTimeSeries::s, mspass::seismic::CoreTimeSeries::set_dt(), mspass::seismic::BasicTimeSeries::set_live(), mspass::seismic::CoreTimeSeries::set_t0(), and mspass::seismic::BasicTimeSeries::set_tref().

is_zerophase()

bool mspass::algorithms::Butterworth::is_zerophase ( ) const

inline

Return true if the filter is defined as a zero phase filter. Returns false if it is minimum phase.

{ return zerophase; };

low_corner()

double mspass::algorithms::Butterworth::low_corner ( ) const

inline

Return the low frequency 3db corner (in Hz).

{ return f3db_lo / dt; };

npoles_high()

int mspass::algorithms::Butterworth::npoles_high ( ) const

inline

Return the number of poles defining the lowpass (highcut) element of the filter.

{ return npoles_hi; };

npoles_low()

int mspass::algorithms::Butterworth::npoles_low ( ) const

inline

Return the number of poles defining the highpass (lowcut) element of the filter.

{ return npoles_lo; };

operator=()

Butterworth & mspass::algorithms::Butterworth::operator=

(

const Butterworth &

parent

)

Standard assignment operator.

                                                             {
  if (this != &parent) {
    use_lo = parent.use_lo;
    use_hi = parent.use_hi;
    zerophase = parent.zerophase;
    f3db_lo = parent.f3db_lo;
    f3db_hi = parent.f3db_hi;
    npoles_lo = parent.npoles_lo;
    npoles_hi = parent.npoles_hi;
    dt = parent.dt;
  }
  return *this;
}

transfer_function()

ComplexArray mspass::algorithms::Butterworth::transfer_function

(

const int

n

)

Return the response of the filter in the frequency domain.

The impulse response of any linear system can always be characterized by either the time domain response to spike signal or the alternative frequency domain version of the same function commonly called the transfer function. This method returns the transfer funtion as a mspass::algorithms::deconvolution::ComplexArray container. Use methods in that object to get amplitude and phase response functions.

Parameters

n

is the number of points that should be used to characterize the transfer function. Note because we are dealing with strictly real valued signals the array returned will be folded at the Nyquist frequency in the standard way of all FFT implementations (current implementation uses the fft in the gnu scientific library that definitely does that).

                                                          {
  CoreTimeSeries imp = this->impulse_response(nfft);
  /* the impulse response function uses a time shift and sets t0 to the
  required position.  A disconnect with any fft is that will introduce an
  undesirable phase shift for many algorithms.  We remove it here using our
  circular shift function */
  int ishift = imp.sample_number(0.0);
  imp.s = circular_shift(imp.s, ishift);
  gsl_fft_complex_wavetable *wavetable = gsl_fft_complex_wavetable_alloc(nfft);
  gsl_fft_complex_workspace *workspace = gsl_fft_complex_workspace_alloc(nfft);
  ComplexArray work(nfft, imp.s);
  gsl_fft_complex_forward(work.ptr(), 1, nfft, wavetable, workspace);
  gsl_fft_complex_wavetable_free(wavetable);
  gsl_fft_complex_workspace_free(workspace);
  return work;
}

References impulse_response(), mspass::algorithms::deconvolution::ComplexArray::ptr(), mspass::seismic::CoreTimeSeries::s, and mspass::seismic::BasicTimeSeries::sample_number().

---

The documentation for this class was generated from the following files:

• /home/runner/work/mspass/mspass/cxx/include/mspass/algorithms/Butterworth.h
• /home/runner/work/mspass/mspass/cxx/src/lib/algorithms/Butterworth.cc