Commit 86dd3cb5 authored by Florian Wittkamp's avatar Florian Wittkamp
Browse files

Merge branch 'import_Seitosh' into 'Seitosh/import'

Import seitosh

@wittkamp @lgass @r 
[libstfinv](https://git.scc.kit.edu/Seitosh/Seitosh/tree/master/src/libs/libstfinv) in [Seitosh](https://git.scc.kit.edu/Seitosh/Seitosh) provides a new feature:

A new taper option for procedures operating in the Fourier domain is implemenetd in libstfinv. This can be used to remove unwanted acausal or coda components of the impulse response of the source time function correction filter. See usage information provided by libstfinv:

    A time domain taper can be applied to the impulse response of the correction
    filter by using option irtap. Four time values are given in units of seconds:
    t1, t2, t3, and t4. They must be in increasing order and (t4-t1) must be
    smaller than the total duration of the time series used to represent signals
    internally. Times value are allowed to be negative. Time series are understood
    to be periodic (due to discrete Fourier transformation). Prior to application
    of the correction filter to the time series passed to the algorithm, the
    correction filter is transformed to the time domain, tapered, and then
    transformed to the Fourier domain again. The values of the taper are:
      0                              if       t <  t1
      0.5-0.5*cos(pi*(t-t1)/(t2-t1)) if t1 <= t <= t2
      1                              if t2 <  t <  t3
      0.5+0.5*cos(pi*(t-t3)/(t3-t4)) if t3 <= t <= t4
      0                              if       t >  t4

Test-cases are set up for [soutifu](https://git.scc.kit.edu/Seitosh/Seitosh/blob/master/src/ts/wf/soutifu.cc) and are described in the accompanying [README.soutifu](https://git.scc.kit.edu/Seitosh/Seitosh/blob/master/src/ts/wf/testcases/README.soutifu).

A vendor update has taken place in [thof/IFOS2D](https://git.scc.kit.edu/thof/IFOS2D) see [issue 5](https://git.scc.kit.edu/thof/IFOS2D/issues/5) there. The code in contrib compiles and installs well. The integration into IFOS2D still is to be tested.

See merge request !3
parents abdeb652 f8fb4695
...@@ -3,7 +3,6 @@ ...@@ -3,7 +3,6 @@
* *
* ---------------------------------------------------------------------------- * ----------------------------------------------------------------------------
* *
* $Id$
* \author Thomas Forbriger * \author Thomas Forbriger
* \date 05/01/2003 * \date 05/01/2003
* *
...@@ -38,8 +37,6 @@ ...@@ -38,8 +37,6 @@
#define TF_POLESNZEROES_H_VERSION \ #define TF_POLESNZEROES_H_VERSION \
"TF_POLESNZEROES_H V1.0 " "TF_POLESNZEROES_H V1.0 "
#define TF_POLESNZEROES_H_CVSID \
"$Id$"
#include<cmath> #include<cmath>
#include<complex> #include<complex>
......
c this is <polesnzeros.inc> c this is <polesnzeros.inc>
cS cS
c ---------------------------------------------------------------------------- c ----------------------------------------------------------------------------
c ($Id$)
c c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt) c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c c
......
this is <README> this is <README>
============================================================================ ============================================================================
SEIFE --- seismic waveform filters SEIFE --- seismic waveform filters
----------------------------------
$Id$
============================================================================ ============================================================================
For compilation instructions see README.1st in the root directory of the tar-ball or For compilation instructions see README.1st in the root directory.
http://gpitrsvn.gpi.uni-karlsruhe.de:8000/TFSoftware/wiki/docs/installation
libseife is a collection of Fortran 77 functions and subroutines for time libseife is a collection of Fortran 77 functions and subroutines for time
series analysis and digital filters. It is a full-grown signal processing series analysis and digital filters. It is a full-grown signal processing
...@@ -37,5 +34,15 @@ Dependencies ...@@ -37,5 +34,15 @@ Dependencies
external dependencies: - external dependencies: -
internal dependencies: - internal dependencies: -
Arrays in C version
-------------------
The proper definition of arrays and parameters to be passed to the C-API of
libseife becomes obvious in function seife_first, which is defined in cseife.c
A call
seife_first(x, n);
will address elements x[0] to x[n-1] in C-array x.
----- END OF README ----- ----- END OF README -----
...@@ -70,7 +70,7 @@ documentation. ...@@ -70,7 +70,7 @@ documentation.
User documentation User documentation
------------------ ------------------
The theory behind the Fourier domain least squares engine is outlined by The theory behind the Fourier domain least squares procedure is outlined by
Lisa Groos (2013, Appendix F, page 146). She also describes a way to find an Lisa Groos (2013, Appendix F, page 146). She also describes a way to find an
approrpiate water level by application of the L-curve criterion (Groos, 2013, approrpiate water level by application of the L-curve criterion (Groos, 2013,
Appendix G, page 148). Appendix G, page 148).
...@@ -87,6 +87,9 @@ tests/onlinehelp Provides access to these texts. Just issue ...@@ -87,6 +87,9 @@ tests/onlinehelp Provides access to these texts. Just issue
to get access. to get access.
A toy example and a step-by-step introduction are provided in subdirectory
src/ts/wf/testcases in README.soutifu
A short descrpition of the library and the accompanying program soutifu is A short descrpition of the library and the accompanying program soutifu is
provided on the OpenTOAST web-page: provided on the OpenTOAST web-page:
http://www.opentoast.de/Data_analysis_code_soutifu_and_libstfinv.php http://www.opentoast.de/Data_analysis_code_soutifu_and_libstfinv.php
......
...@@ -40,7 +40,7 @@ ...@@ -40,7 +40,7 @@
"STFINV_DEBUG_H V1.0" "STFINV_DEBUG_H V1.0"
/*! \brief produce debug output /*! \brief produce debug output
* \ingroup misc_h * \ingroup group_debug
* *
* \param C output will be generated if C == true * \param C output will be generated if C == true
* \param N name of function * \param N name of function
...@@ -54,8 +54,15 @@ ...@@ -54,8 +54,15 @@
std::cerr.flush(); \ std::cerr.flush(); \
} }
/*! \brief report value in a sequence of output operators
* \ingroup group_debug
*
* \param P parameter to dump
*/
#define STFINV_value( P ) #P << "=" << P
/*! \brief report value of expression /*! \brief report value of expression
* \ingroup misc_h * \ingroup group_debug
* *
* \param P parameter to dump * \param P parameter to dump
*/ */
......
...@@ -104,7 +104,7 @@ processing are described in the \ref page_help. ...@@ -104,7 +104,7 @@ processing are described in the \ref page_help.
/*! \brief Engines implemented in libstfinv /*! \brief Engines implemented in libstfinv
\defgroup engines Engines \defgroup group_engines Engines
\todo \todo
A detailed description for implementers is still missing A detailed description for implementers is still missing
...@@ -116,7 +116,7 @@ processing are described in the \ref page_help. ...@@ -116,7 +116,7 @@ processing are described in the \ref page_help.
/*! \brief Tools and utilities used by the libstfinv engines /*! \brief Tools and utilities used by the libstfinv engines
\defgroup tools Internal tools and utilities \defgroup group_tools Internal tools and utilities
\todo \todo
A detailed description for implementers is still missing A detailed description for implementers is still missing
...@@ -136,4 +136,18 @@ processing are described in the \ref page_help. ...@@ -136,4 +136,18 @@ processing are described in the \ref page_help.
\date 04.10.2015 \date 04.10.2015
*/ */
/*======================================================================*/
/*! \brief Debug tools
\defgroup group_debug Debugging module
\todo
Selection of debug statements in most parts of the code is done by bits in
the value of the debug variable. E.g∵ The value 16 (bit 4) selects debugging
of the taper function. This is not yet properly documented.
\date 04.10.2015
*/
// ----- END OF doxygen.txt ----- // ----- END OF doxygen.txt -----
...@@ -94,9 +94,9 @@ ...@@ -94,9 +94,9 @@
Programs using this library will require the following libraries in Programs using this library will require the following libraries in
addition: addition:
- libfourierxx (available from TFSoftware) - libfourierxx (available from Seitosh)
- libfftw3 - libfftw3
- libaff (available from TFSoftware) - libaff (available from Seitosh)
C programs will further be required to link against C programs will further be required to link against
......
...@@ -45,12 +45,12 @@ ...@@ -45,12 +45,12 @@
namespace stfinv { namespace stfinv {
/*! \brief Namespace for internal tools /*! \brief Namespace for internal tools
* \ingroup tools * \ingroup group_tools
*/ */
namespace tools { namespace tools {
/*! strip substring /*! strip substring
* \ingroup tools * \ingroup group_tools
* *
* Strips off first substring up to given delimiter. * Strips off first substring up to given delimiter.
* The string is passed as a reference and will be modified (i.e. the * The string is passed as a reference and will be modified (i.e. the
...@@ -65,14 +65,14 @@ namespace stfinv { ...@@ -65,14 +65,14 @@ namespace stfinv {
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*! \brief A map to store parameters. /*! \brief A map to store parameters.
* \ingroup tools * \ingroup group_tools
*/ */
typedef std::map<std::string,std::string> Tparamap; typedef std::map<std::string,std::string> Tparamap;
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*! \brief Create a parameter map from a parameter string /*! \brief Create a parameter map from a parameter string
* \ingroup tools * \ingroup group_tools
* *
* \param p parameter string * \param p parameter string
* \param delimiter delimiter which separates two parameters * \param delimiter delimiter which separates two parameters
...@@ -86,7 +86,7 @@ namespace stfinv { ...@@ -86,7 +86,7 @@ namespace stfinv {
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*! replace comma by whitespace /*! replace comma by whitespace
* \ingroup tools * \ingroup group_tools
* *
* \param s input string * \param s input string
* \return input string with all commas replaced by whitespace * \return input string with all commas replaced by whitespace
...@@ -96,7 +96,7 @@ namespace stfinv { ...@@ -96,7 +96,7 @@ namespace stfinv {
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*! \brief remove leading and trailing whitespace /*! \brief remove leading and trailing whitespace
* \ingroup tools * \ingroup group_tools
* *
* \param s any string * \param s any string
* \return value a input string with any leading and trailing whitespace * \return value a input string with any leading and trailing whitespace
......
...@@ -47,7 +47,7 @@ ...@@ -47,7 +47,7 @@
namespace stfinv { namespace stfinv {
/*! \brief Fourier domain least squares engine /*! \brief Fourier domain least squares engine
* \ingroup engines * \ingroup group_engines
* *
* \par Concept behind this engine * \par Concept behind this engine
* If * If
......
...@@ -30,7 +30,7 @@ ...@@ -30,7 +30,7 @@
namespace stfinv { namespace stfinv {
/*! \brief Engine to find a finite, causal source time-history in time domain /*! \brief Engine to find a finite, causal source time-history in time domain
* \ingroup engines * \ingroup group_engines
* *
* \par Concept behin this engine * \par Concept behin this engine
* *
......
...@@ -44,7 +44,7 @@ ...@@ -44,7 +44,7 @@
namespace stfinv { namespace stfinv {
/*! \brief Engine to provide a fixed wavelet /*! \brief Engine to provide a fixed wavelet
* \ingroup engines * \ingroup group_engines
*/ */
class STFEngineFixedWavelet: public stfinv::STFFourierDomainEngine { class STFEngineFixedWavelet: public stfinv::STFFourierDomainEngine {
public: public:
......
...@@ -32,13 +32,20 @@ ...@@ -32,13 +32,20 @@
* - 30/09/2011 V1.1 implemented handling of additional time series pairs * - 30/09/2011 V1.1 implemented handling of additional time series pairs
* - 04/10/2011 V1.2 correction in debug message * - 04/10/2011 V1.2 correction in debug message
* - 14/10/2015 V1.3 new end-user usage functions * - 14/10/2015 V1.3 new end-user usage functions
* - 28/06/2016 V1.4 introduce taper parameter to taper correction filter
* response in the time domain
* - 22/07/2016 V1.5 thof:
* - provide separate FFT processor addressing just the
* source time function correction filter
* - implement taper function
* *
* ============================================================================ * ============================================================================
*/ */
#define STFINV_STFINVFOURIER_CC_VERSION \ #define STFINV_STFINVFOURIER_CC_VERSION \
"STFINV_STFINVFOURIER_CC V1.3" "STFINV_STFINVFOURIER_CC V1.5"
#include <sstream> #include <sstream>
#include <cmath>
#include <stfinv/stfinvfourier.h> #include <stfinv/stfinvfourier.h>
#include <stfinv/stfinvfourier_summary_usage.h> #include <stfinv/stfinvfourier_summary_usage.h>
#include <stfinv/stfinvfourier_description_usage.h> #include <stfinv/stfinvfourier_description_usage.h>
...@@ -165,13 +172,21 @@ namespace stfinv { ...@@ -165,13 +172,21 @@ namespace stfinv {
* next integer multiple of the divisor larger than fpad*nsamples. * next integer multiple of the divisor larger than fpad*nsamples.
* *
* \par Workspace * \par Workspace
* Two FFT engines will be created: * Two FFT engines will be created.
* They are only used one-way, since input data is not altered during
* processing input time series only have to be transformed once to Fourier
* domain.
* Results of processing have to be transformed to time domain.
* -# One engine (\c STFFourierDomainEngine::Mfftengineinput) being shared * -# One engine (\c STFFourierDomainEngine::Mfftengineinput) being shared
* by recorded data and synthetic data, because both have to be transformed * by recorded data and synthetic data, because both have to be transformed
* to Fourier domain at once. * to Fourier domain at once.
* Transformation to Fourier domain takes place in
* STFFourierDomainEngine::fftinput.
* -# One engine (\c STFFourierDomainEngine::Mfftengineoutput) being shared * -# One engine (\c STFFourierDomainEngine::Mfftengineoutput) being shared
* by the stf and the convolved synthetics, because both have to be * by the stf and the convolved synthetics, because both have to be
* transformed to the time domain at once. * transformed to the time domain at once.
* Transformation to time domain takes place in
* STFFourierDomainEngine::fftoutput.
*/ */
void STFFourierDomainEngine::initialize() void STFFourierDomainEngine::initialize()
{ {
...@@ -226,6 +241,8 @@ namespace stfinv { ...@@ -226,6 +241,8 @@ namespace stfinv {
Mfftengineoutput=fourier::fft::DRFFTWAFFArrayEngine( Mfftengineoutput=fourier::fft::DRFFTWAFFArrayEngine(
1+this->nreceivers()+this->npairs(), 1+this->nreceivers()+this->npairs(),
nsamples); nsamples);
Mfftenginestf=fourier::fft::DRFFTWAFFArrayEngine(this->stfseries(),
this->stfspec());
// set time shift for STF if requested // set time shift for STF if requested
// ----------------------------------- // -----------------------------------
...@@ -235,27 +252,69 @@ namespace stfinv { ...@@ -235,27 +252,69 @@ namespace stfinv {
} }
Mapplyshift=this->parameterisset("tshift"); Mapplyshift=this->parameterisset("tshift");
// taper definition
// ----------------
{
std::istringstream is(stfinv::tools::commatospace(
this->parameter("irtap","0.,1.,2.,3.")));
is >> Mtt1 >> Mtt2 >> Mtt3 >> Mtt4;
}
Mapplystftaper=this->parameterisset("irtap");
STFINV_assert((Mtt1<=Mtt2) && (Mtt2<=Mtt3) && (Mtt3<=Mtt4),
"ERROR: taper definition times not in increasing order");
STFINV_debug(Mdebug&1, "STFFourierDomainEngine::initialize()", STFINV_debug(Mdebug&1, "STFFourierDomainEngine::initialize()",
"this->nsamples()=" << this->nsamples() << " " << STFINV_value(this->nsamples()) << "\n " <<
"padfactor=" << padfactor << " " << STFINV_value(padfactor) << "\n " <<
"divisor=" << divisor << " " << STFINV_value(divisor) << "\n " <<
"nsamples=" << nsamples << " " << STFINV_value(nsamples) << "\n " <<
"Mapplyshift=" << Mapplyshift << " " << STFINV_value(Mapplyshift) << "\n " <<
"Mtshift=" << Mtshift << " "); STFINV_value(Mtshift) << "\n " <<
STFINV_value(Mapplystftaper) << "\n " <<
STFINV_value(Mtt1) << "\n " <<
STFINV_value(Mtt2) << "\n " <<
STFINV_value(Mtt3) << "\n " <<
STFINV_value(Mtt4));
} // void STFFourierDomainEngine::initialize() } // void STFFourierDomainEngine::initialize()
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*!
* Provide data passed (through the API) by the caller of the library to an
* engine operating in the Fourier domain.
* All input is available as time series data in first place.
* -# Copy time series data to STFFourierDomainEngine::Mfftengineinput
* -# Transform time series data to Fourier domain
* -# Clear output coefficients
*
* This function should be called by the very first statement of the
* exec-function of the derived class (e.g. STFEngineFDLeastSquares::exec).
*/
void STFFourierDomainEngine::fftinput() void STFFourierDomainEngine::fftinput()
{ {
this->getinput(); this->getinput();
Mfftengineinput.r2c(); Mfftengineinput.r2c();
Mfftengineoutput.series()=0.;
} // void STFFourierDomainEngine::fftinput() } // void STFFourierDomainEngine::fftinput()
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*!
* Provide results of Fourier domain operation to the caller of the library.
* -# Apply time domain taper to correction filter, if requested
* -# Apply convolution with correction filter response to all synthetic
* input data
* -# Apply time shift to impulse response of correction filter, if
* requested
* -# Transform Fourier series to time domain
* -# Copy time series data to buffer array accessible through API
*
* This function should be called by the very last statement of the
* exec-function of the derived class (e.g. STFEngineFDLeastSquares::exec).
*/
void STFFourierDomainEngine::fftoutput() void STFFourierDomainEngine::fftoutput()
{ {
if (this->Mapplystftaper) { this->taperstf(); }
this->convolve(); this->convolve();
if (Mapplyshift) { this->stfshift(); } if (Mapplyshift) { this->stfshift(); }
Mfftengineoutput.c2r(); Mfftengineoutput.c2r();
...@@ -296,6 +355,13 @@ namespace stfinv { ...@@ -296,6 +355,13 @@ namespace stfinv {
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
STFFourierDomainEngine::TAseries STFFourierDomainEngine::stfseries() const
{
return(Mfftengineoutput.series(this->nreceivers()));
} // STFFourierDomainEngine::Tspectrum STFFourierDomainEngine::stfseries() const
/*----------------------------------------------------------------------*/
STFFourierDomainEngine::TAspectrum STFFourierDomainEngine::TAspectrum
STFFourierDomainEngine::recordingcoeff(const unsigned int& i) const STFFourierDomainEngine::recordingcoeff(const unsigned int& i) const
{ {
...@@ -353,6 +419,8 @@ namespace stfinv { ...@@ -353,6 +419,8 @@ namespace stfinv {
void STFFourierDomainEngine::getinput() void STFFourierDomainEngine::getinput()
{ {
// clear workspace through reference // clear workspace through reference
// (remark: return value sarray is a reference to an array addressing all
// samples of all time series contained in Mfftengineinput)
TAseries sarray=Mfftengineinput.series(); TAseries sarray=Mfftengineinput.series();
sarray=0.; sarray=0.;
...@@ -379,6 +447,11 @@ namespace stfinv { ...@@ -379,6 +447,11 @@ namespace stfinv {
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*!
* Read out time series data from Mfftengineoutput.
* The time series containers in Mfftengineoutput are expected to be larger
* than those used for the API of the library (because of padding).
*/
void STFFourierDomainEngine::putoutput() void STFFourierDomainEngine::putoutput()
{ {
// scaling factor for source correction filter // scaling factor for source correction filter
...@@ -415,6 +488,18 @@ namespace stfinv { ...@@ -415,6 +488,18 @@ namespace stfinv {
/*----------------------------------------------------------------------*/ /*----------------------------------------------------------------------*/
/*! \brief Convolve all synthetics with the source time function correction
* filter.
*
* The filter response is available in terms of the Fourier coefficients of
* the correction filter as are returned by
* STFFourierDomainEngine::stfspec().
* Convolution takes place in the Fourier domain simply by a multiplication
* of Fourier coefficients.
* Since all Fourier coefficients are expected to be scaled appropriately to
* represent samples for the Fourier integral transform, no scaling has to
* be applied here.
*/
void STFFourierDomainEngine::convolve() void STFFourierDomainEngine::convolve()
{ {
TAspectrum synthetic=this->syntheticspec(); TAspectrum synthetic=this->syntheticspec();
...@@ -483,6 +568,96 @@ namespace stfinv { ...@@ -483,6 +568,96 @@ namespace stfinv {
} // for (unsigned int j=0; j<this->nfreq(); ++j) } // for (unsigned int j=0; j<this->nfreq(); ++j)
} // void STFFourierDomainEngine::stfshift() } // void STFFourierDomainEngine::stfshift()
/*----------------------------------------------------------------------*/
/*!
* Apply a time domain taper to the impulse response of the source time
* function correction filter.
*/
void STFFourierDomainEngine::taperstf()
{
STFINV_debug(Mdebug&16, "STFFourierDomainEngine::taperstf()",
STFINV_value(Mapplystftaper));
if (this->Mapplystftaper)
{
// transform correction filter to time domain
Mfftenginestf.c2r();
Tfftengine::TAseries thestfseries=this->stfseries();
thestfseries *= Mfftenginestf.scale_series(this->dt());
// apply taper
/*
* Concept of application to periodic time series allowing for negative
* values of taper times:
*
* All samples between t2 and t3 remain unaltered. It is reasonable to
* process the taper starting at t3 passing t4 to t1 and then ending at
* t2. Behind this is the concept that time series are implicitely
* periodic in discrete Fourier transform.
*/
// samples in series
int nsamples=thestfseries.size(0);
// duration of time series
double T=this->dt()*nsamples;
// make sure taper definition is not longer than time series.
STFINV_assert((this->Mtt4-this->Mtt1)<T,
"Taper width is larger than length of time series");
// set taper time values under conpect of periodic series
double tt3=this->Mtt3;
double tt4=this->Mtt4;
double tt1=this->Mtt1+T;
double tt2=this->Mtt2+T;
// sample to start at
int l3=int(std::floor(tt3/this->dt()));
int l2=int(std::ceil(tt2/this->dt()));
STFINV_debug(Mdebug&16, "STFFourierDomainEngine::taperstf()",
STFINV_value(T) << "\n " <<
STFINV_value(this->dt()) << "\n " <<
STFINV_value(nsamples) << "\n " <<
STFINV_value(tt3) << "\n " <<
STFINV_value(tt4) << "\n " <<
STFINV_value(tt1) << "\n " <<
STFINV_value(tt2) << "\n " <<
STFINV_value(l3) << "\n " <<
STFINV_value(l2));
for (int l=l3; l<=l2; ++l)
{
// time of sample
double t=l