[seiscomp, scanloc] Install, add .gitignore
This commit is contained in:
157
include/seiscomp/math/conversions.h
Normal file
157
include/seiscomp/math/conversions.h
Normal file
@ -0,0 +1,157 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_MT_CONVERSIONS_H
|
||||
#define SEISCOMP_MATH_MT_CONVERSIONS_H
|
||||
|
||||
|
||||
#include <seiscomp/math/vector3.h>
|
||||
#include <seiscomp/math/tensor.h>
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
struct NODAL_PLANE {
|
||||
double str;
|
||||
double dip;
|
||||
double rake;
|
||||
};
|
||||
|
||||
|
||||
struct AXIS {
|
||||
double str;
|
||||
double dip;
|
||||
double val;
|
||||
int expo;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* Calculates the strike and dip of an axis from a cartesian vector
|
||||
*/
|
||||
template <typename T>
|
||||
void vec2angle(const Math::Vector3<T> &v, T &strike, T &dip);
|
||||
|
||||
|
||||
/**
|
||||
* Calculates the axis of a cartesian vector from strike and dip.
|
||||
*/
|
||||
template <typename T>
|
||||
void angle2vec(T strike, T dip, Math::Vector3<T> &v);
|
||||
|
||||
|
||||
/**
|
||||
* Calculates the fault normal and slip vector from principle axis P and T
|
||||
*/
|
||||
template <typename T>
|
||||
bool pa2nd(const Math::Vector3<T> &t, const Math::Vector3<T> &p, Math::Vector3<T> &n, Math::Vector3<T> &d);
|
||||
|
||||
|
||||
/**
|
||||
* Calculates principle axis P and T from fault normal and slip vector
|
||||
*/
|
||||
template <typename T>
|
||||
bool nd2pa(const Math::Vector3<T> &n, const Math::Vector3<T> &d, Math::Vector3<T> &t, Math::Vector3<T> &p);
|
||||
|
||||
|
||||
/**
|
||||
* Calculates the nodal plane from the normal and slip vector of the fault plane.
|
||||
* The result is in radiants.
|
||||
*/
|
||||
template <typename T>
|
||||
bool nd2np(const Math::Vector3<T> &n, const Math::Vector3<T> &d, NODAL_PLANE &np);
|
||||
|
||||
|
||||
/**
|
||||
* Calculates the the normal and slip vector of the fault plane
|
||||
* from a nodal plane. The nodal plane attributes are expected to be
|
||||
* in radiants.
|
||||
*/
|
||||
template <typename T>
|
||||
bool np2nd(const NODAL_PLANE &np, Math::Vector3<T> &n, Math::Vector3<T> &d);
|
||||
|
||||
|
||||
/**
|
||||
* Calculates a tensor from fault normal and slip vector
|
||||
*/
|
||||
template <typename T>
|
||||
bool nd2tensor(const Math::Vector3<T> &n, const Math::Vector3<T> &d, Math::Tensor2S<T> &t);
|
||||
|
||||
|
||||
/**
|
||||
* Calculates a tensor from a nodal plane
|
||||
*/
|
||||
template <typename T>
|
||||
bool np2tensor(const NODAL_PLANE &np, Math::Tensor2S<T> &t);
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Converts a nodal plane from radiants to degree
|
||||
*/
|
||||
void np2deg(NODAL_PLANE &np);
|
||||
|
||||
|
||||
/**
|
||||
* Computes the double couples from normal and slip vector
|
||||
*/
|
||||
template <typename T>
|
||||
bool nd2dc(const Math::Vector3<T> &n, const Math::Vector3<T> &d, NODAL_PLANE *NP1, NODAL_PLANE *NP2);
|
||||
|
||||
|
||||
template <typename T>
|
||||
void rtp2tensor(T mrr, T mtt, T mpp, T mrt, T mrp, T mtp, Math::Tensor2S<T> &tensor);
|
||||
|
||||
|
||||
template <typename T>
|
||||
double spectral2clvd(Math::Tensor2S<T> &tensor);
|
||||
|
||||
|
||||
template <typename T>
|
||||
void spectral2axis(const Math::Spectral2S<T> &spec,
|
||||
AXIS &t, AXIS &n, AXIS &p, int expo);
|
||||
|
||||
template <typename T>
|
||||
void axis2spectral(const AXIS &t, const AXIS &n, const AXIS &p, int expo,
|
||||
Math::Spectral2S<T> &spec);
|
||||
|
||||
template <typename T>
|
||||
void axis2tensor(const AXIS &t, const AXIS &n, const AXIS &p, int expo,
|
||||
Math::Tensor2S<T> &tensor);
|
||||
|
||||
template <typename T1, typename T2>
|
||||
void spectral2matrix(const Math::Spectral2S<T1> &t, Math::Matrix3<T2> &m);
|
||||
|
||||
template <typename T1, typename T2>
|
||||
bool tensor2matrix(const Math::Tensor2S<T1> &t, Math::Matrix3<T2> &m);
|
||||
|
||||
|
||||
#include "conversions.ipp"
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
258
include/seiscomp/math/conversions.ipp
Normal file
258
include/seiscomp/math/conversions.ipp
Normal file
@ -0,0 +1,258 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
template <typename T>
|
||||
void vec2angle(const Math::Vector3<T> &v, T &strike, T &dip) {
|
||||
if ( fabs(v.x) < 0.0001 && fabs(v.y) < 0.0001 ) {
|
||||
strike = 0;
|
||||
dip = 0;
|
||||
}
|
||||
else {
|
||||
// TODO: Is this correct?
|
||||
if ( v.z < 0 ) {
|
||||
strike = (T)(atan2(-v.y,-v.x));
|
||||
dip = (T)(acos(-v.z));
|
||||
}
|
||||
else {
|
||||
strike = (T)(atan2(v.y,v.x));
|
||||
dip = (T)(acos(v.z));
|
||||
}
|
||||
|
||||
if ( strike < 0 ) strike += (2*M_PI);
|
||||
|
||||
strike = rad2deg(strike);
|
||||
}
|
||||
|
||||
dip = rad2deg(HALF_PI-dip);
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
void angle2vec(T strike, T dip, Math::Vector3<T> &v) {
|
||||
T rdip = deg2rad(dip);
|
||||
T rstrike = deg2rad(strike);
|
||||
|
||||
v.x = cos(rdip)*cos(rstrike);
|
||||
v.y = cos(rdip)*sin(rstrike);
|
||||
v.z = sin(rdip);
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool pa2nd(const Math::Vector3<T> &t, const Math::Vector3<T> &p, Math::Vector3<T> &n, Math::Vector3<T> &d) {
|
||||
n = t + p;
|
||||
n.normalize();
|
||||
|
||||
d = t - p;
|
||||
d.normalize();
|
||||
|
||||
if ( n.z > 0 ) {
|
||||
n *= -1;
|
||||
d *= -1;
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool nd2pa(const Math::Vector3<T> &n, const Math::Vector3<T> &d, Math::Vector3<T> &t, Math::Vector3<T> &p) {
|
||||
p = n - d;
|
||||
p.normalize();
|
||||
if ( p.z < 0 )
|
||||
p *= -1;
|
||||
|
||||
t = n + d;
|
||||
t.normalize();
|
||||
if ( t.z < 0 )
|
||||
t *= -1;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool nd2tensor(const Math::Vector3<T> &n, const Math::Vector3<T> &d, Math::Tensor2S<T> &t) {
|
||||
t._11 = 2*d.x*n.x;
|
||||
t._12 = d.x*n.y+d.y*n.x;
|
||||
t._13 = d.x*n.z+d.z*n.x;
|
||||
t._22 = 2*d.y*n.y;
|
||||
t._23 = d.y*n.z+d.z*n.y;
|
||||
t._33 = 2*d.z*n.z;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool np2tensor(const NODAL_PLANE &np, Math::Tensor2S<T> &t) {
|
||||
Math::Vector3<T> n, d;
|
||||
np2nd(np, n, d);
|
||||
nd2tensor(n, d, t);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool nd2np(const Math::Vector3<T> &n, const Math::Vector3<T> &d, NODAL_PLANE &np) {
|
||||
if ( fabs(n.z) == 1 ) {
|
||||
np.dip = 0;
|
||||
np.str = 0;
|
||||
np.rake = atan2(-d.y,d.x);
|
||||
}
|
||||
else {
|
||||
np.dip = acos(-n.z);
|
||||
np.str = atan2(-n.x, n.y);
|
||||
np.rake = atan2(-d.z/sin(np.dip),d.x*cos(np.str)+d.y*sin(np.str));
|
||||
}
|
||||
|
||||
np.str = (double)fmod(np.str+(2*M_PI),(2*M_PI));
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool np2nd(const NODAL_PLANE &np, Math::Vector3<T> &n, Math::Vector3<T> &d) {
|
||||
T rdip = deg2rad(np.dip), rstr = deg2rad(np.str), rrake = deg2rad(np.rake);
|
||||
|
||||
n.x = -sin(rdip)*sin(rstr);
|
||||
n.y = sin(rdip)*cos(rstr);
|
||||
n.z = -cos(rdip);
|
||||
|
||||
d.x = cos(rrake)*cos(rstr)+cos(rdip)*sin(rrake)*sin(rstr);
|
||||
d.y = cos(rrake)*sin(rstr)-cos(rdip)*sin(rrake)*cos(rstr);
|
||||
d.z = -sin(rdip) *sin(rrake);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool nd2dc(const Math::Vector3<T> &n, const Math::Vector3<T> &d, NODAL_PLANE *NP1, NODAL_PLANE *NP2) {
|
||||
// Compute DC from normal and slip vector
|
||||
nd2np(n, d, *NP1);
|
||||
|
||||
Math::Vector3<T> n1(d);
|
||||
Math::Vector3<T> d1(n);
|
||||
|
||||
// Rotate n-d system be 180 degree
|
||||
if ( n1.z > 0 ) {
|
||||
n1 *= -1;
|
||||
d1 *= -1;
|
||||
}
|
||||
|
||||
// Compute DC from normal (rotated slip) and slip (rotated normal) vector
|
||||
nd2np(n1, d1, *NP2);
|
||||
|
||||
// Convert both nodal planes from radiant to degree
|
||||
np2deg(*NP1);
|
||||
np2deg(*NP2);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
void rtp2tensor(T mrr, T mtt, T mpp, T mrt, T mrp, T mtp, Math::Tensor2S<T> &tensor) {
|
||||
tensor._33 = mrr;
|
||||
tensor._11 = mtt;
|
||||
tensor._22 = mpp;
|
||||
tensor._13 = mrt;
|
||||
tensor._23 = -mrp;
|
||||
tensor._12 = -mtp;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
double spectral2clvd(const Math::Spectral2S<T> &spec) {
|
||||
double mean = (spec.a1 + spec.a2 + spec.a3) / 3;
|
||||
double v[3] = {fabs(spec.a1-mean), fabs(spec.a2-mean), fabs(spec.a3-mean)};
|
||||
std::sort(v, v+3);
|
||||
|
||||
double eps = v[0] / v[2];
|
||||
return 200.0*eps;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
void spectral2axis(const Math::Spectral2S<T> &spec,
|
||||
AXIS &t, AXIS &n, AXIS &p, int expo) {
|
||||
t.val = spec.a1; t.expo = expo;
|
||||
n.val = spec.a2; n.expo = expo;
|
||||
p.val = spec.a3; p.expo = expo;
|
||||
|
||||
vec2angle(spec.n1, t.str, t.dip);
|
||||
vec2angle(spec.n2, n.str, n.dip);
|
||||
vec2angle(spec.n3, p.str, p.dip);
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
void axis2spectral(const AXIS &t, const AXIS &n, const AXIS &p, int expo,
|
||||
Math::Spectral2S<T> &spec) {
|
||||
angle2vec(t.str, t.dip, spec.n1);
|
||||
angle2vec(n.str, n.dip, spec.n2);
|
||||
angle2vec(p.str, p.dip, spec.n3);
|
||||
|
||||
spec.a1 = t.val * pow((T)10, (T)t.expo-expo);
|
||||
spec.a2 = n.val * pow((T)10, (T)n.expo-expo);
|
||||
spec.a3 = p.val * pow((T)10, (T)p.expo-expo);
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
void axis2tensor(const AXIS &t, const AXIS &n, const AXIS &p, int expo,
|
||||
Math::Tensor2S<T> &tensor) {
|
||||
Math::Spectral2S<T> spec;
|
||||
axis2spectral(t, n, p, expo, spec);
|
||||
spec.compose(tensor);
|
||||
}
|
||||
|
||||
|
||||
template <typename T1, typename T2>
|
||||
void spectral2matrix(const Math::Spectral2S<T1> &spec, Math::Matrix3<T2> &m) {
|
||||
Math::Vector3<T1> x3, x1, x2;
|
||||
|
||||
// Convert principle axis to tensor orientation
|
||||
pa2nd(spec.n1, spec.n3, x3, x1);
|
||||
|
||||
// Calculate null axis
|
||||
x2.cross(x3, x1);
|
||||
|
||||
// Slip axis
|
||||
m.setColumn(0, Math::Vector3<T2>(x1.x,x1.y,x1.z));
|
||||
// Dip axis
|
||||
m.setColumn(1, Math::Vector3<T2>(x2.x,x2.y,x2.z));
|
||||
// Fault normal
|
||||
m.setColumn(2, Math::Vector3<T2>(x3.x,x3.y,x3.z));
|
||||
}
|
||||
|
||||
|
||||
template <typename T1, typename T2>
|
||||
bool tensor2matrix(const Math::Tensor2S<T1> &t, Math::Matrix3<T2> &m) {
|
||||
Math::Spectral2S<T1> spec;
|
||||
if ( !spec.spect(t) ) return false;
|
||||
|
||||
spec.sort();
|
||||
spectral2matrix(spec, m);
|
||||
|
||||
return true;
|
||||
}
|
||||
125
include/seiscomp/math/coord.h
Normal file
125
include/seiscomp/math/coord.h
Normal file
@ -0,0 +1,125 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_GEO_COORD_H
|
||||
#define SEISCOMP_MATH_GEO_COORD_H
|
||||
|
||||
#include <seiscomp/core/baseobject.h>
|
||||
#include <string>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Geo {
|
||||
|
||||
|
||||
template<typename T>
|
||||
struct Coord : public Core::BaseObject {
|
||||
DECLARE_SERIALIZATION;
|
||||
|
||||
typedef T ValueType;
|
||||
|
||||
Coord();
|
||||
Coord(T lat_, T lon_);
|
||||
|
||||
void set(T lat_, T lon_);
|
||||
|
||||
T latitude() const { return lat; }
|
||||
T longitude() const { return lon; }
|
||||
|
||||
bool operator==(const Coord<T> &other) const;
|
||||
bool operator!=(const Coord<T> &other) const;
|
||||
|
||||
T lat;
|
||||
T lon;
|
||||
};
|
||||
|
||||
typedef Coord<float> CoordF;
|
||||
typedef Coord<double> CoordD;
|
||||
|
||||
|
||||
template<typename T>
|
||||
class NamedCoord : public Coord<T> {
|
||||
public:
|
||||
NamedCoord();
|
||||
NamedCoord(const std::string& name, T lat_, T lon_);
|
||||
|
||||
~NamedCoord();
|
||||
|
||||
using Coord<T>::set;
|
||||
void set(const std::string& name, T lat_, T lon_);
|
||||
|
||||
public:
|
||||
void setName(const std::string& name);
|
||||
const std::string& name() const;
|
||||
|
||||
void serialize(Core::BaseObject::Archive& ar);
|
||||
|
||||
private:
|
||||
std::string _name;
|
||||
};
|
||||
|
||||
|
||||
typedef NamedCoord<float> NamedCoordF;
|
||||
typedef NamedCoord<double> NamedCoordD;
|
||||
|
||||
|
||||
template<typename T>
|
||||
class City : public NamedCoord<T> {
|
||||
public:
|
||||
City();
|
||||
City(const std::string& name,
|
||||
T lat_, T lon_, size_t population);
|
||||
|
||||
City(const std::string& name, const std::string& countryID,
|
||||
T lat_, T lon_, size_t population);
|
||||
|
||||
~City();
|
||||
|
||||
public:
|
||||
//! Population in thousands
|
||||
void setPopulation(double population);
|
||||
double population() const;
|
||||
|
||||
void setCountryID(const std::string &);
|
||||
const std::string &countryID() const;
|
||||
|
||||
void setCategory(std::string &);
|
||||
const std::string &category() const;
|
||||
|
||||
void serialize(Core::BaseObject::Archive& ar);
|
||||
|
||||
private:
|
||||
std::string _countryID;
|
||||
double _population;
|
||||
std::string _category;
|
||||
};
|
||||
|
||||
|
||||
typedef City<float> CityF;
|
||||
typedef City<double> CityD;
|
||||
|
||||
|
||||
} // of ns Geo
|
||||
} // of ns Math
|
||||
} // of ns Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
93
include/seiscomp/math/decomp.ipp
Normal file
93
include/seiscomp/math/decomp.ipp
Normal file
@ -0,0 +1,93 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
#include<vector>
|
||||
#include<math.h>
|
||||
|
||||
#define VPVS(sigma) sqrt((2.-2.*(sigma))/(1.-2.*(sigma)))
|
||||
|
||||
namespace Seiscomp
|
||||
{
|
||||
|
||||
namespace Math
|
||||
{
|
||||
|
||||
namespace Filtering
|
||||
{
|
||||
|
||||
template<typename TYPE>
|
||||
int rotate(std::vector<TYPE> &f1, std::vector<TYPE> &f2, double phi)
|
||||
{
|
||||
double m11 = cos(phi*M_PI/180.), m12 = sin(phi*M_PI/180.),
|
||||
m21 = -m12, m22 = m11, x, y;
|
||||
|
||||
if (f1.size() != f2.size())
|
||||
throw AlignmentError("rotate(): vector's not aligned");
|
||||
|
||||
TYPE *p1 = &f1[0], *p2 = &f2[0];
|
||||
int n = f1.size();
|
||||
|
||||
while (n--) {
|
||||
x = (*p1)*m11 + (*p2)*m12;
|
||||
y = (*p1)*m21 + (*p2)*m22;
|
||||
(*p1++) = (TYPE) x;
|
||||
(*p2++) = (TYPE) y;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
int decompose(std::vector<TYPE> &f1, std::vector<TYPE> &f2,
|
||||
double p, double vs, double sigma)
|
||||
{
|
||||
double vp = vs*VPVS(sigma), qa, qb;
|
||||
double m11, m12, m21, m22, x, y;
|
||||
|
||||
if (p >= 1./vp)
|
||||
return -1;
|
||||
|
||||
if (f1.size() != f2.size())
|
||||
throw AlignmentError("decompose(): vector's not aligned");
|
||||
|
||||
TYPE *p1 = &f1[0], *p2 = &f2[0];
|
||||
int n = f1.size();
|
||||
|
||||
qa = sqrt (1./(vp*vp)-p*p);
|
||||
qb = sqrt (1./(vs*vs)-p*p);
|
||||
m11 = -(2*vs*vs*p*p-1.)/(vp*qa);
|
||||
m12 = 2.*p*vs*vs/vp;
|
||||
m21 = -2.*p*vs;
|
||||
m22 = (1.-2.*vs*vs*p*p)/(vs*qb);
|
||||
|
||||
while (n--) {
|
||||
x = (*p1)*m11 + (*p2)*m12;
|
||||
y = (*p1)*m21 + (*p2)*m22;
|
||||
(*p1++) = (TYPE) x;
|
||||
(*p2++) = (TYPE) y;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
} // namespace Seiscomp::Math::Filter
|
||||
|
||||
} // namespace Seiscomp::Math
|
||||
|
||||
} // namespace Seiscomp
|
||||
90
include/seiscomp/math/fft.h
Normal file
90
include/seiscomp/math/fft.h
Normal file
@ -0,0 +1,90 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FFT_H
|
||||
#define SEISCOMP_MATH_FFT_H
|
||||
|
||||
|
||||
#include <complex>
|
||||
#include <vector>
|
||||
#include <seiscomp/core/typedarray.h>
|
||||
#include <seiscomp/math/math.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
typedef std::complex<double> Complex;
|
||||
typedef std::vector<Math::Complex> ComplexArray;
|
||||
|
||||
|
||||
template <typename T>
|
||||
void fft(ComplexArray &spec, int n, const T *data);
|
||||
|
||||
template <typename T>
|
||||
void fft(ComplexArray &spec, const std::vector<T> &data) {
|
||||
fft(spec, static_cast<int>(data.size()), &data[0]);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void fft(ComplexArray &spec, const TypedArray<T> &data) {
|
||||
fft(spec, data.impl());
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void fft(Seiscomp::ComplexDoubleArray &spec, const TypedArray<T> &data) {
|
||||
fft(spec.impl(), data.impl());
|
||||
}
|
||||
|
||||
inline void fft(Seiscomp::ComplexDoubleArray &spec, const DoubleArray &data) {
|
||||
fft(spec.impl(), data.impl());
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
void ifft(int n, T *out, ComplexArray &spec);
|
||||
|
||||
template <typename T>
|
||||
void ifft(std::vector<T> &out, ComplexArray &spec) {
|
||||
ifft(static_cast<int>(out.size()), &out[0], spec);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void ifft(TypedArray<T> &out, ComplexArray &spec) {
|
||||
ifft(out.impl(), &out[0], spec);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void ifft(TypedArray<T> &out, Seiscomp::ComplexDoubleArray &spec) {
|
||||
ifft(out.impl(), &out[0], spec.impl());
|
||||
}
|
||||
|
||||
inline void ifft(DoubleArray &out, Seiscomp::ComplexDoubleArray &spec) {
|
||||
ifft(out.impl(), spec.impl());
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
182
include/seiscomp/math/filter.h
Normal file
182
include/seiscomp/math/filter.h
Normal file
@ -0,0 +1,182 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_FILTERING_FILTER_H
|
||||
#define SEISCOMP_FILTERING_FILTER_H
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/core/baseobject.h>
|
||||
#include <seiscomp/core/interfacefactory.h>
|
||||
#include <seiscomp/core/typedarray.h>
|
||||
#include <seiscomp/core/datetime.h>
|
||||
#include <seiscomp/core.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
class SC_SYSTEM_CORE_API AlignmentError : public std::exception {
|
||||
public:
|
||||
AlignmentError(const char *txt) { _txt = txt; }
|
||||
const char* what() const throw() { return _txt; }
|
||||
private:
|
||||
const char *_txt;
|
||||
};
|
||||
|
||||
// virtual base class that all filter classes should be derived from
|
||||
|
||||
template<typename TYPE>
|
||||
class SC_SYSTEM_CORE_API InPlaceFilter : public Core::BaseObject {
|
||||
public:
|
||||
virtual ~InPlaceFilter() { }
|
||||
|
||||
public:
|
||||
/**
|
||||
* Sets the start time of the data to be filtered.
|
||||
* This is not handled in the base class but can be used
|
||||
* in derived classes to apply filters depending on
|
||||
* this information. This function should be called before
|
||||
* apply is called the first time.
|
||||
*/
|
||||
virtual void setStartTime(const Core::Time &time) {}
|
||||
|
||||
/**
|
||||
* Sets the streamID of the data to be filtered.
|
||||
* This is not handled in the base class but can be used
|
||||
* in derived classes to apply filters depending on
|
||||
* this information. This function should be called before
|
||||
* apply is called the first time.
|
||||
*/
|
||||
virtual void setStreamID(const std::string &net,
|
||||
const std::string &sta,
|
||||
const std::string &loc,
|
||||
const std::string &cha) {}
|
||||
|
||||
virtual void setSamplingFrequency(double fsamp) = 0;
|
||||
|
||||
/**
|
||||
* Sets the filter parameters
|
||||
* @param n The number of given parameters
|
||||
* @param params The parameter list.
|
||||
* @return Positive value: the number of required parameters or,
|
||||
* negative value: a value error at parameter position abs(return)-1
|
||||
*/
|
||||
virtual int setParameters(int n, const double *params) = 0;
|
||||
|
||||
virtual void apply(int n, TYPE *inout) = 0;
|
||||
|
||||
void apply(std::vector<TYPE> &f) {
|
||||
apply(f.size(), &f[0]);
|
||||
}
|
||||
|
||||
void apply(TypedArray<TYPE> &arr) {
|
||||
apply(arr.size(), arr.typedData());
|
||||
}
|
||||
|
||||
// EXPERIMENTAL
|
||||
virtual void handleGap(int n=0) {}
|
||||
|
||||
//! Creates a new filter instance having the same
|
||||
//! type as the called instance. The configuration
|
||||
//! parameters will be copied, too but not the internal
|
||||
//! filter state depending on former input data
|
||||
virtual InPlaceFilter<TYPE>* clone() const = 0;
|
||||
|
||||
//! Creates a new filter by name. The user is responsible to
|
||||
//! release the memory
|
||||
//! pointed to by the return value.
|
||||
//! \return The corresponding filter or nullptr
|
||||
static InPlaceFilter<TYPE> *Create(const std::string &strFilter,
|
||||
std::string *strError = nullptr);
|
||||
};
|
||||
|
||||
|
||||
DEFINE_TEMPLATE_INTERFACE_FACTORY(InPlaceFilter);
|
||||
|
||||
#define REGISTER_INPLACE_FILTER(Class, Service) \
|
||||
Seiscomp::Core::Generic::InterfaceFactory<Seiscomp::Math::Filtering::InPlaceFilter<float>, Class<float> > __##Class##FloatInterfaceFactory__(Service); \
|
||||
Seiscomp::Core::Generic::InterfaceFactory<Seiscomp::Math::Filtering::InPlaceFilter<double>, Class<double> > __##Class##DoubleInterfaceFactory__(Service)
|
||||
|
||||
#define REGISTER_INPLACE_FILTER2(Class, FactoryPrefix, Service) \
|
||||
Seiscomp::Core::Generic::InterfaceFactory<Seiscomp::Math::Filtering::InPlaceFilter<float>, Class<float> > __##Class##FactoryPrefix##FloatInterfaceFactory__(Service); \
|
||||
Seiscomp::Core::Generic::InterfaceFactory<Seiscomp::Math::Filtering::InPlaceFilter<double>, Class<double> > __##Class##FactoryPrefix##DoubleInterfaceFactory__(Service)
|
||||
|
||||
#define INSTANTIATE_INPLACE_FILTER(Class, APIDef) \
|
||||
template class APIDef Class<float>; \
|
||||
template class APIDef Class<double>
|
||||
|
||||
|
||||
//! Dummy filter that does nothing but can be used in more complex
|
||||
//! filters like Op2Filter.
|
||||
template<typename TYPE>
|
||||
class SC_SYSTEM_CORE_API SelfFilter : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
SelfFilter() {}
|
||||
|
||||
public:
|
||||
virtual void setSamplingFrequency(double fsamp) {}
|
||||
virtual int setParameters(int n, const double *params) { return 0; }
|
||||
|
||||
virtual void apply(int n, TYPE *inout) {}
|
||||
|
||||
virtual InPlaceFilter<TYPE>* clone() const { return new SelfFilter(); }
|
||||
};
|
||||
|
||||
|
||||
template <class TYPE>
|
||||
int minmax(std::vector<TYPE> const &f, int i1, int i2, int *imax, TYPE *fmax);
|
||||
|
||||
template <class TYPE>
|
||||
int find_max(std::vector<TYPE> const &f, int i1, int i2, int *imax, TYPE *fmax);
|
||||
|
||||
template <typename TYPE>
|
||||
int rotate(std::vector<TYPE> &f1, std::vector<TYPE> &f2, double phi);
|
||||
|
||||
template<typename TYPE>
|
||||
int decompose(std::vector<TYPE> &f1, std::vector<TYPE> &f2,
|
||||
double p, double vs, double sigma);
|
||||
|
||||
template<typename TYPE>
|
||||
void hilbert_transform(std::vector<TYPE> &trace, int direction=1);
|
||||
template<typename TYPE>
|
||||
void envelope(std::vector<TYPE> &trace);
|
||||
|
||||
SC_SYSTEM_CORE_API long next_power_of_2(long);
|
||||
|
||||
template<typename TYPE>
|
||||
void cosRamp(std::vector<TYPE> &ramp, TYPE f1, TYPE f2);
|
||||
|
||||
template <typename T>
|
||||
void cosRamp(size_t n, T *inout, size_t istart, size_t iend, size_t estart, size_t eend);
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filter
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#include<seiscomp/math/hilbert.ipp>
|
||||
#include<seiscomp/math/minmax.ipp>
|
||||
#include<seiscomp/math/decomp.ipp>
|
||||
|
||||
|
||||
#endif
|
||||
53
include/seiscomp/math/filter/abs.h
Normal file
53
include/seiscomp/math/filter/abs.h
Normal file
@ -0,0 +1,53 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_ABS_H
|
||||
#define SEISCOMP_MATH_FILTER_ABS_H
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename T>
|
||||
class AbsFilter : public InPlaceFilter<T> {
|
||||
public:
|
||||
AbsFilter();
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
void apply(int n, T *inout) override;
|
||||
|
||||
InPlaceFilter<T> *clone() const override;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
68
include/seiscomp/math/filter/average.h
Normal file
68
include/seiscomp/math/filter/average.h
Normal file
@ -0,0 +1,68 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_RCA_H
|
||||
#define SEISCOMP_MATH_FILTER_RCA_H
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class Average : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
Average(double timeSpan /*sec*/ = 1.0, double fsamp = 0.0);
|
||||
|
||||
public:
|
||||
void setLength(double timeSpan);
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
// Resets the filter values
|
||||
void reset();
|
||||
|
||||
private:
|
||||
double _timeSpan;
|
||||
double _fsamp;
|
||||
double _oocount;
|
||||
int _sampleCount;
|
||||
int _index;
|
||||
double _lastSum;
|
||||
bool _firstSample;
|
||||
std::vector<TYPE> _buffer;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
|
||||
177
include/seiscomp/math/filter/biquad.h
Normal file
177
include/seiscomp/math/filter/biquad.h
Normal file
@ -0,0 +1,177 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_BIQUAD_H
|
||||
#define SEISCOMP_MATH_FILTER_BIQUAD_H
|
||||
|
||||
|
||||
#include <vector>
|
||||
#include <string>
|
||||
#include <ostream>
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
namespace IIR {
|
||||
|
||||
|
||||
/**
|
||||
* Coefficients of a Biquad.
|
||||
*
|
||||
* Note the order of the coefficients passed to the constructor and the
|
||||
* setter: the numerator coefficients (the b's) followed by the denominator
|
||||
* coefficients (the a's).
|
||||
*
|
||||
* Also note that the coeffients are normalized and a0 is *always* assumed
|
||||
* to be 1. This is not checked.
|
||||
*/
|
||||
struct BiquadCoefficients {
|
||||
BiquadCoefficients(double b0 = 0, double b1 = 0, double b2 = 0,
|
||||
double a0 = 1, double a1 = 0, double a2 = 0);
|
||||
BiquadCoefficients(BiquadCoefficients const &bq);
|
||||
|
||||
void set(double b0, double b1, double b2,
|
||||
double a0, double a1, double a2);
|
||||
|
||||
// filter coefficients
|
||||
double b0, b1, b2; // numerator coefficients
|
||||
double a0, a1, a2; // denominator coefficients (a0==1)
|
||||
};
|
||||
|
||||
|
||||
typedef std::vector<BiquadCoefficients> Biquads;
|
||||
|
||||
|
||||
std::ostream &operator<<(std::ostream &os, const BiquadCoefficients &biq);
|
||||
|
||||
|
||||
/**
|
||||
* Template class to implement the filter interface for a single set of
|
||||
* biquad coefficients.
|
||||
*/
|
||||
template<typename TYPE>
|
||||
class Biquad : public InPlaceFilter<TYPE> {
|
||||
// ------------------------------------------------------------------
|
||||
// X'truction
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
Biquad(double b0 = 0, double b1 = 0, double b2 = 0,
|
||||
double a0 = 1, double a1 = 0, double a2 = 0);
|
||||
Biquad(const BiquadCoefficients &bq);
|
||||
Biquad(const Biquad &bq);
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Public methods
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
// reset the filter by erasing memory of past samples
|
||||
void reset();
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// InplaceFilter interface
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
|
||||
// Implement InPlaceFilter interface with default values
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Public members
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
// filter coefficients
|
||||
BiquadCoefficients coefficients;
|
||||
// filter memory
|
||||
double v1, v2;
|
||||
|
||||
};
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class BiquadCascade : public InPlaceFilter<TYPE> {
|
||||
// ------------------------------------------------------------------
|
||||
// X'truction
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
BiquadCascade();
|
||||
BiquadCascade(const BiquadCascade &other);
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Public methods
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
// number of biquads comprising the cascade
|
||||
int size() const;
|
||||
|
||||
// resets the filter, i.e. erases the filter memory
|
||||
void reset();
|
||||
|
||||
// append a single biquad to this cascade
|
||||
void append(const Biquad<TYPE> &biq);
|
||||
|
||||
void set(const Biquads &biquads);
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// InplaceFilter interface
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
void setSamplingFrequency(double /*fsamp*/) override {}
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Protected members
|
||||
// ------------------------------------------------------------------
|
||||
protected:
|
||||
std::vector< Biquad<TYPE> > _biq;
|
||||
|
||||
template <typename T>
|
||||
friend std::ostream &operator<<(std::ostream &os, const BiquadCascade<T> &b);
|
||||
};
|
||||
|
||||
|
||||
template <typename T>
|
||||
std::ostream &operator<<(std::ostream &os, const BiquadCascade<T> &b);
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering::IIR
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
136
include/seiscomp/math/filter/biquad.ipp
Normal file
136
include/seiscomp/math/filter/biquad.ipp
Normal file
@ -0,0 +1,136 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
// This file is included by "biquad.cpp"
|
||||
|
||||
template<typename TYPE>
|
||||
Biquad<TYPE>::Biquad(double b0, double b1, double b2,
|
||||
double a0, double a1, double a2)
|
||||
: InPlaceFilter<TYPE>()
|
||||
, coefficients(b0, b1, b2, a0, a1, a2) {
|
||||
v1 = v2 = 0.;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
Biquad<TYPE>::Biquad(const BiquadCoefficients &coeff)
|
||||
: InPlaceFilter<TYPE>()
|
||||
, coefficients(coeff) {
|
||||
v1 = v2 = 0.;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
Biquad<TYPE>::Biquad(const Biquad<TYPE> &other)
|
||||
: InPlaceFilter<TYPE>()
|
||||
, coefficients(other.coefficients) {
|
||||
v1 = v2 = 0.;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
void Biquad<TYPE>::apply(int n, TYPE *inout) {
|
||||
// This is the direct form 2 implementation according to
|
||||
// https://en.wikipedia.org/wiki/Digital_biquad_filter#Direct_form_2
|
||||
TYPE *ff = inout;
|
||||
for ( int i = 0; i < n; ++i ) {
|
||||
|
||||
// a0 is assumed to be 1
|
||||
double v0 = ff[i] - coefficients.a1*v1 - coefficients.a2*v2;
|
||||
ff[i] = TYPE(coefficients.b0*v0 + coefficients.b1*v1 + coefficients.b2*v2);
|
||||
v2 = v1; v1 = v0;
|
||||
}
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
InPlaceFilter<TYPE>* Biquad<TYPE>::clone() const {
|
||||
return new Biquad<TYPE>(coefficients);
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
void Biquad<TYPE>::reset() {
|
||||
v1 = v2 = 0.;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
void Biquad<TYPE>::setSamplingFrequency(double fsamp) {}
|
||||
|
||||
template<typename TYPE>
|
||||
int Biquad<TYPE>::setParameters(int n, const double *params) {
|
||||
if ( n != 6 ) return 6;
|
||||
|
||||
reset();
|
||||
coefficients.set(params[0], params[1], params[2],
|
||||
params[3], params[4], params[5]);
|
||||
|
||||
return n;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
BiquadCascade<TYPE>::BiquadCascade() {}
|
||||
|
||||
template<typename TYPE>
|
||||
BiquadCascade<TYPE>::BiquadCascade(const BiquadCascade &other) {
|
||||
_biq = other._biq;
|
||||
reset();
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
int BiquadCascade<TYPE>::size() const { return _biq.size(); }
|
||||
|
||||
template<typename TYPE>
|
||||
void BiquadCascade<TYPE>::apply(int n, TYPE *inout) {
|
||||
for ( Biquad<TYPE> &biq : _biq )
|
||||
biq.apply(n, inout);
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
int BiquadCascade<TYPE>::setParameters(int /*n*/, const double* /*params*/) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
InPlaceFilter<TYPE>* BiquadCascade<TYPE>::clone() const {
|
||||
return new BiquadCascade<TYPE>(*this);
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
void BiquadCascade<TYPE>::reset() {
|
||||
for ( Biquad<TYPE> &biq : _biq )
|
||||
biq.reset();
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
void BiquadCascade<TYPE>::append(Biquad<TYPE> const &biq) {
|
||||
_biq.push_back(biq);
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
void BiquadCascade<TYPE>::set(const Biquads &biquads) {
|
||||
_biq.clear();
|
||||
for ( const BiquadCoefficients &biq : biquads )
|
||||
_biq.push_back(biq);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
std::ostream &operator<<(std::ostream &os, const BiquadCascade<T> &b) {
|
||||
int i = 0;
|
||||
for ( const Biquad<T> &biq : b._biq )
|
||||
os << "Biquad #" << ++i << std::endl << biq.coefficients;
|
||||
return os;
|
||||
}
|
||||
|
||||
99
include/seiscomp/math/filter/bpenv.h
Normal file
99
include/seiscomp/math/filter/bpenv.h
Normal file
@ -0,0 +1,99 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_FILTERING_BPENV_H
|
||||
#define SEISCOMP_FILTERING_BPENV_H
|
||||
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
#include <seiscomp/math/filter/butterworth.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
// Band pass and envelope filter
|
||||
//
|
||||
// This is a recursive band pass and envelope filter combination. Strictly
|
||||
// speaking the envelope would require frequency-domain filtering because it
|
||||
// requires the Hilbert transform of the signal. We cannot compute the Hilbert
|
||||
// transform using time-domain recursive filters. By applying a rather narrow
|
||||
// band pass to the data, however, we can compute an approximate Hilbert
|
||||
// transform. This is because for a narrow-band signal the Hilbert transform
|
||||
// will be the scaled time derivative. This actually works quite well in
|
||||
// practice as long as the band pass is not too wide and has relatively steep
|
||||
// corners to block energy from outside the passband. A band width of one
|
||||
// octave and a filter order of four was found to work well in most cases.
|
||||
// It is therefore highly recommended to stick to these defaults and only set
|
||||
// the center frequency.
|
||||
//
|
||||
// This filter may be used in a filter chain. Its name is BPENV and it accepts
|
||||
// up to three arguments:
|
||||
// - center frequency in Hz (must be specified)
|
||||
// - bandwidth in octaves (default: 1)
|
||||
// - filter order (default: 4)
|
||||
//
|
||||
// Author:
|
||||
// Joachim Saul, GFZ Potsdam
|
||||
// saul@gfz-potsdam.de
|
||||
|
||||
|
||||
template <typename T>
|
||||
class BandPassEnvelope : public InPlaceFilter<T> {
|
||||
public:
|
||||
BandPassEnvelope(double centerFrequency=1, double bandwidthOctaves=1, int order=4, double fsamp=0);
|
||||
|
||||
// Apply filter to data vector in place
|
||||
void apply(int ndata, T *data) override;
|
||||
|
||||
// Apply a full reset
|
||||
void reset();
|
||||
|
||||
// Set the sampling frequency in Hz. Allows delayed
|
||||
// initialization when the data arrive.
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<T> *clone() const override;
|
||||
|
||||
|
||||
private:
|
||||
double _centerFrequency;
|
||||
double _bandwidthOctaves;
|
||||
int _order;
|
||||
double _samplingFrequency;
|
||||
double _K, _yp;
|
||||
|
||||
using FilterPtr = typename Core::SmartPointer<Math::Filtering::InPlaceFilter<T>>::Impl;
|
||||
FilterPtr _bp;
|
||||
|
||||
// Flag to indicate some pending initialization
|
||||
bool _afterReset;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
131
include/seiscomp/math/filter/butterworth.h
Normal file
131
include/seiscomp/math/filter/butterworth.h
Normal file
@ -0,0 +1,131 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_FILTERING_IIR_BUTTERWORTH_H
|
||||
#define SEISCOMP_FILTERING_IIR_BUTTERWORTH_H
|
||||
|
||||
|
||||
#include <seiscomp/math/filter/biquad.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
namespace IIR {
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
class ButterworthLowpass : public BiquadCascade<TYPE> {
|
||||
public:
|
||||
ButterworthLowpass(int order = 3, double fmax = 0.7, double fsamp=0);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
private:
|
||||
int _order;
|
||||
double _fmax, _fsamp;
|
||||
};
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
class ButterworthHighpass : public BiquadCascade<TYPE> {
|
||||
public:
|
||||
ButterworthHighpass(int order = 3, double fmin = 2.0, double fsamp=0);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
private:
|
||||
int _order;
|
||||
double _fmin, _fsamp;
|
||||
};
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
class ButterworthBandpass : public BiquadCascade<TYPE> {
|
||||
public:
|
||||
ButterworthBandpass(int order = 3, double fmin = 0.7, double fmax = 2.0,
|
||||
double fsamp=0);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
|
||||
protected:
|
||||
// configuration
|
||||
int _order;
|
||||
double _fmin, _fmax, _fsamp;
|
||||
};
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
class ButterworthBandstop : public BiquadCascade<TYPE> {
|
||||
public:
|
||||
ButterworthBandstop(int order = 3, double fmin = 0.7, double fmax = 2.0,
|
||||
double fsamp=0);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
protected:
|
||||
// configuration
|
||||
int _order;
|
||||
double _fmin, _fmax, _fsamp;
|
||||
};
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
class ButterworthHighLowpass : public BiquadCascade<TYPE> {
|
||||
public:
|
||||
ButterworthHighLowpass(int order = 3, double fmin = 0.7,
|
||||
double fmax = 2.0, double fsamp = 0);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
private:
|
||||
int _order;
|
||||
double _fmin, _fmax, _fsamp;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering::IIR
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
264
include/seiscomp/math/filter/butterworth.ipp
Normal file
264
include/seiscomp/math/filter/butterworth.ipp
Normal file
@ -0,0 +1,264 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
// This file is included by "butterworth.cpp"
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
ButterworthLowpass<TYPE>::ButterworthLowpass(int order, double fmax, double fsamp)
|
||||
: _order(order), _fmax(fmax), _fsamp(0.0) {
|
||||
if ( fsamp ) {
|
||||
setSamplingFrequency(fsamp);
|
||||
}
|
||||
}
|
||||
|
||||
template<class TYPE>
|
||||
void ButterworthLowpass<TYPE>::setSamplingFrequency(double fsamp) {
|
||||
if (_fsamp == fsamp) {
|
||||
return;
|
||||
}
|
||||
|
||||
double fmax = _fmax;
|
||||
|
||||
if ( fmax < 0 ) {
|
||||
fmax = -fmax * fsamp;
|
||||
}
|
||||
|
||||
BiquadCascade<TYPE>::set(init_bw_biquads(_order, 0, fmax, fsamp, BUTTERWORTH_LOWPASS));
|
||||
_fsamp = fsamp;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
int ButterworthLowpass<TYPE>::setParameters(int n, const double *params) {
|
||||
if ( n != 2 ) {
|
||||
return 2;
|
||||
}
|
||||
|
||||
int order = (int)params[0];
|
||||
|
||||
if ( order <= 0 ) {
|
||||
return -1;
|
||||
}
|
||||
|
||||
_order = order;
|
||||
_fmax = params[1];
|
||||
|
||||
return 2;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
InPlaceFilter<TYPE>* ButterworthLowpass<TYPE>::clone() const {
|
||||
return new ButterworthLowpass<TYPE>(_order, _fmax, _fsamp);
|
||||
}
|
||||
|
||||
template<class TYPE>
|
||||
ButterworthHighpass<TYPE>::ButterworthHighpass(int order, double fmin, double fsamp)
|
||||
: _order(order), _fmin(fmin), _fsamp(0.0) {
|
||||
if ( fsamp )
|
||||
setSamplingFrequency(fsamp);
|
||||
}
|
||||
|
||||
template<class TYPE>
|
||||
void ButterworthHighpass<TYPE>::setSamplingFrequency(double fsamp) {
|
||||
if ( _fsamp == fsamp ) {
|
||||
return;
|
||||
}
|
||||
|
||||
double fmin = _fmin;
|
||||
|
||||
if ( fmin < 0 ) {
|
||||
fmin = -fmin * fsamp;
|
||||
}
|
||||
|
||||
BiquadCascade<TYPE>::set(init_bw_biquads(_order, fmin, 0, fsamp, BUTTERWORTH_HIGHPASS));
|
||||
_fsamp = fsamp;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
int ButterworthHighpass<TYPE>::setParameters(int n, const double *params) {
|
||||
if ( n != 2 ) {
|
||||
return 2;
|
||||
}
|
||||
|
||||
int order = (int)params[0];
|
||||
|
||||
if ( order <= 0 ) {
|
||||
return -1;
|
||||
}
|
||||
|
||||
_order = order;
|
||||
_fmin = params[1];
|
||||
|
||||
return 2;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
InPlaceFilter<TYPE>* ButterworthHighpass<TYPE>::clone() const {
|
||||
return new ButterworthHighpass<TYPE>(_order, _fmin, _fsamp);
|
||||
}
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
ButterworthBandpass<TYPE>::ButterworthBandpass(int order, double fmin, double fmax,
|
||||
double fsamp)
|
||||
: _order(order), _fmin(fmin), _fmax(fmax), _fsamp(0.0){
|
||||
if ( fsamp ) {
|
||||
setSamplingFrequency(fsamp);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
void ButterworthBandpass<TYPE>::setSamplingFrequency(double fsamp) {
|
||||
if ( _fsamp == fsamp ) {
|
||||
return;
|
||||
}
|
||||
|
||||
double fmin = _fmin;
|
||||
double fmax = _fmax;
|
||||
|
||||
if ( fmin < 0 ) {
|
||||
fmin = -fmin * fsamp;
|
||||
}
|
||||
|
||||
if ( fmax < 0 ) {
|
||||
fmax = -fmax * fsamp;
|
||||
}
|
||||
|
||||
BiquadCascade<TYPE>::set(init_bw_biquads(_order, fmin, fmax, fsamp, BUTTERWORTH_BANDPASS));
|
||||
_fsamp = fsamp;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
int ButterworthBandpass<TYPE>::setParameters(int n, const double *params) {
|
||||
if ( n != 3 ) {
|
||||
return 3;
|
||||
}
|
||||
|
||||
_order = (int)params[0];
|
||||
_fmin = params[1];
|
||||
_fmax = params[2];
|
||||
|
||||
return n;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
InPlaceFilter<TYPE>* ButterworthBandpass<TYPE>::clone() const {
|
||||
return new ButterworthBandpass<TYPE>(_order, _fmin, _fmax, _fsamp);
|
||||
}
|
||||
|
||||
template<class TYPE>
|
||||
ButterworthBandstop<TYPE>::ButterworthBandstop(int order, double fmin, double fmax,
|
||||
double fsamp)
|
||||
: _order(order), _fmin(fmin), _fmax(fmax), _fsamp(0.0) {
|
||||
if ( fsamp ) {
|
||||
setSamplingFrequency(fsamp);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
void ButterworthBandstop<TYPE>::setSamplingFrequency(double fsamp) {
|
||||
if ( _fsamp == fsamp ) {
|
||||
return;
|
||||
}
|
||||
|
||||
double fmin = _fmin;
|
||||
double fmax = _fmax;
|
||||
|
||||
if ( fmin < 0 ) {
|
||||
fmin = -fmin * fsamp;
|
||||
}
|
||||
|
||||
if ( fmax < 0 ) {
|
||||
fmax = -fmax * fsamp;
|
||||
}
|
||||
|
||||
BiquadCascade<TYPE>::set(init_bw_biquads(_order, fmin, fmax, fsamp, BUTTERWORTH_BANDSTOP));
|
||||
_fsamp = fsamp;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
int ButterworthBandstop<TYPE>::setParameters(int n, const double *params) {
|
||||
if ( n != 3 ) {
|
||||
return 3;
|
||||
}
|
||||
|
||||
_order = (int)params[0];
|
||||
_fmin = params[1];
|
||||
_fmax = params[2];
|
||||
|
||||
return n;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
InPlaceFilter<TYPE>* ButterworthBandstop<TYPE>::clone() const {
|
||||
return new ButterworthBandstop<TYPE>(_order, _fmin, _fmax, _fsamp);
|
||||
}
|
||||
|
||||
template<class TYPE>
|
||||
ButterworthHighLowpass<TYPE>::ButterworthHighLowpass(int order, double fmin,
|
||||
double fmax, double fsamp)
|
||||
: _order(order), _fmin(fmin), _fmax(fmax), _fsamp(0.0) {
|
||||
if ( fsamp ) {
|
||||
setSamplingFrequency(fsamp);
|
||||
}
|
||||
}
|
||||
|
||||
template<class TYPE>
|
||||
void ButterworthHighLowpass<TYPE>::setSamplingFrequency(double fsamp) {
|
||||
if ( _fsamp == fsamp ) {
|
||||
return;
|
||||
}
|
||||
|
||||
double fmin = _fmin;
|
||||
double fmax = _fmax;
|
||||
|
||||
if ( fmin < 0 ) {
|
||||
fmin = -fmin * fsamp;
|
||||
}
|
||||
|
||||
if ( fmax < 0 ) {
|
||||
fmax = -fmax * fsamp;
|
||||
}
|
||||
|
||||
BiquadCascade<TYPE>::set(init_bw_biquads(_order, fmin, fmax, fsamp, BUTTERWORTH_HIGHLOWPASS));
|
||||
_fsamp = fsamp;
|
||||
}
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
int ButterworthHighLowpass<TYPE>::setParameters(int n, const double *params) {
|
||||
if ( n != 3 ) {
|
||||
return 3;
|
||||
}
|
||||
|
||||
_order = (int)params[0];
|
||||
_fmin = params[1];
|
||||
_fmax = params[2];
|
||||
|
||||
return n;
|
||||
}
|
||||
|
||||
|
||||
template<class TYPE>
|
||||
InPlaceFilter<TYPE>* ButterworthHighLowpass<TYPE>::clone() const {
|
||||
return new ButterworthHighLowpass<TYPE>(_order, _fmin, _fmax, _fsamp);
|
||||
}
|
||||
102
include/seiscomp/math/filter/chainfilter.h
Normal file
102
include/seiscomp/math/filter/chainfilter.h
Normal file
@ -0,0 +1,102 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_CHAINFILTER
|
||||
#define SEISCOMP_MATH_CHAINFILTER
|
||||
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class ChainFilter : public InPlaceFilter<TYPE> {
|
||||
// ------------------------------------------------------------------
|
||||
// X'truction
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
ChainFilter();
|
||||
~ChainFilter();
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Interface
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
//! Adds a filter to the chain. The ownership of the filter to
|
||||
//! be added goes the ChainFilter instance which will delete
|
||||
//! the filter.
|
||||
bool add(InPlaceFilter<TYPE> *filter);
|
||||
|
||||
//! Removes and deletes the filter at a certain position in
|
||||
//! the chain.
|
||||
bool remove(size_t pos);
|
||||
|
||||
//! Removes the filter at position pos and returns it. The
|
||||
//! filter instance will not be deleted and the ownership goes
|
||||
//! to the caller.
|
||||
InPlaceFilter<TYPE>* take(size_t pos);
|
||||
|
||||
//! Returns the index of a filter in the chain or -1 if not
|
||||
//! found
|
||||
size_t indexOf(InPlaceFilter<TYPE> *filter) const;
|
||||
|
||||
//! Returns the number of filters in the chain
|
||||
size_t filterCount() const;
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Derived filter interface
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
void apply(int n, TYPE *inout) override;
|
||||
|
||||
void setStartTime(const Core::Time &time) override;
|
||||
|
||||
void setStreamID(const std::string &net,
|
||||
const std::string &sta,
|
||||
const std::string &loc,
|
||||
const std::string &cha) override;
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Private members
|
||||
// ------------------------------------------------------------------
|
||||
private:
|
||||
typedef std::vector<InPlaceFilter<TYPE>*> FilterChain;
|
||||
FilterChain _filters;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filter
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
55
include/seiscomp/math/filter/const.h
Normal file
55
include/seiscomp/math/filter/const.h
Normal file
@ -0,0 +1,55 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_CONST_H
|
||||
#define SEISCOMP_MATH_FILTER_CONST_H
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename T>
|
||||
class ConstFilter : public InPlaceFilter<T> {
|
||||
public:
|
||||
ConstFilter(T c = 1);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
void apply(int n, T *inout) override;
|
||||
|
||||
InPlaceFilter<T>* clone() const override;
|
||||
|
||||
private:
|
||||
T _const;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
58
include/seiscomp/math/filter/cutoff.h
Normal file
58
include/seiscomp/math/filter/cutoff.h
Normal file
@ -0,0 +1,58 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_CUTOFF_H
|
||||
#define SEISCOMP_MATH_FILTER_CUTOFF_H
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class CutOff : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
CutOff(TYPE threshold = 0);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
private:
|
||||
TYPE _threshold;
|
||||
TYPE _samples[2];
|
||||
int _outstanding;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
59
include/seiscomp/math/filter/duration.h
Normal file
59
include/seiscomp/math/filter/duration.h
Normal file
@ -0,0 +1,59 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_DURATION_H
|
||||
#define SEISCOMP_MATH_FILTER_DURATION_H
|
||||
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class Duration : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
Duration(TYPE triggerOn = 3, TYPE triggerOff = 1.5);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
private:
|
||||
TYPE _triggerOn;
|
||||
TYPE _triggerOff;
|
||||
double _dt;
|
||||
int _triggerCount{0};
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
66
include/seiscomp/math/filter/iirdifferentiate.h
Normal file
66
include/seiscomp/math/filter/iirdifferentiate.h
Normal file
@ -0,0 +1,66 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_IIRDIFFERENTIATE_H
|
||||
#define SEISCOMP_IIRDIFFERENTIATE_H
|
||||
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template <typename T>
|
||||
class IIRDifferentiate : public InPlaceFilter<T> {
|
||||
public:
|
||||
IIRDifferentiate(double fsamp = 0);
|
||||
IIRDifferentiate(const IIRDifferentiate<T> &other);
|
||||
|
||||
|
||||
public:
|
||||
void reset();
|
||||
|
||||
|
||||
// InPlaceFilter interface
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
void apply(int n, T *inout) override;
|
||||
|
||||
InPlaceFilter<T> *clone() const override;
|
||||
|
||||
private:
|
||||
T _v1;
|
||||
T _fsamp;
|
||||
bool _init;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
70
include/seiscomp/math/filter/iirintegrate.h
Normal file
70
include/seiscomp/math/filter/iirintegrate.h
Normal file
@ -0,0 +1,70 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_IIRINTEGRATE_H
|
||||
#define SEISCOMP_IIRINTEGRATE_H
|
||||
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template <typename T>
|
||||
class IIRIntegrate : public InPlaceFilter<T> {
|
||||
public:
|
||||
IIRIntegrate(double a = 0, double fsamp = 0);
|
||||
IIRIntegrate(const IIRIntegrate<T> &other);
|
||||
|
||||
|
||||
public:
|
||||
void reset();
|
||||
|
||||
|
||||
// InPlaceFilter interface
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
void apply(int n, T *inout) override;
|
||||
|
||||
InPlaceFilter<T> *clone() const override;
|
||||
|
||||
private:
|
||||
void init(double a);
|
||||
|
||||
private:
|
||||
double _ia0, _ia1, _ia2;
|
||||
|
||||
double _a0, _a1, _a2;
|
||||
double _b0, _b1, _b2;
|
||||
T _v1, _v2;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
67
include/seiscomp/math/filter/median.h
Normal file
67
include/seiscomp/math/filter/median.h
Normal file
@ -0,0 +1,67 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_MEDIAN_H
|
||||
#define SEISCOMP_MATH_FILTER_MEDIAN_H
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class Median : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
Median(double timeSpan /*sec*/ = 1.0, double fsamp = 0.0);
|
||||
|
||||
public:
|
||||
void setLength(double timeSpan);
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
// Resets the filter values
|
||||
void reset();
|
||||
|
||||
private:
|
||||
double _timeSpan;
|
||||
double _fsamp;
|
||||
size_t _sampleCount;
|
||||
size_t _index;
|
||||
bool _firstSample;
|
||||
std::vector<TYPE> _buffer; // ring buffer of last _sampleCount samples
|
||||
std::vector<TYPE> _sorted; // sorted version of ring buffer
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
|
||||
105
include/seiscomp/math/filter/minmax.h
Normal file
105
include/seiscomp/math/filter/minmax.h
Normal file
@ -0,0 +1,105 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_MINMAX_H
|
||||
#define SEISCOMP_MATH_FILTER_MINMAX_H
|
||||
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
/**
|
||||
* @brief The MinMax class is the base class for either the Min or Max
|
||||
* filter. Each output sample holds the minimum/maximum of all prior
|
||||
* samples within the configured time window.
|
||||
*/
|
||||
template<typename TYPE>
|
||||
class MinMax : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
MinMax(double timeSpan /*sec*/ = 1.0, double fsamp = 0.0);
|
||||
|
||||
|
||||
public:
|
||||
/**
|
||||
* @brief Sets the length of the minmax time window.
|
||||
* @param timeSpan Length in seconds
|
||||
*/
|
||||
void setLength(double timeSpan);
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// Resets the filter values
|
||||
void reset();
|
||||
|
||||
|
||||
protected:
|
||||
double _timeSpan;
|
||||
double _fsamp;
|
||||
int _sampleCount;
|
||||
int _index;
|
||||
bool _firstSample;
|
||||
std::vector<TYPE> _buffer;
|
||||
};
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class Min : public MinMax<TYPE> {
|
||||
public:
|
||||
Min(double timeSpan /*sec*/ = 1.0, double fsamp = 0.0);
|
||||
|
||||
public:
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
private:
|
||||
TYPE _minimum;
|
||||
};
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class Max : public MinMax<TYPE> {
|
||||
public:
|
||||
Max(double timeSpan /*sec*/ = 1.0, double fsamp = 0.0);
|
||||
|
||||
public:
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
private:
|
||||
TYPE _maximum;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
77
include/seiscomp/math/filter/op2filter.h
Normal file
77
include/seiscomp/math/filter/op2filter.h
Normal file
@ -0,0 +1,77 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_OP2FILTER
|
||||
#define SEISCOMP_MATH_OP2FILTER
|
||||
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
class Op2Filter : public InPlaceFilter<TYPE> {
|
||||
// ------------------------------------------------------------------
|
||||
// X'truction
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
Op2Filter(InPlaceFilter<TYPE> *op1, InPlaceFilter<TYPE> *op2, double fsamp=0.0);
|
||||
~Op2Filter();
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Derived filter interface
|
||||
// ------------------------------------------------------------------
|
||||
public:
|
||||
void apply(int n, TYPE *inout) override;
|
||||
|
||||
void setStartTime(const Core::Time &time) override;
|
||||
void setStreamID(const std::string &net,
|
||||
const std::string &sta,
|
||||
const std::string &loc,
|
||||
const std::string &cha) override;
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, double const *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
|
||||
// ------------------------------------------------------------------
|
||||
// Private members
|
||||
// ------------------------------------------------------------------
|
||||
private:
|
||||
InPlaceFilter<TYPE> *_op1;
|
||||
InPlaceFilter<TYPE> *_op2;
|
||||
};
|
||||
|
||||
|
||||
#include <seiscomp/math/filter/op2filter.ipp>
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filter
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
88
include/seiscomp/math/filter/op2filter.ipp
Normal file
88
include/seiscomp/math/filter/op2filter.ipp
Normal file
@ -0,0 +1,88 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
Op2Filter<TYPE, OPERATION>::Op2Filter(InPlaceFilter<TYPE> *op1, InPlaceFilter<TYPE> *op2, double fsamp) {
|
||||
_op1 = op1;
|
||||
_op2 = op2;
|
||||
|
||||
if ( fsamp )
|
||||
setSamplingFrequency(fsamp);
|
||||
}
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
Op2Filter<TYPE, OPERATION>::~Op2Filter() {
|
||||
if ( _op1 ) delete _op1;
|
||||
if ( _op2 ) delete _op2;
|
||||
}
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
void Op2Filter<TYPE, OPERATION>::apply(int n, TYPE *inout) {
|
||||
OPERATION<TYPE> op;
|
||||
for ( int i = 0; i < n; ++i ) {
|
||||
TYPE v1 = *inout, v2 = *inout;
|
||||
_op1->apply(1, &v1);
|
||||
_op2->apply(1, &v2);
|
||||
|
||||
*inout = op(v1, v2);
|
||||
|
||||
++inout;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
InPlaceFilter<TYPE>* Op2Filter<TYPE, OPERATION>::clone() const {
|
||||
return new Op2Filter<TYPE, OPERATION>(_op1->clone(), _op2->clone());
|
||||
}
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
void Op2Filter<TYPE, OPERATION>::setStartTime(const Core::Time &time) {
|
||||
if ( _op1 ) _op1->setStartTime(time);
|
||||
if ( _op2 ) _op2->setStartTime(time);
|
||||
}
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
void Op2Filter<TYPE, OPERATION>::setStreamID(const std::string &net,
|
||||
const std::string &sta,
|
||||
const std::string &loc,
|
||||
const std::string &cha) {
|
||||
if ( _op1 ) _op1->setStreamID(net, sta, loc, cha);
|
||||
if ( _op2 ) _op2->setStreamID(net, sta, loc, cha);
|
||||
}
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
void Op2Filter<TYPE, OPERATION>::setSamplingFrequency(double fsamp) {
|
||||
if ( _op1 ) _op1->setSamplingFrequency(fsamp);
|
||||
if ( _op2 ) _op2->setSamplingFrequency(fsamp);
|
||||
}
|
||||
|
||||
|
||||
template<typename TYPE, template <class T> class OPERATION>
|
||||
int Op2Filter<TYPE, OPERATION>::setParameters(int n, double const *params) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
84
include/seiscomp/math/filter/random.h
Normal file
84
include/seiscomp/math/filter/random.h
Normal file
@ -0,0 +1,84 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_RANDOM_H
|
||||
#define SEISCOMP_MATH_FILTER_RANDOM_H
|
||||
|
||||
#include <random>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class RandomUniform : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
RandomUniform(double minValue /*sec*/ = -1.0, double maxValue /*sec*/ = 1.0);
|
||||
|
||||
public:
|
||||
void setLength(double timeSpan);
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
// Resets the filter values
|
||||
void reset();
|
||||
|
||||
private:
|
||||
double _minimumValue{-1.0};
|
||||
double _maximumValue{1.0};
|
||||
};
|
||||
|
||||
template<typename TYPE>
|
||||
class RandomNormal : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
RandomNormal(double meanValue /*sec*/ = 0.0, double stdDev /*sec*/ = 1.0);
|
||||
|
||||
public:
|
||||
void setLength(double timeSpan);
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
// Resets the filter values
|
||||
void reset();
|
||||
|
||||
private:
|
||||
double _meanValue{0.0};
|
||||
double _standardDev{1.0};
|
||||
};
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
|
||||
84
include/seiscomp/math/filter/rmhp.h
Normal file
84
include/seiscomp/math/filter/rmhp.h
Normal file
@ -0,0 +1,84 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_FILTERING_RMHP_H
|
||||
#define SEISCOMP_FILTERING_RMHP_H
|
||||
|
||||
|
||||
#include<vector>
|
||||
#include<seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class RunningMean : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
RunningMean(double windowLength=0, double fsamp=0.0);
|
||||
~RunningMean() {}
|
||||
|
||||
|
||||
public:
|
||||
void setLength(double windowLength) {
|
||||
_windowLength = windowLength;
|
||||
}
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// resets the filter, i.e. erases the filter memory
|
||||
void reset();
|
||||
|
||||
|
||||
protected:
|
||||
double _windowLength, _samplingFrequency;
|
||||
int _windowLengthI, _sampleCount;
|
||||
double _average;
|
||||
};
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class RunningMeanHighPass : public RunningMean<TYPE> {
|
||||
public:
|
||||
RunningMeanHighPass(double windowLength=0, double fsamp=0.0);
|
||||
~RunningMeanHighPass() {}
|
||||
|
||||
|
||||
public:
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
218
include/seiscomp/math/filter/seismometers.h
Normal file
218
include/seiscomp/math/filter/seismometers.h
Normal file
@ -0,0 +1,218 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_SEISMOMETERS_H
|
||||
#define SEISCOMP_SEISMOMETERS_H
|
||||
|
||||
|
||||
#include <complex>
|
||||
#include <seiscomp/math/filter/biquad.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
enum GroundMotion { Displacement, Velocity, Acceleration };
|
||||
|
||||
namespace SeismometerResponse {
|
||||
|
||||
typedef std::complex<double> Pole;
|
||||
typedef std::complex<double> Zero;
|
||||
|
||||
struct FAP {
|
||||
FAP() {}
|
||||
FAP(double f, double a, double p)
|
||||
: frequency(f), amplitude(a), phaseAngle(p) {}
|
||||
|
||||
bool operator<(const FAP &other) const {
|
||||
return frequency < other.frequency;
|
||||
}
|
||||
|
||||
double frequency; //! Frequency in Hz
|
||||
double amplitude; //! Amplitude
|
||||
double phaseAngle; //! Phase angle in degree
|
||||
};
|
||||
|
||||
typedef std::vector<Pole> Poles;
|
||||
typedef std::vector<Zero> Zeros;
|
||||
typedef std::vector<FAP> FAPs;
|
||||
|
||||
class PolesAndZeros {
|
||||
public:
|
||||
PolesAndZeros();
|
||||
|
||||
PolesAndZeros(const Poles &poles, const Zeros &zeros,
|
||||
double norm);
|
||||
|
||||
PolesAndZeros(const PolesAndZeros &other);
|
||||
|
||||
public:
|
||||
Poles poles;
|
||||
Zeros zeros;
|
||||
double norm;
|
||||
};
|
||||
|
||||
class WoodAnderson : public PolesAndZeros {
|
||||
public:
|
||||
// Gutenberg(1935) -> gain=2800, T0=0.8s, h=0.8
|
||||
// Uhrhammer and Collins (1990) -> gain=2080, T0=0.8s, h=0.7
|
||||
//
|
||||
// Note that use of the Uhrhammer and Collins (1990) version
|
||||
// is recommended by the IASPEI Magnitude Working Group
|
||||
// recommendations of 2011 September 9.
|
||||
struct Config {
|
||||
Config(double gain=2080, double T0=0.8, double h=0.7) :
|
||||
gain(gain), T0(T0), h(h) {}
|
||||
|
||||
double gain, T0, h;
|
||||
};
|
||||
WoodAnderson(GroundMotion input, Config config=Config());
|
||||
};
|
||||
|
||||
|
||||
class Seismometer5sec : public PolesAndZeros {
|
||||
public:
|
||||
Seismometer5sec(GroundMotion input);
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
// Responses that soon need to be implemented are (at least)
|
||||
// * STS2
|
||||
// * STS1
|
||||
// * Kirnos
|
||||
// * corresponding restitution filters
|
||||
// * generic filters that read PAZ from file/database/etc.
|
||||
|
||||
namespace IIR {
|
||||
|
||||
template <typename T>
|
||||
class Filter : public SeismometerResponse::PolesAndZeros,
|
||||
public Math::Filtering::InPlaceFilter<T> {
|
||||
public:
|
||||
Filter();
|
||||
|
||||
Filter(const SeismometerResponse::Poles &poles,
|
||||
const SeismometerResponse::Zeros &zeros, double norm);
|
||||
|
||||
Filter(const Filter &other);
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
void apply(int n, T *inout) override;
|
||||
|
||||
Math::Filtering::InPlaceFilter<T> *clone() const override;
|
||||
|
||||
private:
|
||||
Math::Filtering::IIR::BiquadCascade<T> _cascade;
|
||||
};
|
||||
|
||||
|
||||
template <typename T>
|
||||
class WWSSN_SP_Filter : public Filter<T> {
|
||||
public:
|
||||
WWSSN_SP_Filter(GroundMotion input=Velocity);
|
||||
WWSSN_SP_Filter(const WWSSN_SP_Filter &other);
|
||||
|
||||
public:
|
||||
int setParameters(int n, const double *params) override;
|
||||
Math::Filtering::InPlaceFilter<T> *clone() const override;
|
||||
|
||||
void setInput(GroundMotion input);
|
||||
};
|
||||
|
||||
|
||||
template <typename T>
|
||||
class WWSSN_LP_Filter : public Filter<T> {
|
||||
public:
|
||||
WWSSN_LP_Filter(GroundMotion input=Velocity);
|
||||
WWSSN_LP_Filter(const WWSSN_LP_Filter &other);
|
||||
|
||||
public:
|
||||
int setParameters(int n, const double *params) override;
|
||||
Math::Filtering::InPlaceFilter<T> *clone() const override;
|
||||
|
||||
void setInput(GroundMotion input);
|
||||
};
|
||||
|
||||
|
||||
template <typename T>
|
||||
class WoodAndersonFilter : public Filter<T> {
|
||||
public:
|
||||
WoodAndersonFilter(GroundMotion input=Velocity, SeismometerResponse::WoodAnderson::Config config = SeismometerResponse::WoodAnderson::Config());
|
||||
WoodAndersonFilter(const WoodAndersonFilter &other);
|
||||
|
||||
public:
|
||||
int setParameters(int n, const double *params) override;
|
||||
Math::Filtering::InPlaceFilter<T> *clone() const override;
|
||||
|
||||
void setInput(GroundMotion input,
|
||||
SeismometerResponse::WoodAnderson::Config config = SeismometerResponse::WoodAnderson::Config());
|
||||
};
|
||||
|
||||
|
||||
template <typename T>
|
||||
class GenericSeismometer : public Filter<T> {
|
||||
public:
|
||||
GenericSeismometer(double cornerPeriod=1.,
|
||||
GroundMotion input=Velocity);
|
||||
GenericSeismometer(const GenericSeismometer &other);
|
||||
|
||||
public:
|
||||
int setParameters(int n, const double *params) override;
|
||||
Math::Filtering::InPlaceFilter<T> *clone() const override;
|
||||
|
||||
void setInput(GroundMotion input);
|
||||
|
||||
private:
|
||||
double _cornerPeriod;
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
class Seismometer5secFilter : public Filter<T> {
|
||||
public:
|
||||
Seismometer5secFilter(GroundMotion input=Velocity);
|
||||
Seismometer5secFilter(const Seismometer5secFilter &other);
|
||||
|
||||
public:
|
||||
int setParameters(int n, const double *params) override;
|
||||
Math::Filtering::InPlaceFilter<T> *clone() const override;
|
||||
|
||||
void setInput(GroundMotion input);
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
74
include/seiscomp/math/filter/sr.h
Normal file
74
include/seiscomp/math/filter/sr.h
Normal file
@ -0,0 +1,74 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_SR_H
|
||||
#define SEISCOMP_MATH_FILTER_SR_H
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename T>
|
||||
class SamplingRate : public InPlaceFilter<T> {
|
||||
public:
|
||||
SamplingRate() = default;
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
void apply(int n, T *inout) override;
|
||||
|
||||
InPlaceFilter<T> *clone() const override;
|
||||
|
||||
private:
|
||||
T _fsamp{0};
|
||||
};
|
||||
|
||||
|
||||
template<typename T>
|
||||
class SamplingTime : public InPlaceFilter<T> {
|
||||
public:
|
||||
SamplingTime() = default;
|
||||
|
||||
public:
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
void apply(int n, T *inout) override;
|
||||
|
||||
InPlaceFilter<T> *clone() const override;
|
||||
|
||||
private:
|
||||
T _dt{0};
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
|
||||
186
include/seiscomp/math/filter/stalta.h
Normal file
186
include/seiscomp/math/filter/stalta.h
Normal file
@ -0,0 +1,186 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_PICKER_STALTA_H
|
||||
#define SEISCOMP_PICKER_STALTA_H
|
||||
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
// Recursive STA/LTA filter
|
||||
//
|
||||
// Very simple and doesn't need a lot of memory for buffering.
|
||||
|
||||
template<typename TYPE>
|
||||
class STALTA : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
STALTA(double lenSTA=2, double lenLTA=50, double fsamp=1.);
|
||||
|
||||
public:
|
||||
// Apply the picker in place to the (previously filtered) data.
|
||||
void apply(int ndata, TYPE *data) override;
|
||||
|
||||
void reset();
|
||||
|
||||
// Set the sampling frequency in Hz. Allows delayed
|
||||
// initialization when the data arrive.
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
protected:
|
||||
// length of STA and LTA windows in seconds
|
||||
double _lenSTA, _lenLTA;
|
||||
|
||||
// length of STA and LTA windows in samples
|
||||
int _numSTA, _numLTA;
|
||||
|
||||
// number of samples processed since init/reset
|
||||
int _sampleCount;
|
||||
|
||||
// The *current* and continuously updated STA and LTA values
|
||||
double _STA, _LTA; // must be double
|
||||
|
||||
// sampling frequency in Hz
|
||||
double _fsamp;
|
||||
};
|
||||
|
||||
|
||||
|
||||
|
||||
// Another recursive STA/LTA filter
|
||||
//
|
||||
// In this version the LTA is not updated during an event,
|
||||
// making it more sensitive to secondary signals.
|
||||
|
||||
template<typename TYPE>
|
||||
class STALTA2 : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
STALTA2(double lenSTA=2, double lenLTA=50,
|
||||
double eventON=3., double eventOFF=1.,
|
||||
double fsamp=1.);
|
||||
|
||||
public:
|
||||
// Apply the picker in place to the (previously filtered) data.
|
||||
void apply(int ndata, TYPE *data) override;
|
||||
|
||||
void reset();
|
||||
|
||||
// Set the sampling frequency in Hz. Allows delayed
|
||||
// initialization when the data arrive.
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
protected:
|
||||
// length of STA and LTA windows in seconds
|
||||
double _lenSTA, _lenLTA;
|
||||
|
||||
// length of STA and LTA windows in samples
|
||||
int _numSTA, _numLTA;
|
||||
|
||||
// number of samples processed since init/reset
|
||||
int _sampleCount;
|
||||
|
||||
// The *current* and continuously updated STA and LTA values
|
||||
double _STA, _LTA; // must be double
|
||||
|
||||
// sampling frequency in Hz
|
||||
double _fsamp;
|
||||
|
||||
// thresholds to declare an event on and off
|
||||
double _eventOn, _eventOff;
|
||||
|
||||
// This value controls how much the LTA will be updated during
|
||||
// an event. Current values are only 0 and 1. With 0 meaning
|
||||
// completely frozen LTA during the event, 1 meaning the LTA
|
||||
// is fully updated.
|
||||
double _updateLTA;
|
||||
};
|
||||
|
||||
|
||||
// Classical, non-recursive STA/LTA filter
|
||||
//
|
||||
// As simple as it can get, but also needs a lot more memory.
|
||||
|
||||
template<typename TYPE>
|
||||
class STALTA_Classic : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
STALTA_Classic(
|
||||
double lenSTA = 2.,
|
||||
double lenLTA = 50.,
|
||||
double fsamp = 1.);
|
||||
|
||||
public:
|
||||
// Apply the picker in place to the (previously filtered) data.
|
||||
void apply(int ndata, TYPE *data) override;
|
||||
|
||||
void reset();
|
||||
|
||||
// Set the sampling frequency in Hz. Allows delayed
|
||||
// initialization when the data arrive.
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
protected:
|
||||
// length of STA and LTA windows in seconds
|
||||
double _lenSTA, _lenLTA;
|
||||
|
||||
// length of STA and LTA windows in samples
|
||||
int _numSTA, _numLTA;
|
||||
|
||||
// number of samples processed since init/reset
|
||||
int _sampleCount;
|
||||
|
||||
// The *current* and continuously updated STA and LTA values
|
||||
double _STA, _LTA; // must be double
|
||||
|
||||
// For this non-recursive filter we need to save all samples
|
||||
// that contribute to the STA and LTA. Note that the space
|
||||
// requirement is not optimized because the STA buffer is
|
||||
// redundant. This is for simplicity and because the STA
|
||||
// buffer length is small compared to the LTA buffer length.
|
||||
std::vector<double> _sta_buffer, _lta_buffer;
|
||||
|
||||
// sampling frequency in Hz
|
||||
double _fsamp;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filter
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
67
include/seiscomp/math/filter/sum.h
Normal file
67
include/seiscomp/math/filter/sum.h
Normal file
@ -0,0 +1,67 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_FILTER_SUM_H
|
||||
#define SEISCOMP_MATH_FILTER_SUM_H
|
||||
|
||||
#include <vector>
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class Sum : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
Sum(double timeSpan /*sec*/ = 1.0, double fsamp = 0.0);
|
||||
|
||||
public:
|
||||
void setLength(double timeSpan);
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
// Resets the filter values
|
||||
void reset();
|
||||
|
||||
private:
|
||||
double _timeSpan;
|
||||
double _fsamp;
|
||||
int _sampleCount;
|
||||
int _index;
|
||||
double _lastSum;
|
||||
bool _firstSample;
|
||||
std::vector<TYPE> _buffer;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
69
include/seiscomp/math/filter/taper.h
Normal file
69
include/seiscomp/math/filter/taper.h
Normal file
@ -0,0 +1,69 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_FILTERING_TAPER_H
|
||||
#define SEISCOMP_FILTERING_TAPER_H
|
||||
|
||||
|
||||
#include <seiscomp/math/filter.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Filtering {
|
||||
|
||||
|
||||
template<typename TYPE>
|
||||
class InitialTaper : public InPlaceFilter<TYPE> {
|
||||
public:
|
||||
InitialTaper(double taperLength=0, TYPE offset=0, double fsamp=0);
|
||||
~InitialTaper() {}
|
||||
|
||||
void setLength(double taperLength, TYPE offset=0) {
|
||||
_taperLength = taperLength;
|
||||
_offset = offset;
|
||||
}
|
||||
|
||||
// apply filter to data vector **in*place**
|
||||
void apply(int n, TYPE *inout) override;
|
||||
|
||||
InPlaceFilter<TYPE> *clone() const override;
|
||||
|
||||
void setSamplingFrequency(double fsamp) override;
|
||||
|
||||
int setParameters(int n, const double *params) override;
|
||||
|
||||
// resets the filter, i.e. erases the filter memory
|
||||
void reset() { _sampleCount = 0; }
|
||||
|
||||
|
||||
private:
|
||||
double _taperLength, _samplingFrequency;
|
||||
int _taperLengthI, _sampleCount;
|
||||
TYPE _offset;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Filtering
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
311
include/seiscomp/math/geo.h
Normal file
311
include/seiscomp/math/geo.h
Normal file
@ -0,0 +1,311 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_GEO
|
||||
#define SEISCOMP_MATH_GEO
|
||||
|
||||
#include <seiscomp/core.h>
|
||||
#include <seiscomp/math/coord.h>
|
||||
#include <vector>
|
||||
|
||||
namespace Seiscomp
|
||||
{
|
||||
|
||||
namespace Math
|
||||
{
|
||||
|
||||
namespace Geo
|
||||
{
|
||||
|
||||
/**
|
||||
* For two points (lat1, lon1) and (lat2, lon2),
|
||||
* the angular distance in degrees is returned.
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
double delta(double lat1, double lon1, double lat2, double lon2);
|
||||
|
||||
/**
|
||||
* For two points (lat1, lon1) and (lat2, lon2),
|
||||
* the angular distance 'dist' in degrees,
|
||||
* the azimuth 'azi1' (azimuth of point 2 seen from point 1) and
|
||||
* the azimuth 'azi2' (azimuth of point 1 seen from point 2)
|
||||
* are computed.
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
void delazi(double lat1, double lon1, double lat2, double lon2,
|
||||
double *out_dist = nullptr, double *out_azi1 = nullptr,
|
||||
double *out_azi2 = nullptr);
|
||||
|
||||
|
||||
/**
|
||||
* For two points (lat1, lon1) and (lat2, lon2),
|
||||
* the angular distance 'dist' in degrees,
|
||||
* the azimuth 'azi1' (azimuth of point 2 seen from point 1) and
|
||||
* the azimuth 'azi2' (azimuth of point 1 seen from point 2)
|
||||
* are computed.
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
void delazi_wgs84(double lat1, double lon1, double lat2, double lon2,
|
||||
double *out_dist, double *out_azi1, double *out_azi2);
|
||||
|
||||
|
||||
/**
|
||||
* Computes the coordinates (lat, lon) of the point which
|
||||
* is at an azimuth of 'azi' and a distance of 'dist' as seen
|
||||
* from the point (lat0, lon0).
|
||||
* -> lat, lon
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
void delandaz2coord(double dist, double azi, double lat0, double lon0,
|
||||
double *out_lat, double *out_lon);
|
||||
|
||||
/**
|
||||
* Compute the intersection points between two small circles.
|
||||
* Returns the number of intersection points found; 0 means circles don't
|
||||
* intersect at all (or their centers coincide), 1 means both intersection
|
||||
* points coincide, i.e. the circles just "touch" at one point, 2 is a
|
||||
* normal intersection
|
||||
*
|
||||
* @returns lat, lon of intersection points my the way of latx1, lonx1, latx2,
|
||||
* lonx2
|
||||
*
|
||||
* XXX Corrently the case 1 (circles "touch") is not implemented well. Best
|
||||
* way is to test if in case 2 the intersection points are very close.
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
int scxsc(double lat1, double lon1, double r1,
|
||||
double lat2, double lon2, double r2,
|
||||
double *latx1, double *lonx1, double *latx2, double *lonx2,
|
||||
double epsilon=0);
|
||||
|
||||
/**
|
||||
* "Draws" a small circle
|
||||
*
|
||||
* Around a center point (lat0, lon0), compute n equally spaced points on a
|
||||
* small circle with specified radius (in degrees).
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
int scdraw(double lat0, double lon0, double radius,
|
||||
int n, double *lat, double *lon);
|
||||
|
||||
/**
|
||||
* Definition of degree to km conversion constant is based on the
|
||||
* Mean Radius of the Three Semi-Axes as defined by WGS84
|
||||
* (https://nsgreg.nga.mil/doc/view?i=4085)
|
||||
*
|
||||
* - Semi-major Axis (Equatorial Radius of the Earth) [a] : 6378137.0m
|
||||
* - Flattening Factor of the Earth [1/f] : 298.257223563
|
||||
* - Semi-major Axis (Polar Radius of the Earth) [b] : 6356752.314245179m
|
||||
* b = a * (1 - f)
|
||||
* - Mean Radius of the Three Semi-Axes [r_1] : 6371008.77141506m
|
||||
* r_1 = a * (1 - f/3) = (2a + b) / 3
|
||||
* - Kilometer per Degree [km_deg] : 111.195079734632km
|
||||
* km_deg = r_1 * Pi / (180 * 1000)
|
||||
*
|
||||
* Note:
|
||||
* Prior to SeisComP 6 (API 16.0.0) this constant was set to
|
||||
* 111.1329149013519096
|
||||
*/
|
||||
#define KM_OF_DEGREE 111.195079734632
|
||||
|
||||
|
||||
template<typename T>
|
||||
T deg2km(T deg) { return deg * (T)KM_OF_DEGREE; }
|
||||
|
||||
template<typename T>
|
||||
T km2deg(T km) { return km / (T)KM_OF_DEGREE; }
|
||||
|
||||
/**
|
||||
* Converts X, Y, Z coordinates to LTP coordinates using WGS-84 constants for
|
||||
* the ellipsoid.
|
||||
*
|
||||
* For reference, The X, Y, Z coordinate system has origin at the mass center
|
||||
* of the Earth, with Z axis along the spin axis of the Earth (+ at North pole).
|
||||
* The +X axis emerges from the Earth at the equator-prime meridian
|
||||
* intersection. The +Y axis defines a right-hand coordinate system, emerging
|
||||
* from the Earth at +90 degrees longitude on the equator.
|
||||
|
||||
* The Local Tangential Plane (LTP) coordinates are in latitude/longitude/
|
||||
* altitude where North latitude is + and East longitude is +. Altitude is in
|
||||
* meters above the reference ellipsoid (which may be either above or below the
|
||||
* local mean sea level).
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
void xyz2ltp(const double x, const double y, const double z,
|
||||
double *lat, double *lon, double *alt);
|
||||
|
||||
/**
|
||||
* Converts LTP coordinates to X, Y, Z coordinates using WGS-84 constants for
|
||||
* the ellipsoid.
|
||||
*
|
||||
* For reference, The X, Y, Z coordinate system has origin at the mass center
|
||||
* of the Earth, with Z axis along the spin axis of the Earth (+ at North pole).
|
||||
* The +X axis emerges from the Earth at the equator-prime meridian
|
||||
* intersection. The +Y axis defines a right-hand coordinate system, emerging
|
||||
* from the Earth at +90 degrees longitude on the equator.
|
||||
|
||||
* The Local Tangential Plane (LTP) coordinates are in latitude/longitude/
|
||||
* altitude where North latitude is + and East longitude is +. Altitude is in
|
||||
* meters above the reference ellipsoid (which may be either above or below the
|
||||
* local mean sea level).
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
void ltp2xyz(double lat, double lon, double alt,
|
||||
double *x, double *y, double *z);
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Finds the nearest hotspot of an array of hotspots regarding a given
|
||||
* location.
|
||||
*
|
||||
* Returns the distance in degree and the azimuth seen from the location.
|
||||
*/
|
||||
SC_SYSTEM_CORE_API
|
||||
const NamedCoordD*
|
||||
nearestHotspot(double lat, double lon, double maxDist,
|
||||
int nCoords, const NamedCoordD* coordArray,
|
||||
double *dist, double *azi);
|
||||
|
||||
SC_SYSTEM_CORE_API
|
||||
const NamedCoordD*
|
||||
nearestHotspot(double lat, double lon, double maxDist,
|
||||
const std::vector<NamedCoordD>& coords,
|
||||
double *dist, double *azi);
|
||||
|
||||
|
||||
SC_SYSTEM_CORE_API
|
||||
const CityD*
|
||||
nearestCity(double lat, double lon,
|
||||
double maxDist, double minPopulation,
|
||||
int nCities, const CityD* cityArray,
|
||||
double *dist, double *azi);
|
||||
|
||||
SC_SYSTEM_CORE_API
|
||||
const CityD*
|
||||
nearestCity(double lat, double lon,
|
||||
double maxDist, double minPopulation,
|
||||
const std::vector<CityD>& cities,
|
||||
double *dist, double *azi);
|
||||
|
||||
SC_SYSTEM_CORE_API
|
||||
const CityD*
|
||||
largestCity(double lat, double lon,
|
||||
double maxDist,
|
||||
const std::vector<CityD>& cities,
|
||||
double *dist, double *azi);
|
||||
|
||||
|
||||
class SC_SYSTEM_CORE_API PositionInterpolator {
|
||||
public:
|
||||
/**
|
||||
* Constructs an interpolator between lat1,lon1 and lat2,lon2.
|
||||
* Between both positions 'steps-2' new equidistanct positions
|
||||
* are going to be calculated.
|
||||
* @param lat1 Latitude of position 1
|
||||
* @param lon1 Longitude of position 1
|
||||
* @param lat2 Latitude of position 2
|
||||
* @param lon2 Longitude of position 2
|
||||
* @param steps Number of steps including position 1 and position 2
|
||||
*/
|
||||
PositionInterpolator(double lat1, double lon1,
|
||||
double lat2, double lon2,
|
||||
int steps);
|
||||
|
||||
/**
|
||||
* Same as above but the number of steps depends on the
|
||||
* distance between both positions
|
||||
* steps = int(|pos1-pos2|*stepsDistScale)
|
||||
* @param lat1 Latitude of position 1
|
||||
* @param lon1 Longitude of position 1
|
||||
* @param lat2 Latitude of position 2
|
||||
* @param lon2 Longitude of position 2
|
||||
* @param stepsDistScale The distance scale factor to calculate
|
||||
* the number of steps
|
||||
*/
|
||||
PositionInterpolator(double lat1, double lon1,
|
||||
double lat2, double lon2,
|
||||
double stepsDistScale);
|
||||
|
||||
public:
|
||||
PositionInterpolator& operator++();
|
||||
PositionInterpolator operator++(int);
|
||||
|
||||
bool end() const;
|
||||
|
||||
//! Returns the distance between pos1 and pos2
|
||||
//! given to constructor
|
||||
double overallDistance() const;
|
||||
|
||||
//! Returns the current interpolation distance from pos1
|
||||
double distance() const;
|
||||
//! Returns the azimuth of pos2 seen from pos1
|
||||
double azimuth() const;
|
||||
|
||||
//! Returns the current interpolated latitude
|
||||
double latitude() const;
|
||||
//! Returns the current interpolated longitude
|
||||
double longitude() const;
|
||||
|
||||
private:
|
||||
void update();
|
||||
|
||||
private:
|
||||
double _lat1, _lon1;
|
||||
double _lat2, _lon2;
|
||||
double _dist;
|
||||
double _azimuth;
|
||||
double _stepDist;
|
||||
double _currentDist;
|
||||
double _lat, _lon;
|
||||
|
||||
int _nSteps;
|
||||
int _currentStep;
|
||||
};
|
||||
|
||||
|
||||
inline double PositionInterpolator::overallDistance() const {
|
||||
return _dist;
|
||||
}
|
||||
|
||||
inline double PositionInterpolator::distance() const {
|
||||
return _currentDist;
|
||||
}
|
||||
|
||||
inline double PositionInterpolator::azimuth() const {
|
||||
return _azimuth;
|
||||
}
|
||||
|
||||
inline double PositionInterpolator::latitude() const {
|
||||
return _lat;
|
||||
}
|
||||
|
||||
inline double PositionInterpolator::longitude() const {
|
||||
return _lon;
|
||||
}
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Geo
|
||||
|
||||
} // namespace Seiscomp::Math
|
||||
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
109
include/seiscomp/math/hilbert.ipp
Normal file
109
include/seiscomp/math/hilbert.ipp
Normal file
@ -0,0 +1,109 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#include<math.h>
|
||||
|
||||
// Disable hilbert transformation for now
|
||||
#if 0
|
||||
|
||||
#include<fftw3.h>
|
||||
|
||||
namespace Seiscomp
|
||||
{
|
||||
namespace Math
|
||||
{
|
||||
namespace Filtering
|
||||
{
|
||||
|
||||
|
||||
// low-level, real in-place hilbert transform
|
||||
//
|
||||
// The trace length need not be a power of 2
|
||||
template<typename TYPE>
|
||||
void hilbert_transform(std::vector<TYPE> &trace,
|
||||
int direction=1)
|
||||
{
|
||||
TYPE *f = &trace[0];
|
||||
int ndata=trace.size(),
|
||||
nfft = next_power_of_2(ndata),
|
||||
n2 = nfft/2;
|
||||
|
||||
// temporary trace for real FFTW
|
||||
std::vector<double> temp(nfft);
|
||||
double *g = &temp[0];
|
||||
// copy real array and pad with zeros
|
||||
for (int i=0; i<ndata; i++)
|
||||
g[i] = f[i];
|
||||
for (int i=ndata; i<nfft; i++)
|
||||
g[i] = 0.;
|
||||
|
||||
fftw_plan plan_fwd, plan_inv;
|
||||
plan_fwd = fftw_plan_r2r_1d (nfft, g, g, FFTW_R2HC, FFTW_ESTIMATE);
|
||||
plan_inv = fftw_plan_r2r_1d (nfft, g, g, FFTW_HC2R, FFTW_ESTIMATE);
|
||||
|
||||
fftw_execute(plan_fwd); // forward FFT
|
||||
|
||||
// perform actual Hilbert transform in frequency domain
|
||||
if (direction>0)
|
||||
{
|
||||
for (int i=1; i<n2; i++)
|
||||
{
|
||||
double x = g[nfft-i];
|
||||
g[nfft-i] = -g[i];
|
||||
g[i] = x;
|
||||
}
|
||||
}
|
||||
else {
|
||||
for (int i=1; i<n2; i++) {
|
||||
double x = g[nfft-i];
|
||||
g[nfft-i] = g[i];
|
||||
g[i] = -x;
|
||||
}
|
||||
}
|
||||
|
||||
fftw_execute(plan_inv); // inverse FFT
|
||||
|
||||
// scale and copy back to real array
|
||||
double q = 1./nfft;
|
||||
for (int i=0; i<ndata; i++)
|
||||
f[i] = q * g[i];
|
||||
|
||||
fftw_destroy_plan(plan_fwd);
|
||||
fftw_destroy_plan(plan_inv);
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
void envelope(std::vector<TYPE> &trace)
|
||||
{
|
||||
std::vector<TYPE> copy(trace.begin(), trace.end());
|
||||
hilbert_transform(copy);
|
||||
int nsamp = trace.size();
|
||||
TYPE *f1 = &trace[0], *f2 = ©[0];
|
||||
for(int i=0; i<nsamp; i++)
|
||||
f1[i] = (TYPE) ::sqrt( f1[i]*f1[i] + f2[i]*f2[i] );
|
||||
}
|
||||
|
||||
} // namespace Seiscomp::Math::Filter
|
||||
|
||||
} // namespace Seiscomp::Math
|
||||
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
80
include/seiscomp/math/math.h
Normal file
80
include/seiscomp/math/math.h
Normal file
@ -0,0 +1,80 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_MATH_H
|
||||
#define SEISCOMP_MATH_MATH_H
|
||||
|
||||
|
||||
#include <seiscomp/core.h>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
|
||||
|
||||
#ifndef deg2rad
|
||||
#define deg2rad(d) (M_PI*(d)/180.0)
|
||||
#endif
|
||||
#ifndef rad2deg
|
||||
#define rad2deg(d) (180.0*(d)/M_PI)
|
||||
#endif
|
||||
|
||||
#ifndef HALF_PI
|
||||
#define HALF_PI (M_PI/2)
|
||||
#endif
|
||||
|
||||
|
||||
#if defined(_MSC_VER)
|
||||
#include <float.h>
|
||||
#endif
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
using Complex = std::complex<double>;
|
||||
using Fraction = std::pair<int, int>;
|
||||
|
||||
|
||||
#if defined(WIN32)
|
||||
template <typename T>
|
||||
inline bool isNaN(T v) { return _isnan(v)!=0; }
|
||||
#elif defined(__SUNPRO_CC) || defined(__sun)
|
||||
template <typename T>
|
||||
inline bool isNaN(T v) { return isnan(v)!=0; }
|
||||
#else
|
||||
template <typename T>
|
||||
inline bool isNaN(T v) { return std::isnan(v)!=0; }
|
||||
#endif
|
||||
|
||||
/**
|
||||
* @brief Converts a double to a fraction containing numerator and
|
||||
* denominator. If the double value exceeds the numeric integer limits
|
||||
* numeric_limits<int>::max()/1 or numeric_limits<int>::min()/1 is returned.
|
||||
* @param value The double number to be converted.
|
||||
* @return Fraction containing numerator and denominator.
|
||||
*/
|
||||
Fraction double2frac(double value);
|
||||
|
||||
|
||||
} // namespace Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
|
||||
#endif
|
||||
88
include/seiscomp/math/matrix3.h
Normal file
88
include/seiscomp/math/matrix3.h
Normal file
@ -0,0 +1,88 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_MATRIX3_H
|
||||
#define SEISCOMP_MATH_MATRIX3_H
|
||||
|
||||
|
||||
#include <seiscomp/math/vector3.h>
|
||||
#include <string.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename T>
|
||||
struct Matrix3 {
|
||||
Matrix3() {}
|
||||
Matrix3(const Matrix3<T> &other) {
|
||||
memcpy(d, other.d, sizeof(d));
|
||||
}
|
||||
|
||||
Matrix3<T> &identity();
|
||||
|
||||
Matrix3<T> &setRow(int row, const Vector3<T> &v);
|
||||
Matrix3<T> &setColumn(int col, const Vector3<T> &v);
|
||||
|
||||
//! Returns a copy of the r-th row
|
||||
Vector3<T> row(int r) const;
|
||||
|
||||
//! Returns a copy of the c-th column
|
||||
Vector3<T> column(int c) const;
|
||||
|
||||
Matrix3<T> &loadRotateX(T theta);
|
||||
Matrix3<T> &loadRotateY(T theta);
|
||||
Matrix3<T> &loadRotateZ(T theta);
|
||||
|
||||
Matrix3<T> &mult(const Matrix3<T> &a, const Matrix3<T> &b);
|
||||
|
||||
Vector3<T> &transform(Vector3<T> &dst, const Vector3<T> &v) const;
|
||||
Vector3<T> &invTransform(Vector3<T> &dst, const Vector3<T> &v) const;
|
||||
|
||||
operator T *() { return (T*)this; }
|
||||
operator const T *() const { return (const T*)this; }
|
||||
|
||||
Vector3<T> operator*(const Vector3<T> &v) const;
|
||||
|
||||
// Coefficients
|
||||
union {
|
||||
T d[3][3];
|
||||
|
||||
struct {
|
||||
T _11, _12, _13;
|
||||
T _21, _22, _23;
|
||||
T _31, _32, _33;
|
||||
} c;
|
||||
};
|
||||
};
|
||||
|
||||
|
||||
typedef Matrix3<float> Matrix3f;
|
||||
typedef Matrix3<double> Matrix3d;
|
||||
|
||||
|
||||
#include <seiscomp/math/matrix3.ipp>
|
||||
|
||||
|
||||
} // namespace Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
49
include/seiscomp/math/matrix3.ipp
Normal file
49
include/seiscomp/math/matrix3.ipp
Normal file
@ -0,0 +1,49 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
#define _mat3_elem(row,column) d[row][column]
|
||||
|
||||
template <typename T>
|
||||
inline Vector3<T> &Matrix3<T>::transform(Vector3<T> &dst, const Vector3<T> &v) const {
|
||||
for ( int r = 0; r < 3; ++r )
|
||||
dst[r] = _mat3_elem(r,0)*v[0] +
|
||||
_mat3_elem(r,1)*v[1] +
|
||||
_mat3_elem(r,2)*v[2];
|
||||
|
||||
return dst;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
inline Vector3<T> &Matrix3<T>::invTransform(Vector3<T> &dst, const Vector3<T> &v) const {
|
||||
for ( int r = 0; r < 3; ++r )
|
||||
dst[r] = _mat3_elem(0,r)*v[0] +
|
||||
_mat3_elem(1,r)*v[1] +
|
||||
_mat3_elem(2,r)*v[2];
|
||||
|
||||
return dst;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
inline Vector3<T> Matrix3<T>::operator*(const Vector3<T> &v) const {
|
||||
Vector3<T> r;
|
||||
transform(r, v);
|
||||
return r;
|
||||
}
|
||||
83
include/seiscomp/math/mean.h
Normal file
83
include/seiscomp/math/mean.h
Normal file
@ -0,0 +1,83 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_MEAN_H
|
||||
#define SEISCOMP_MATH_MEAN_H
|
||||
|
||||
#include<vector>
|
||||
#include <seiscomp/core/typedarray.h>
|
||||
#include <seiscomp/core.h>
|
||||
|
||||
namespace Seiscomp
|
||||
{
|
||||
namespace Math
|
||||
{
|
||||
namespace Statistics
|
||||
{
|
||||
|
||||
SC_SYSTEM_CORE_API double mean(const DoubleArray &);
|
||||
SC_SYSTEM_CORE_API double mean(const std::vector<double> &);
|
||||
SC_SYSTEM_CORE_API double mean(int n, const double *);
|
||||
|
||||
SC_SYSTEM_CORE_API double median(const DoubleArray &);
|
||||
SC_SYSTEM_CORE_API double median(const std::vector<double> &);
|
||||
SC_SYSTEM_CORE_API double median(int n, const double *);
|
||||
|
||||
SC_SYSTEM_CORE_API double fractile(const DoubleArray &, double x);
|
||||
SC_SYSTEM_CORE_API double fractile(const std::vector<double> &, double x);
|
||||
SC_SYSTEM_CORE_API double fractile(int n, const double *, double x);
|
||||
|
||||
SC_SYSTEM_CORE_API bool computeTrimmedMean(int n, const double *f, double percent, double &value, double &stdev, double *weights = nullptr);
|
||||
SC_SYSTEM_CORE_API bool computeTrimmedMean(int n, const double *f, double &value, double &stdev, double *weights = nullptr);
|
||||
SC_SYSTEM_CORE_API bool computeTrimmedMean(const std::vector<double> &v, double percent, double &value, double &stdev, std::vector<double> *weights = nullptr);
|
||||
|
||||
SC_SYSTEM_CORE_API bool computeMedianTrimmedMean(int n, const double *f, double distance, double &value, double &stdev, double *weights = nullptr);
|
||||
SC_SYSTEM_CORE_API bool computeMedianTrimmedMean(const std::vector<double> &v, double distance, double &value, double &stdev, std::vector<double> *weights = nullptr);
|
||||
|
||||
SC_SYSTEM_CORE_API bool computeMean(const std::vector<double> &v, double &value, double &stdev);
|
||||
|
||||
SC_SYSTEM_CORE_API bool average(int n, const double *values, const double *weights, double &value, double &stdev);
|
||||
SC_SYSTEM_CORE_API bool average(const std::vector<double> &values, const std::vector<double> &weights, double &value, double &stdev);
|
||||
|
||||
SC_SYSTEM_CORE_API void computeLinearTrend(const std::vector<float> &data, double &m, double &n);
|
||||
SC_SYSTEM_CORE_API void computeLinearTrend(const std::vector<double> &data, double &m, double &n);
|
||||
SC_SYSTEM_CORE_API void computeLinearTrend(int cnt, const float *data, double &m, double &n);
|
||||
SC_SYSTEM_CORE_API void computeLinearTrend(int cnt, const double *data, double &m, double &n);
|
||||
|
||||
SC_SYSTEM_CORE_API void detrend(std::vector<float> &data, double m, double n);
|
||||
SC_SYSTEM_CORE_API void detrend(std::vector<double> &data, double m, double n);
|
||||
SC_SYSTEM_CORE_API void detrend(int cnt, float *data, double m, double n);
|
||||
SC_SYSTEM_CORE_API void detrend(int cnt, double *data, double m, double n);
|
||||
|
||||
// From Rosenberger, J.L., and Gasko, M. (1983). "Comparing Location
|
||||
// Estimators: Trimmed Means, Medians, and Trimean", in Understanding
|
||||
// Robust and Exploratory Data Analysis, ed. Hoaglin, D.C., Mosteller,
|
||||
// F., and Tukey, J.W., p. 297-338, John Wiley, NY
|
||||
|
||||
// XXX deprecated XXX
|
||||
SC_SYSTEM_CORE_API double trimmedMean(const DoubleArray &, double percent=25);
|
||||
SC_SYSTEM_CORE_API double trimmedMean(const std::vector<double> &f, double percent=25);
|
||||
SC_SYSTEM_CORE_API double trimmedMean(int n, const double *f, double percent=25);
|
||||
|
||||
} // namespace Statistics
|
||||
} // namespace Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
133
include/seiscomp/math/minmax.ipp
Normal file
133
include/seiscomp/math/minmax.ipp
Normal file
@ -0,0 +1,133 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
template <class TYPE>
|
||||
int
|
||||
minmax(int n, const TYPE *f, int i1, int i2, TYPE *fmin, TYPE *fmax)
|
||||
{
|
||||
int i;
|
||||
/*
|
||||
if(i2<i1) {
|
||||
i=i1; i1=i2; i2=i;
|
||||
}
|
||||
*/
|
||||
if(i1<0) i1=0;
|
||||
if(i2>n) i2=n;
|
||||
|
||||
*fmin = *fmax = f[i1];
|
||||
f += i1+1;
|
||||
for (i=i1+1; i<i2; i++, f++) {
|
||||
if(*f > *fmax) *fmax = *f;
|
||||
if(*f < *fmin) *fmin = *f;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
template <class TYPE>
|
||||
int
|
||||
find_max(int n, const TYPE *f, int i1, int i2, int *imax, TYPE *fmax)
|
||||
{
|
||||
int i;
|
||||
|
||||
if(i1<0) i1=0;
|
||||
if(i2>n) i2=n;
|
||||
|
||||
*fmax = f[i1];
|
||||
*imax = i1;
|
||||
f += i1+1;
|
||||
for (i=i1+1; i<i2; i++, f++) {
|
||||
if(*f > *fmax) {
|
||||
*imax = i;
|
||||
*fmax = *f;
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
template <typename TYPE>
|
||||
int
|
||||
find_absmax(int n, const TYPE *f, int i1, int i2, TYPE offset=0)
|
||||
{
|
||||
int i;
|
||||
|
||||
if(i1<0) i1=0;
|
||||
if(i2>n) i2=n;
|
||||
|
||||
double fmax = fabs(f[i1]-offset);
|
||||
int imax = i1;
|
||||
|
||||
for (i=i1+1; i<i2; i++) {
|
||||
double ff = fabs(f[i]-offset);
|
||||
if(ff > fmax) {
|
||||
imax = i;
|
||||
fmax = ff;
|
||||
}
|
||||
}
|
||||
|
||||
return imax;
|
||||
}
|
||||
|
||||
template <class TYPE>
|
||||
void
|
||||
find_minmax(int &lmin, int &lmax, int n, const TYPE *f, int i1, int i2, TYPE offset=0)
|
||||
{
|
||||
int i;
|
||||
|
||||
if(i1<0) i1=0;
|
||||
if(i2>n) i2=n;
|
||||
|
||||
TYPE fmax = (f[i1]-offset);
|
||||
TYPE fmin = fmax;
|
||||
|
||||
lmin = i1;
|
||||
lmax = i1;
|
||||
|
||||
for (i=i1+1; i<i2; i++) {
|
||||
TYPE ff = f[i]-offset;
|
||||
if (ff > fmax) {
|
||||
lmax = i;
|
||||
fmax = ff;
|
||||
}
|
||||
else if (ff < fmin) {
|
||||
lmin = i;
|
||||
fmin = ff;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class TYPE>
|
||||
int
|
||||
minmax(std::vector<TYPE> const &f, int i1, int i2, TYPE *fmin, TYPE *fmax)
|
||||
{
|
||||
const TYPE *f0 = &f[0];
|
||||
int n = f.size();
|
||||
return minmax(n, f0, i1, i2, fmin, fmax);
|
||||
}
|
||||
|
||||
template <class TYPE>
|
||||
int
|
||||
find_max(std::vector<TYPE> const &f, int i1, int i2, int *imax, TYPE *fmax)
|
||||
{
|
||||
const TYPE *f0 = &f[0];
|
||||
int n = f.size();
|
||||
return find_max(n, f0, i1, i2, imax, fmax);
|
||||
}
|
||||
78
include/seiscomp/math/misc.ipp
Normal file
78
include/seiscomp/math/misc.ipp
Normal file
@ -0,0 +1,78 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
// namespace is Seiscomp::Filter
|
||||
|
||||
#include<math.h>
|
||||
|
||||
#include<vector>
|
||||
#include<algorithm>
|
||||
using namespace std;
|
||||
|
||||
template<typename TYPE>
|
||||
static vector<TYPE> mytransform(vector<TYPE> const &f, double (*func)(double))
|
||||
{
|
||||
vector<TYPE> y(f.size());
|
||||
transform( f.begin() , f.end() , y.begin() , func) ;
|
||||
return y;
|
||||
}
|
||||
|
||||
template<typename TYPE>
|
||||
vector<TYPE> sin(vector<TYPE> const &f) { return mytransform(f, ::sin); }
|
||||
template<typename TYPE>
|
||||
vector<TYPE> cos(vector<TYPE> const &f) { return mytransform(f, ::cos); }
|
||||
template<typename TYPE>
|
||||
vector<TYPE> tan(vector<TYPE> const &f) { return mytransform(f, ::tan); }
|
||||
template<typename TYPE>
|
||||
vector<TYPE> sqrt(vector<TYPE> const &f) { return mytransform(f, ::sqrt); }
|
||||
template<typename TYPE>
|
||||
vector<TYPE> log(vector<TYPE> const &f) { return mytransform(f, ::log); }
|
||||
template<typename TYPE>
|
||||
vector<TYPE> log10(vector<TYPE> const &f) { return mytransform(f, ::log10); }
|
||||
template<typename TYPE>
|
||||
vector<TYPE> exp(vector<TYPE> const &f) { return mytransform(f, ::exp); }
|
||||
|
||||
template<typename TYPE>
|
||||
static vector<TYPE> arange(TYPE xmax)
|
||||
{
|
||||
int size = int(xmax);
|
||||
size += int(ceil(xmax-size));
|
||||
vector<TYPE> y(size);
|
||||
TYPE *yy = &y[0];
|
||||
for(int i=0; i<size; i++)
|
||||
yy[i] = TYPE(i);
|
||||
return y;
|
||||
}
|
||||
|
||||
#include<iostream>
|
||||
|
||||
int main()
|
||||
{
|
||||
for (int count=0; count<5000; count++)
|
||||
{
|
||||
// vector<double> x(10000, 0.5), y;
|
||||
// y = sin(x);
|
||||
// cerr << x[4] << " " << y[4] << endl; break;
|
||||
vector<double> y;
|
||||
y = arange(10.1);
|
||||
for (int i=0; i<y.size(); i++)
|
||||
cerr << i << " " << y[i]/3 << endl;
|
||||
break;
|
||||
}
|
||||
}
|
||||
82
include/seiscomp/math/polygon.h
Normal file
82
include/seiscomp/math/polygon.h
Normal file
@ -0,0 +1,82 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_GEO_POLYGON_H
|
||||
#define SEISCOMP_MATH_GEO_POLYGON_H
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
#include <utility>
|
||||
#include <seiscomp/math/coord.h>
|
||||
|
||||
|
||||
namespace Seiscomp
|
||||
{
|
||||
|
||||
namespace Math
|
||||
{
|
||||
|
||||
namespace Geo
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
class Polygon : public std::vector< Coord<T> > {
|
||||
public:
|
||||
Polygon();
|
||||
~Polygon();
|
||||
|
||||
/**
|
||||
* Returns whether a location lies inside the polygon
|
||||
* or not.
|
||||
* @param p The location
|
||||
* @return True, if the location lies inside, else false
|
||||
*/
|
||||
bool operator&(const Coord<T>& c);
|
||||
|
||||
void addVertex(double lat, double lon);
|
||||
void addCoord(const Coord<T>& c);
|
||||
|
||||
bool pointInPolygon(const T& lat, const T& lon) const;
|
||||
bool pointInPolygon(const Coord<T>& c) const;
|
||||
|
||||
size_t vertexCount() const;
|
||||
|
||||
const Coord<T>& vertex(int i) const;
|
||||
Coord<T>& vertex(int i);
|
||||
|
||||
void print() const;
|
||||
|
||||
};
|
||||
|
||||
|
||||
typedef Polygon<float> PolygonF;
|
||||
typedef Polygon<double> PolygonD;
|
||||
|
||||
|
||||
} // of ns Geo
|
||||
} // of ns Math
|
||||
} // of ns Seiscomp
|
||||
|
||||
|
||||
template <typename T>
|
||||
std::ostream& operator<<(std::ostream& os, const Seiscomp::Math::Geo::Polygon<T>& poly);
|
||||
|
||||
|
||||
#endif
|
||||
76
include/seiscomp/math/restitution/fft.h
Normal file
76
include/seiscomp/math/restitution/fft.h
Normal file
@ -0,0 +1,76 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_RESTITUTION_FFT_H
|
||||
#define SEISCOMP_MATH_RESTITUTION_FFT_H
|
||||
|
||||
|
||||
#include <complex>
|
||||
#include <vector>
|
||||
#include <seiscomp/math/math.h>
|
||||
#include <seiscomp/math/restitution/types.h>
|
||||
#include <seiscomp/math/restitution/transferfunction.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Restitution {
|
||||
|
||||
|
||||
// Transforms a time series into the spectra, deconvolves it with the
|
||||
// transfer function and transforms the spectra back into the time series.
|
||||
// The spectra is tapered before min_freq and after max_freq. min_freq or
|
||||
// max_freq has to be greater than 0 otherwise the tapering on the
|
||||
// corresponding end is disabled.
|
||||
template <typename T>
|
||||
bool transformFFT(int n, T *inout, double fsamp, const FFT::TransferFunction *tf,
|
||||
double cutoff, double min_freq, double max_freq);
|
||||
|
||||
template <typename T>
|
||||
bool transformFFT(std::vector<T> &inout, double fsamp, const FFT::TransferFunction *tf,
|
||||
double cutoff, double min_freq, double max_freq) {
|
||||
return transformFFT(inout.size(), &inout[0], fsamp, tf, cutoff, min_freq, max_freq);
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
bool transformFFT(int n, T *inout, double fsamp, int n_poles, SeismometerResponse::Pole *poles,
|
||||
int n_zeros, SeismometerResponse::Zero *zeros, double norm,
|
||||
double cutoff, double min_freq, double max_freq) {
|
||||
FFT::PolesAndZeros paz(n_poles, poles, n_zeros, zeros, norm);
|
||||
return transformFFT(n, inout, fsamp, &paz, cutoff, min_freq, max_freq);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
bool transformFFT(std::vector<T> &inout, double fsamp, const Poles &poles,
|
||||
const Zeros &zeros, double norm, double cutoff,
|
||||
double min_freq, double max_freq) {
|
||||
return transformFFT(inout.size(), &inout[0], poles.size(), &poles[0],
|
||||
zeros.size(), &zeros[0], norm, cutoff, min_freq, max_freq);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
117
include/seiscomp/math/restitution/td.h
Normal file
117
include/seiscomp/math/restitution/td.h
Normal file
@ -0,0 +1,117 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
#ifndef SEIS_SIGNAL_TDRESTITUTION_H
|
||||
#define SEIS_SIGNAL_TDRESTITUTION_H
|
||||
|
||||
#include <math.h>
|
||||
#include <seiscomp/math/filter.h>
|
||||
#include <seiscomp/math/filter/butterworth.h>
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Restitution {
|
||||
|
||||
// subroutines to compute parameters for the recursive filter
|
||||
|
||||
// from seismometer eigenperiod T0 and damping parameter h bool
|
||||
bool coefficients_from_T0_h(double fsamp, double gain, double T0, double h, double *c0, double *c1, double *c2);
|
||||
|
||||
// from the two seismometer eigenperiods T1 and T2
|
||||
bool coefficients_from_T1_T2(double fsamp, double gain, double T1, double T2, double *c0, double *c1, double *c2);
|
||||
|
||||
template<typename TYPE>
|
||||
class TimeDomain : public Filtering::InPlaceFilter<TYPE> {
|
||||
public:
|
||||
TimeDomain();
|
||||
~TimeDomain();
|
||||
|
||||
// configuration
|
||||
void setBandpass(int order, double fmin, double fmax);
|
||||
void setSamplingFrequency(double fsamp);
|
||||
void setCoefficients(double c0, double c1, double c2);
|
||||
|
||||
virtual int setParameters(int n, const double *params);
|
||||
|
||||
virtual void reset() {}
|
||||
virtual void apply(int n, TYPE *inout);
|
||||
virtual std::string print() const;
|
||||
|
||||
protected:
|
||||
virtual void init();
|
||||
|
||||
protected:
|
||||
// configuration
|
||||
double c0, c1, c2;
|
||||
double fsamp, dt;
|
||||
double gain;
|
||||
|
||||
private:
|
||||
// filter configuration
|
||||
int order;
|
||||
double fmin, fmax;
|
||||
// temp variables
|
||||
double y0, y1, y2, a1, a2;
|
||||
double cumsum1, cumsum2;
|
||||
|
||||
Filtering::IIR::ButterworthHighLowpass<TYPE> *bandpass;
|
||||
};
|
||||
|
||||
template<typename TYPE>
|
||||
class TimeDomain_from_T0_h: public TimeDomain<TYPE> {
|
||||
public:
|
||||
TimeDomain_from_T0_h(double T0, double h, double gain, double fsamp=0);
|
||||
|
||||
void setBandpass(int order, double fmin, double fmax);
|
||||
virtual std::string print() const;
|
||||
|
||||
Filtering::InPlaceFilter<TYPE>* clone() const;
|
||||
|
||||
protected:
|
||||
virtual void init();
|
||||
|
||||
private:
|
||||
// configuration
|
||||
double T0, h;
|
||||
};
|
||||
|
||||
template<typename TYPE>
|
||||
class TimeDomain_from_T1_T2: public TimeDomain<TYPE> {
|
||||
public:
|
||||
TimeDomain_from_T1_T2(double T1, double T2, double gain, double fsamp=0);
|
||||
|
||||
void setBandpass(int order, double fmin, double fmax);
|
||||
virtual std::string print() const;
|
||||
|
||||
Filtering::InPlaceFilter<TYPE>* clone() const;
|
||||
|
||||
protected:
|
||||
virtual void init();
|
||||
|
||||
private:
|
||||
// configuration
|
||||
double T1, T2;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Seiscomp::Math::Restitution
|
||||
} // namespace Seiscomp::Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
139
include/seiscomp/math/restitution/transferfunction.h
Normal file
139
include/seiscomp/math/restitution/transferfunction.h
Normal file
@ -0,0 +1,139 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_RESTITUTION_TRANSFERFUNCTION_H
|
||||
#define SEISCOMP_MATH_RESTITUTION_TRANSFERFUNCTION_H
|
||||
|
||||
|
||||
#include <seiscomp/core/baseobject.h>
|
||||
#include <seiscomp/math/math.h>
|
||||
#include <seiscomp/math/filter/seismometers.h>
|
||||
#include <seiscomp/math/restitution/types.h>
|
||||
#include <vector>
|
||||
#include <cstdlib>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Restitution {
|
||||
namespace FFT {
|
||||
|
||||
|
||||
DEFINE_SMARTPOINTER(TransferFunction);
|
||||
|
||||
class TransferFunction : public Core::BaseObject {
|
||||
public:
|
||||
//! Returns a transfer function proxy object that acts as
|
||||
//! the product of 'this' and 'a'.
|
||||
//! Both input transfer function objects are referenced
|
||||
//! through a pointer and not managed by the proxy.
|
||||
TransferFunction *operator*(const TransferFunction &a) const;
|
||||
|
||||
//! Returns a transfer function proxy object that acts as
|
||||
//! the quotient of 'this' and 'a'.
|
||||
//! Both input transfer function objects are referenced
|
||||
//! through a pointer and not managed by the proxy.
|
||||
TransferFunction *operator/(const TransferFunction &a) const;
|
||||
|
||||
//! Evaluates the transfer function at nodes x and returns the
|
||||
//! result in out. Out must have enough space for n samples.
|
||||
void evaluate(Complex *out, int n, const double *x) const {
|
||||
evaluate_(out, n, x);
|
||||
}
|
||||
|
||||
void evaluate(std::vector<Complex> &out, const std::vector<double> &x) const {
|
||||
out.resize(x.size());
|
||||
evaluate_(&out[0], (int)x.size(), &x[0]);
|
||||
}
|
||||
|
||||
//! Devides the spectra by the evaluated nodes of the transfer function
|
||||
//! for {startFreq, startFreq + 1*df, ..., startFreq + (n-1)*df}
|
||||
void deconvolve(int n, Complex *spec, double startFreq, double df) const {
|
||||
deconvolve_(n, spec, startFreq, df);
|
||||
}
|
||||
|
||||
//! Convenience wrapper using a vector as output
|
||||
void deconvolve(std::vector<Complex> &spec, double startFreq, double df) const {
|
||||
deconvolve_((int)spec.size(), &spec[0], startFreq, df);
|
||||
}
|
||||
|
||||
//! Multiplies the spectra by the evaluated nodes of the transfer function
|
||||
//! for {startFreq, startFreq + 1*df, ..., startFreq + (n-1)*df}
|
||||
void convolve(int n, Complex *spec, double startFreq, double df) const {
|
||||
convolve_(n, spec, startFreq, df);
|
||||
}
|
||||
|
||||
//! Convenience wrapper using a vector as output
|
||||
void convolve(std::vector<Complex> &spec, double startFreq, double df) const {
|
||||
convolve_((int)spec.size(), &spec[0], startFreq, df);
|
||||
}
|
||||
|
||||
|
||||
protected:
|
||||
//! The implementations
|
||||
virtual void evaluate_(Complex *out, int n, const double *x) const = 0;
|
||||
virtual void deconvolve_(int n, Complex *spec, double startFreq, double df) const;
|
||||
virtual void convolve_(int n, Complex *spec, double startFreq, double df) const;
|
||||
};
|
||||
|
||||
|
||||
class PolesAndZeros : public TransferFunction {
|
||||
public:
|
||||
PolesAndZeros(const SeismometerResponse::PolesAndZeros &polesAndZeros);
|
||||
PolesAndZeros(int n_poles, Pole *poles, int n_zeros, Zero *zeros, double k, int addZeros = 0);
|
||||
|
||||
protected:
|
||||
void evaluate_(Complex *out, int n, const double *x) const;
|
||||
void deconvolve_(int n, Complex *spec, double startFreq, double df) const;
|
||||
void convolve_(int n, Complex *spec, double startFreq, double df) const;
|
||||
|
||||
public:
|
||||
SeismometerResponse::PolesAndZeros paz;
|
||||
};
|
||||
|
||||
|
||||
/**
|
||||
* @brief The ResponseList class evaluates a FAP response list. Phase angles
|
||||
* are expected in degrees.
|
||||
*/
|
||||
class ResponseList : public TransferFunction {
|
||||
public:
|
||||
ResponseList(const SeismometerResponse::FAPs &faps, int addZeros = 0);
|
||||
ResponseList(int n_tuples, const SeismometerResponse::FAP *faps, int addZeros = 0);
|
||||
|
||||
protected:
|
||||
void evaluate_(Complex *out, int n, const double *x) const;
|
||||
void deconvolve_(int n, Complex *spec, double startFreq, double df) const;
|
||||
void convolve_(int n, Complex *spec, double startFreq, double df) const;
|
||||
|
||||
public:
|
||||
SeismometerResponse::FAPs faps;
|
||||
int nZeros;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
43
include/seiscomp/math/restitution/types.h
Normal file
43
include/seiscomp/math/restitution/types.h
Normal file
@ -0,0 +1,43 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_RESTITUTION_TYPES_H
|
||||
#define SEISCOMP_MATH_RESTITUTION_TYPES_H
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
namespace Restitution {
|
||||
|
||||
|
||||
typedef std::complex<double> Pole;
|
||||
typedef std::complex<double> Zero;
|
||||
|
||||
typedef std::vector<Pole> Poles;
|
||||
typedef std::vector<Zero> Zeros;
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
262
include/seiscomp/math/tensor.h
Normal file
262
include/seiscomp/math/tensor.h
Normal file
@ -0,0 +1,262 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_TENSOR_H
|
||||
#define SEISCOMP_MATH_TENSOR_H
|
||||
|
||||
|
||||
#include <seiscomp/math/matrix3.h>
|
||||
#include <stdio.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename T>
|
||||
class Tensor2S;
|
||||
|
||||
|
||||
/**
|
||||
* Non-Symmetric second rank tensor (gradient & rotation tensors)
|
||||
*/
|
||||
template <typename T>
|
||||
class Tensor2N {
|
||||
public:
|
||||
Tensor2N() { assign(0.0); }
|
||||
|
||||
//! Equality operator
|
||||
void assign(Tensor2N<T> &other);
|
||||
//! Initialize
|
||||
void assign(T k);
|
||||
|
||||
//! Make unit tensor
|
||||
void unit();
|
||||
|
||||
//! Add scaled tensor
|
||||
void sum(const Tensor2N<T> &other, T k);
|
||||
|
||||
//! Multiply by a scalar
|
||||
void scale(T k);
|
||||
|
||||
//! Determinant of gradient
|
||||
T det() const;
|
||||
|
||||
//! Make inverse tensor
|
||||
void inverse(const Tensor2N<T>& src);
|
||||
|
||||
void derive(T *dX1,
|
||||
T *dX2,
|
||||
T *dX3,
|
||||
T *disp,
|
||||
int * inodes,
|
||||
int n);
|
||||
|
||||
//! Spin tensor
|
||||
void spin (Tensor2N<T> &U);
|
||||
|
||||
//! FD = [det(F)^-1/3] * F
|
||||
void deviator(Tensor2N<T> &F);
|
||||
|
||||
//! Product of non-symmetric
|
||||
//! tensor with symmetric tensor
|
||||
void product(const Tensor2N<T> &N,
|
||||
const Tensor2S<T> &S);
|
||||
|
||||
//! Product of two non-symmetric
|
||||
//! tensors
|
||||
void product(const Tensor2N<T> &A,
|
||||
const Tensor2N<T> &B);
|
||||
|
||||
//! Performs polar decomposition of the deformation gradient
|
||||
//! tensor into the right stretch tensor and the rotation tensor
|
||||
void polarDecomp(Tensor2N<T> &F);
|
||||
|
||||
//! Hughes-Winget rotation increment
|
||||
void HughesWinget(Tensor2N<T> &U);
|
||||
|
||||
|
||||
public:
|
||||
T _11, _12, _13;
|
||||
T _21, _22, _23;
|
||||
T _31, _32, _33;
|
||||
};
|
||||
|
||||
|
||||
typedef Tensor2N<float> Tensor2Nf;
|
||||
typedef Tensor2N<double> Tensor2Nd;
|
||||
|
||||
|
||||
|
||||
/**
|
||||
* Symmetric second rank tensor (stress & strain tensors)
|
||||
*/
|
||||
template <typename T>
|
||||
class Tensor2S {
|
||||
public:
|
||||
Tensor2S() { assign(0.0); }
|
||||
Tensor2S(T __11, T __12, T __13,
|
||||
T __22, T __23,
|
||||
T __33);
|
||||
|
||||
//! Equality operator
|
||||
void assign(const Tensor2S<T> &A);
|
||||
//! Initialize
|
||||
void assign(T k = 0.0);
|
||||
|
||||
//! Make spheric tensor
|
||||
void spheric(T k = 1.0);
|
||||
|
||||
//! Add scaled tensor
|
||||
void sum(const Tensor2S<T> &A, T k = 1.0);
|
||||
|
||||
//! Multiply by a scalar
|
||||
void scale(T k);
|
||||
|
||||
//! Strain tensor
|
||||
void strain(const Tensor2N<T> &U);
|
||||
|
||||
//! Strain deviator
|
||||
void deviator(const Tensor2N<T> &U);
|
||||
|
||||
//! Green-Lagrange strain tensor
|
||||
void GLStrain(const Tensor2N<T> &F);
|
||||
|
||||
//! Left Cauchy-Green deformation tensor
|
||||
void leftCG(const Tensor2N<T> &F);
|
||||
|
||||
//! Right Cauchy-Green deformation tensor
|
||||
void rightCG(const Tensor2N<T> &F);
|
||||
|
||||
//! Square of symmetric tensor C^2 = C^t*C
|
||||
void square(const Tensor2S<T> &C);
|
||||
|
||||
//! Stress Push-Forward (PF)
|
||||
void pshFrwd(const Tensor2N<T> &F,
|
||||
const Tensor2S<T> &S);
|
||||
|
||||
//! Spatial to unrotated tensor transf.
|
||||
void unrotate(const Matrix3<T> &R);
|
||||
|
||||
//! Spatial to unrotated tensor transf.
|
||||
void unrotate(const Matrix3<T> &R,
|
||||
const Tensor2S<T> &S);
|
||||
|
||||
//! Unrotated to spatial transf.
|
||||
void rotate(const Matrix3<T> &R);
|
||||
|
||||
//! Unrotated to spatial tensor transf.
|
||||
void rotate(const Matrix3<T> &R,
|
||||
const Tensor2S<T> &U);
|
||||
|
||||
//! Trace of tensor
|
||||
T I1() const;
|
||||
|
||||
//! Second invariant of tensor
|
||||
T I2() const;
|
||||
|
||||
//! Mean stress
|
||||
T mean() const;
|
||||
|
||||
//! Euclidean norm of tensor
|
||||
T norm() const;
|
||||
|
||||
//! Euclidean norm of deviator
|
||||
T devnorm() const;
|
||||
|
||||
//! Make deviator & return trace
|
||||
T dtrace();
|
||||
|
||||
//! Make deviator & return mean value
|
||||
T dmean();
|
||||
|
||||
//! Augment deviator with spheric tensor
|
||||
void augment(T p);
|
||||
|
||||
//! Update deviatoric stress
|
||||
void update(Tensor2S<T> &d,
|
||||
T G);
|
||||
|
||||
//! Prints tensor to file
|
||||
void print(FILE * fl);
|
||||
|
||||
//! apply Jacoby rotation
|
||||
void jacoby(T &pp,
|
||||
T &qq,
|
||||
T &pq,
|
||||
T &rp,
|
||||
T &rq);
|
||||
|
||||
//! compute eigenvalues
|
||||
bool eigenval(T tol = 1e-15,
|
||||
int itmax = 30);
|
||||
|
||||
public:
|
||||
T _11, _12, _13;
|
||||
T _22, _23;
|
||||
T _33;
|
||||
};
|
||||
|
||||
|
||||
typedef Tensor2S<float> Tensor2Sf;
|
||||
typedef Tensor2S<double> Tensor2Sd;
|
||||
|
||||
|
||||
/**
|
||||
* Spectral representation of symmetric tensor
|
||||
*/
|
||||
template <typename T>
|
||||
class Spectral2S {
|
||||
public:
|
||||
Spectral2S() { assign(0.0); }
|
||||
|
||||
//! initialize
|
||||
void assign(T k);
|
||||
|
||||
//! spectral decomposition
|
||||
bool spect(const Tensor2S<T> &A,
|
||||
T atol = 1e-12,
|
||||
int itmax = 50);
|
||||
//! sort principal values in decending order & permute principal directions
|
||||
void sort();
|
||||
|
||||
//! sort absolute principal values in decending order & permute principal directions
|
||||
void absSort();
|
||||
|
||||
//! compose symmetric tensor from spectral representation
|
||||
void compose(Tensor2S<T> &A);
|
||||
|
||||
//! compute norm of spectrally decomposed tensor
|
||||
T norm();
|
||||
|
||||
public:
|
||||
T a1, a2, a3; // principal values of tensor
|
||||
Vector3<T> n1, n2, n3; // principal directions of tensor
|
||||
};
|
||||
|
||||
|
||||
typedef Spectral2S<float> Spectral2Sf;
|
||||
typedef Spectral2S<double> Spectral2Sd;
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
75
include/seiscomp/math/vector3.h
Normal file
75
include/seiscomp/math/vector3.h
Normal file
@ -0,0 +1,75 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_VECTOR3_H
|
||||
#define SEISCOMP_MATH_VECTOR3_H
|
||||
|
||||
|
||||
#include <seiscomp/math/math.h>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename T>
|
||||
struct Vector3 {
|
||||
Vector3() {}
|
||||
Vector3(const Vector3<T> &other) : x(other.x), y(other.y), z(other.z) {}
|
||||
Vector3(T x_, T y_, T z_) : x(x_), y(y_), z(z_) {}
|
||||
|
||||
T length() const;
|
||||
T dot(const Vector3<T> &v) const;
|
||||
|
||||
Vector3<T> &cross(const Vector3<T> &a, const Vector3<T> &b);
|
||||
Vector3<T> &normalize();
|
||||
|
||||
Vector3<T> &operator=(const Vector3<T> &other);
|
||||
Vector3<T> operator*(T scale) const;
|
||||
Vector3<T> &operator*=(T scale);
|
||||
T operator*(const Vector3<T> &other) const;
|
||||
|
||||
operator T *() { return (T*)this; }
|
||||
operator const T *() const { return (const T*)this; }
|
||||
|
||||
Vector3<T> &operator+=(const Vector3<T> &other);
|
||||
Vector3<T> &operator-=(const Vector3<T> &other);
|
||||
|
||||
Vector3<T> operator+(const Vector3<T> &other) const;
|
||||
Vector3<T> operator-(const Vector3<T> &other) const;
|
||||
|
||||
Vector3<T> &fromAngles(T radAzimuth, T radDip);
|
||||
Vector3<T> &toAngles(T &radAzimuth, T &radDip);
|
||||
|
||||
T x,y,z;
|
||||
};
|
||||
|
||||
|
||||
typedef Vector3<float> Vector3f;
|
||||
typedef Vector3<double> Vector3d;
|
||||
|
||||
|
||||
#include <seiscomp/math/vector3.ipp>
|
||||
|
||||
|
||||
} // namespace Math
|
||||
} // namespace Seiscomp
|
||||
|
||||
#endif
|
||||
39
include/seiscomp/math/vector3.ipp
Normal file
39
include/seiscomp/math/vector3.ipp
Normal file
@ -0,0 +1,39 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
template <typename T>
|
||||
inline T Vector3<T>::length() const {
|
||||
return sqrt(x*x + y*y + z*z);
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
inline T Vector3<T>::dot(const Vector3<T> &v) const {
|
||||
return x*v.x + y*v.y + z*v.z;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
inline Vector3<T> &Vector3<T>::operator=(const Vector3<T> &other) {
|
||||
x = other.x;
|
||||
y = other.y;
|
||||
z = other.z;
|
||||
|
||||
return *this;
|
||||
}
|
||||
135
include/seiscomp/math/windowfunc.h
Normal file
135
include/seiscomp/math/windowfunc.h
Normal file
@ -0,0 +1,135 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_WINDOWFUNC_H
|
||||
#define SEISCOMP_MATH_WINDOWFUNC_H
|
||||
|
||||
|
||||
#include <seiscomp/core/baseobject.h>
|
||||
#include <seiscomp/core/typedarray.h>
|
||||
#include <seiscomp/core.h>
|
||||
|
||||
#include <vector>
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename TYPE>
|
||||
class SC_SYSTEM_CORE_API WindowFunc : public Core::BaseObject {
|
||||
public:
|
||||
virtual ~WindowFunc();
|
||||
|
||||
|
||||
public:
|
||||
/**
|
||||
* @brief Applies the window function to the given data
|
||||
* @param n Number of samples
|
||||
* @param inout The data vector where each sample is multiplied with
|
||||
* the respective sample of the window function.
|
||||
* @param width The width of the window function. The default is 0.5
|
||||
* which means 50% at either side. 0.1 would mean that the
|
||||
* left half of the window function is applied on 10% of
|
||||
* the left portion of the data vector and the right half
|
||||
* of the window function is applied on the right 10% of
|
||||
* the data vector. The value is clipped into range [0,0.5].
|
||||
*/
|
||||
void apply(int n, TYPE *inout, double width = 0.5) const;
|
||||
|
||||
/**
|
||||
* @brief Applies the window function to the given data
|
||||
* @param inout The data vector where each sample is multiplied with
|
||||
* the respective sample of the window function.
|
||||
* @param width The width of the window function. The default is 0.5
|
||||
* which means 50% at either side. 0.1 would mean that the
|
||||
* left half of the window function is applied on 10% of
|
||||
* the left portion of the data vector and the right half
|
||||
* of the window function is applied on the right 10% of
|
||||
* the data vector. The value is clipped into range [0,0.5].
|
||||
*/
|
||||
void apply(std::vector<TYPE> &inout, double width = 0.5) const;
|
||||
|
||||
/**
|
||||
* @brief Applies the window function to the given data
|
||||
* @param inout The data array where each sample is multiplied with
|
||||
* the respective sample of the window function.
|
||||
* @param width The width of the window function. The default is 0.5
|
||||
* which means 50% at either side. 0.1 would mean that the
|
||||
* left half of the window function is applied on 10% of
|
||||
* the left portion of the data vector and the right half
|
||||
* of the window function is applied on the right 10% of
|
||||
* the data vector. The value is clipped into range [0,0.5].
|
||||
*/
|
||||
void apply(TypedArray<TYPE> &inout, double width = 0.5) const;
|
||||
|
||||
/**
|
||||
* @brief Applies the window function to the given data
|
||||
* @param inout The data array where each sample is multiplied with
|
||||
* the respective sample of the window function.
|
||||
* @param width The width of the window function. The default is 0.5
|
||||
* which means 50% at either side. 0.1 would mean that the
|
||||
* left half of the window function is applied on 10% of
|
||||
* the left portion of the data vector and the right half
|
||||
* of the window function is applied on the right 10% of
|
||||
* the data vector. The value is clipped into range [0,0.5].
|
||||
*/
|
||||
void apply(TypedArray<TYPE> *inout, double width = 0.5) const;
|
||||
|
||||
//! Apply methods for non-symmetric window lengths.
|
||||
void apply(int n, TYPE *inout, double left, double right) const;
|
||||
void apply(std::vector<TYPE> &inout, double left, double right) const;
|
||||
void apply(TypedArray<TYPE> &inout, double left, double right) const;
|
||||
void apply(TypedArray<TYPE> *inout, double left, double right) const;
|
||||
|
||||
|
||||
protected:
|
||||
/**
|
||||
* @brief Applies the window function to the given data. This method has
|
||||
* to be implemented by derived classes. It is called by all
|
||||
* apply variants.
|
||||
* @param n Number of samples
|
||||
* @param inout The data vector where each sample is multiplied with
|
||||
* the respective sample of the window function.
|
||||
* @param left The width of the window function at the left side.
|
||||
* The default is 0.5 which means 50% at the left side.
|
||||
* 0.1 would mean that the left half of the window function
|
||||
* is applied on 10% of the left portion of the data vector.
|
||||
* The value is clipped into range [0,0.5].
|
||||
* @param right The width of the window function at the right side
|
||||
* respectively.
|
||||
*/
|
||||
virtual void process(int n, TYPE *inout,
|
||||
double left = 0.5,
|
||||
double right = 0.5) const = 0;
|
||||
|
||||
|
||||
private:
|
||||
void checkAndProcess(int n, TYPE *inout,
|
||||
double left = 0.5,
|
||||
double right = 0.5) const;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
43
include/seiscomp/math/windows/bartlett.h
Normal file
43
include/seiscomp/math/windows/bartlett.h
Normal file
@ -0,0 +1,43 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_WINDOWS_BARTLETT_H
|
||||
#define SEISCOMP_MATH_WINDOWS_BARTLETT_H
|
||||
|
||||
|
||||
#include "../windowfunc.h"
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename TYPE>
|
||||
class SC_SYSTEM_CORE_API BartlettWindow : public WindowFunc<TYPE> {
|
||||
protected:
|
||||
virtual void process(int n, TYPE *inout, double left = 0.5, double right = 0.5) const override;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
51
include/seiscomp/math/windows/blackman.h
Normal file
51
include/seiscomp/math/windows/blackman.h
Normal file
@ -0,0 +1,51 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_WINDOWS_BLACKMAN_H
|
||||
#define SEISCOMP_MATH_WINDOWS_BLACKMAN_H
|
||||
|
||||
|
||||
#include "../windowfunc.h"
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename TYPE>
|
||||
class SC_SYSTEM_CORE_API BlackmanWindow : public WindowFunc<TYPE> {
|
||||
public:
|
||||
BlackmanWindow(TYPE alpha = 0.16);
|
||||
|
||||
|
||||
protected:
|
||||
virtual void process(int n, TYPE *inout, double left = 0.5, double right = 0.5) const override;
|
||||
|
||||
|
||||
private:
|
||||
TYPE _alpha;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
43
include/seiscomp/math/windows/cosine.h
Normal file
43
include/seiscomp/math/windows/cosine.h
Normal file
@ -0,0 +1,43 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_WINDOWS_COSINE_H
|
||||
#define SEISCOMP_MATH_WINDOWS_COSINE_H
|
||||
|
||||
|
||||
#include "../windowfunc.h"
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename TYPE>
|
||||
class SC_SYSTEM_CORE_API CosineWindow : public WindowFunc<TYPE> {
|
||||
protected:
|
||||
virtual void process(int n, TYPE *inout, double left = 0.5, double right = 0.5) const override;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
43
include/seiscomp/math/windows/hamming.h
Normal file
43
include/seiscomp/math/windows/hamming.h
Normal file
@ -0,0 +1,43 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_WINDOWS_HAMMING_H
|
||||
#define SEISCOMP_MATH_WINDOWS_HAMMING_H
|
||||
|
||||
|
||||
#include "../windowfunc.h"
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename TYPE>
|
||||
class SC_SYSTEM_CORE_API HammingWindow : public WindowFunc<TYPE> {
|
||||
protected:
|
||||
virtual void process(int n, TYPE *inout, double left = 0.5, double right = 0.5) const override;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
43
include/seiscomp/math/windows/hann.h
Normal file
43
include/seiscomp/math/windows/hann.h
Normal file
@ -0,0 +1,43 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_WINDOWS_HANN_H
|
||||
#define SEISCOMP_MATH_WINDOWS_HANN_H
|
||||
|
||||
|
||||
#include "../windowfunc.h"
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename TYPE>
|
||||
class SC_SYSTEM_CORE_API HannWindow : public WindowFunc<TYPE> {
|
||||
protected:
|
||||
virtual void process(int n, TYPE *inout, double left = 0.5, double right = 0.5) const override;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
43
include/seiscomp/math/windows/triangle.h
Normal file
43
include/seiscomp/math/windows/triangle.h
Normal file
@ -0,0 +1,43 @@
|
||||
/***************************************************************************
|
||||
* Copyright (C) gempa GmbH *
|
||||
* All rights reserved. *
|
||||
* Contact: gempa GmbH (seiscomp-dev@gempa.de) *
|
||||
* *
|
||||
* GNU Affero General Public License Usage *
|
||||
* This file may be used under the terms of the GNU Affero *
|
||||
* Public License version 3.0 as published by the Free Software Foundation *
|
||||
* and appearing in the file LICENSE included in the packaging of this *
|
||||
* file. Please review the following information to ensure the GNU Affero *
|
||||
* Public License version 3.0 requirements will be met: *
|
||||
* https://www.gnu.org/licenses/agpl-3.0.html. *
|
||||
* *
|
||||
* Other Usage *
|
||||
* Alternatively, this file may be used in accordance with the terms and *
|
||||
* conditions contained in a signed written agreement between you and *
|
||||
* gempa GmbH. *
|
||||
***************************************************************************/
|
||||
|
||||
|
||||
#ifndef SEISCOMP_MATH_WINDOWS_TRIANGLE_H
|
||||
#define SEISCOMP_MATH_WINDOWS_TRIANGLE_H
|
||||
|
||||
|
||||
#include "../windowfunc.h"
|
||||
|
||||
|
||||
namespace Seiscomp {
|
||||
namespace Math {
|
||||
|
||||
|
||||
template <typename TYPE>
|
||||
class SC_SYSTEM_CORE_API TriangleWindow : public WindowFunc<TYPE> {
|
||||
protected:
|
||||
virtual void process(int n, TYPE *inout, double left = 0.5, double right = 0.5) const override;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user