Commit 245aa343 authored by thomas.forbriger's avatar thomas.forbriger

libstfinv [MERGE][FEATURE]: make taper available on master

Merge branch 'libstfinv_fdlsq' into master

A new taper option for procedures operating in the Fourier domain is
implemenetd in libstfinv. This commit makes the new feature available on
master.
parents fd9e7c84 3fb21daf
......@@ -40,7 +40,7 @@
"STFINV_DEBUG_H V1.0"
/*! \brief produce debug output
* \ingroup misc_h
* \ingroup group_debug
*
* \param C output will be generated if C == true
* \param N name of function
......@@ -54,8 +54,15 @@
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
* \ingroup misc_h
* \ingroup group_debug
*
* \param P parameter to dump
*/
......
......@@ -104,7 +104,7 @@ processing are described in the \ref page_help.
/*! \brief Engines implemented in libstfinv
\defgroup engines Engines
\defgroup group_engines Engines
\todo
A detailed description for implementers is still missing
......@@ -116,7 +116,7 @@ processing are described in the \ref page_help.
/*! \brief Tools and utilities used by the libstfinv engines
\defgroup tools Internal tools and utilities
\defgroup group_tools Internal tools and utilities
\todo
A detailed description for implementers is still missing
......@@ -136,4 +136,18 @@ processing are described in the \ref page_help.
\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 -----
......@@ -45,12 +45,12 @@
namespace stfinv {
/*! \brief Namespace for internal tools
* \ingroup tools
* \ingroup group_tools
*/
namespace tools {
/*! strip substring
* \ingroup tools
* \ingroup group_tools
*
* Strips off first substring up to given delimiter.
* The string is passed as a reference and will be modified (i.e. the
......@@ -65,14 +65,14 @@ namespace stfinv {
/*----------------------------------------------------------------------*/
/*! \brief A map to store parameters.
* \ingroup tools
* \ingroup group_tools
*/
typedef std::map<std::string,std::string> Tparamap;
/*----------------------------------------------------------------------*/
/*! \brief Create a parameter map from a parameter string
* \ingroup tools
* \ingroup group_tools
*
* \param p parameter string
* \param delimiter delimiter which separates two parameters
......@@ -86,7 +86,7 @@ namespace stfinv {
/*----------------------------------------------------------------------*/
/*! replace comma by whitespace
* \ingroup tools
* \ingroup group_tools
*
* \param s input string
* \return input string with all commas replaced by whitespace
......@@ -96,7 +96,7 @@ namespace stfinv {
/*----------------------------------------------------------------------*/
/*! \brief remove leading and trailing whitespace
* \ingroup tools
* \ingroup group_tools
*
* \param s any string
* \return value a input string with any leading and trailing whitespace
......
......@@ -47,7 +47,7 @@
namespace stfinv {
/*! \brief Fourier domain least squares engine
* \ingroup engines
* \ingroup group_engines
*
* \par Concept behind this engine
* If
......
......@@ -30,7 +30,7 @@
namespace stfinv {
/*! \brief Engine to find a finite, causal source time-history in time domain
* \ingroup engines
* \ingroup group_engines
*
* \par Concept behin this engine
*
......
......@@ -44,7 +44,7 @@
namespace stfinv {
/*! \brief Engine to provide a fixed wavelet
* \ingroup engines
* \ingroup group_engines
*/
class STFEngineFixedWavelet: public stfinv::STFFourierDomainEngine {
public:
......
......@@ -32,13 +32,20 @@
* - 30/09/2011 V1.1 implemented handling of additional time series pairs
* - 04/10/2011 V1.2 correction in debug message
* - 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 \
"STFINV_STFINVFOURIER_CC V1.3"
"STFINV_STFINVFOURIER_CC V1.5"
#include <sstream>
#include <cmath>
#include <stfinv/stfinvfourier.h>
#include <stfinv/stfinvfourier_summary_usage.h>
#include <stfinv/stfinvfourier_description_usage.h>
......@@ -165,13 +172,21 @@ namespace stfinv {
* next integer multiple of the divisor larger than fpad*nsamples.
*
* \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
* by recorded data and synthetic data, because both have to be transformed
* to Fourier domain at once.
* Transformation to Fourier domain takes place in
* STFFourierDomainEngine::fftinput.
* -# One engine (\c STFFourierDomainEngine::Mfftengineoutput) being shared
* by the stf and the convolved synthetics, because both have to be
* transformed to the time domain at once.
* Transformation to time domain takes place in
* STFFourierDomainEngine::fftoutput.
*/
void STFFourierDomainEngine::initialize()
{
......@@ -226,6 +241,8 @@ namespace stfinv {
Mfftengineoutput=fourier::fft::DRFFTWAFFArrayEngine(
1+this->nreceivers()+this->npairs(),
nsamples);
Mfftenginestf=fourier::fft::DRFFTWAFFArrayEngine(this->stfseries(),
this->stfspec());
// set time shift for STF if requested
// -----------------------------------
......@@ -234,28 +251,70 @@ namespace stfinv {
is >> Mtshift;
}
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()",
"this->nsamples()=" << this->nsamples() << " " <<
"padfactor=" << padfactor << " " <<
"divisor=" << divisor << " " <<
"nsamples=" << nsamples << " " <<
"Mapplyshift=" << Mapplyshift << " " <<
"Mtshift=" << Mtshift << " ");
STFINV_value(this->nsamples()) << "\n " <<
STFINV_value(padfactor) << "\n " <<
STFINV_value(divisor) << "\n " <<
STFINV_value(nsamples) << "\n " <<
STFINV_value(Mapplyshift) << "\n " <<
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()
/*----------------------------------------------------------------------*/
/*!
* 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()
{
this->getinput();
Mfftengineinput.r2c();
Mfftengineoutput.series()=0.;
} // 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()
{
if (this->Mapplystftaper) { this->taperstf(); }
this->convolve();
if (Mapplyshift) { this->stfshift(); }
Mfftengineoutput.c2r();
......@@ -296,6 +355,13 @@ namespace stfinv {
/*----------------------------------------------------------------------*/
STFFourierDomainEngine::TAseries STFFourierDomainEngine::stfseries() const
{
return(Mfftengineoutput.series(this->nreceivers()));
} // STFFourierDomainEngine::Tspectrum STFFourierDomainEngine::stfseries() const
/*----------------------------------------------------------------------*/
STFFourierDomainEngine::TAspectrum
STFFourierDomainEngine::recordingcoeff(const unsigned int& i) const
{
......@@ -353,6 +419,8 @@ namespace stfinv {
void STFFourierDomainEngine::getinput()
{
// 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();
sarray=0.;
......@@ -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()
{
// scaling factor for source correction filter
......@@ -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()
{
TAspectrum synthetic=this->syntheticspec();
......@@ -483,6 +568,96 @@ namespace stfinv {
} // for (unsigned int j=0; j<this->nfreq(); ++j)
} // 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*this->dt();
// index to series
int rl = l%nsamples;
if (rl < 0) { rl += nsamples; }
STFINV_assert(((rl >= thestfseries.f(0)) &&
(rl <= thestfseries.l(0))),
"Index out of range. Internal programming error! "
"Report as bug!")
// taper factor
double factor=1.;
if ( (t>=tt3) && (t<=tt4) && (tt4>tt3) )
{
factor=0.5+0.5*cos(M_PI*(t-tt3)/(tt4-tt3));
}
else if ( (t>tt4) && (t<tt1) )
{
factor=0.;
}
else if ( (t>=tt1) && (t<=tt2) && (tt2>tt1) )
{
factor=0.5-0.5*cos(M_PI*(t-tt1)/(tt2-tt1));
}
// apply to series
thestfseries(rl) = thestfseries(rl)*factor;
STFINV_debug(Mdebug&16, "STFFourierDomainEngine::taperstf()",
STFINV_value(t) << " " <<
STFINV_value(l) << " " <<
STFINV_value(rl) << " " <<
STFINV_value(factor));
}
// transform correction filter to Fourier domain
Mfftenginestf.r2c();
this->stfspec() *= Mfftenginestf.scale_spectrum(this->dt());
} // if (this->Mapplystftaper)
} // void STFFourierDomainEngine::taperstf()
} // namespace stfinv
/* ----- END OF stfinvfourier.cc ----- */
......@@ -30,6 +30,9 @@
* - 08/05/2011 V1.0 Thomas Forbriger
* - 30/09/2011 V1.1 implemented handling of additional time series pairs
* - 14/10/2015 V1.2 new end-user usage functions
* - 28/06/2016 V1.3 provide time domain tapering of filter response
* - 22/07/2016 V1.4 provide separate FFT processor addressing just the
* source time function correction filter
*
* ============================================================================
*/
......@@ -38,7 +41,7 @@
#ifndef STFINV_STFINVFOURIER_H_VERSION
#define STFINV_STFINVFOURIER_H_VERSION \
"STFINV_STFINVFOURIER_H V1.2"
"STFINV_STFINVFOURIER_H V1.4"
#include <stfinv/stfinvbase.h>
#include <aff/array.h>
......@@ -51,12 +54,19 @@ namespace stfinv {
* This is just a base class.
* The constructor is protected and should only be called from a derived
* class.
* The intention of this class is to provide all processing steps common to
* all engines operating in the Fourier domain in one single base class.
* The individual engines of different Fourier domain approaches then need
* not reimplement these steps.
* They essentially need only provide a specific exec-function (e.g.
* STFEngineFDLeastSquares::exec).
*
* This class maintains a workspace for Fourier transforms.
* It provides the FFT from input signals to the workspace through a member
* functions as well as the convolution of the synthetic data with a given
* source wavelet Fourier transform and a subsequent FFT to time domain for
* the convolved synthetics as well as the source correction filter separately.
* the convolved synthetics as well as the source correction filter
* separately.
*
* \par What STFFourierDomainEngine does for you
* All derived classes call STFFourierDomainEngine::fftinput prior to
......@@ -69,13 +79,27 @@ namespace stfinv {
* \par
* When processing has finished, the derived classes should call
* STFFourierDomainEngine::fftoutput.
* This function convolves the synthetic data with the source correction
* This function first applies a time domain taper to the correction filter
* impulse response if requested (STFFourierDomainEngine::taperstf).
* Then it convolves the synthetic data with the source correction
* filter (STFFourierDomainEngine::convolve).
* Then it applies a time shift to the source correction filter if requested
* (STFFourierDomainEngine::stfshift).
* If requested it applies a time shift to the source correction filter
* as a next step (STFFourierDomainEngine::stfshift).
* The convolved synthetics as well as the source correction filter then are
* transformed to time domain and written to the users workspace
* ((STFFourierDomainEngine::putoutput).
* (STFFourierDomainEngine::putoutput).
*
* \par
* This should take place in the exec-function (e.g.
* STFEngineFDLeastSquares::exec) of the derived class.
* I.e. the first statement in the exec function is a call to function
* STFFourierDomainEngine::fftinput of the base class and the very last
* statement is a call to function STFFourierDomainEngine::fftoutput of the
* base class.
* This also guarantees that STFFourierDomainEngine::fftoutput is only
* called once per derived correction filter response.
* This is necessary, since otherwise the taper function and the time shift
* would be applied twice to the impulse response.
*
* \par Layout of Fourier transform arrays
* The workspace for the Fourier transform engine is initialized by
......@@ -119,6 +143,12 @@ namespace stfinv {
//! \brief return name of engine
virtual const char* name() const;
protected:
/*! \name Access and control functions to be used by derived classes.
*
* These functions are part of the interface implemented in
* STFFourierDomainEngine.
*/
//@{
/*! \brief copy input signals to workspace and
* transform input workspace to Fourier domain
*/
......@@ -150,7 +180,14 @@ namespace stfinv {
double frequency(const unsigned int& i) const;
//! \brief return number of frequencies in use
unsigned int nfreq() const;
//@}
private:
/*! \name Internal processing control functions of.
*
* These functions are part of the interface implemented in
* STFFourierDomainEngine.
*/
//@{
//! \brief initialize work space
void initialize();
/*! \brief copy input time series for recorded data and synthetics
......@@ -167,6 +204,12 @@ namespace stfinv {
/*! \brief apply time shift to stf prior to FFT to time domain
*/
void stfshift();
/*! \brief apply a time domain taper to the correction filter response.
*/
void taperstf();
//! \brief return reference to time series container of stf
TAseries stfseries() const;
//@}
// member data
// -----------
......@@ -223,12 +266,43 @@ namespace stfinv {
* M is returned by function npairs().
*/
Tfftengine Mfftengineoutput;
/*! \brief FFT processor for source time function correction filter
*
* This uses a reference to the source time function correction filter
* data in Mfftengineoutput. It is used in cases, where this data has to
* be transformed alone (like in STFFourierDomainEngine::taperstf).
*
* \note
* This processor does not maintain a separate data space.
* It rather operates on a reference to data space also maintained by
* Mfftengineoutput.
*/
Tfftengine Mfftenginestf;
/*! \brief time shift to be applied to STF in order to expose
* acausal parts
*/
double Mtshift;
//! \brief true if shift must be applied
bool Mapplyshift;
/*! \brief true if time domain taper should be applied to filter
* response.
*/
bool Mapplystftaper;
/*! \brief time values defining taper.
*
* All samples at times
* - t<Mtt1 will be set to zero.
* - Mtt1<=t<=Mtt2 will be scaled by
* 0.5-0.5*cos(pi*(t-Mtt1)/(Mtt2-Mtt1)).
* - Mtt2<t<Mtt3 will remain unaltered.
* - Mtt3<=t<=Mtt4 will be scaled by
* 0.5+0.5*cos(pi*(t-Mtt3)/(Mtt4-Mtt3)).
* - t>Mtt4 will be set to zero.
*
* @{
*/
double Mtt1, Mtt2, Mtt3, Mtt4;
//!@}
}; // class STFFourierDomainEngine
}
......
......@@ -45,7 +45,7 @@
namespace stfinv {
/*! \brief Engine to apply a scalar factor
* \ingroup engines
* \ingroup group_engines
*
* \par Concept behin this engine
* This engine convolves the synthetic data with a discrete delta pulse so
......
......@@ -45,7 +45,7 @@
namespace stfinv {
/*! \brief Normalization engine
* \ingroup engines
* \ingroup group_engines
*
* \par Motivation
* On the down-side of Fourier domain least squares as is implemented
......
/*! \file tools.cc
* \brief tools and utilities (implementation)
*
* \ingroup tools
* \ingroup group_tools
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
......
/*! \file tools.h
* \brief tools and utilities (prototypes)
*
* \ingroup tools
* \ingroup group_tools
* ----------------------------------------------------------------------------
*
* \author Thomas Forbriger
......@@ -48,7 +48,7 @@ namespace stfinv {
namespace tools {
/*! \brief function to compare doubles
* \ingroup tools
* \ingroup group_tools
* \param a a double value
* \param b a double value
* \param eps relative residual allowed for \c a and \c b
......@@ -60,7 +60,7 @@ namespace stfinv {
/* ---------------------------------------------------------------------- */
/*! \brief report engine identifier
* \ingroup tools
* \ingroup group_tools
* \param C class to report ID and oneline description
* \param os output stream to send output to
*/
......@@ -78,7 +78,7 @@ namespace stfinv {
/* ---------------------------------------------------------------------- */
/*! \brief report engine identifier with heading
* \ingroup tools
* \ingroup group_tools
* \param C class to report ID and oneline description
* \param os output stream to send output to
*/
......
......@@ -4,11 +4,12 @@
# Procedures in the Fourier domain
# --------------------------------
Options and parameters in common for procedures in the Fourier domain:
fpow2 use power of two for number of coefficients
fdiv=d use integer multiple of d for number of coefficients
fpad=f padding factor
tshift=d delay source correction filter wavelet by d (in seconds)
in order to expose acausal components
fpow2 use power of two for number of coefficients
fdiv=d use integer multiple of d for number of coefficients
fpad=f padding factor
tshift=d delay source correction filter wavelet by d (in seconds)
in order to expose acausal components
irtap=t1,t2,t3,t4 taper impulse response of correction filter