...
 
Commits (48)
this is <README>
============================================================================
AFF --- A container for numbers
-------------------------------
$Id$
============================================================================
For compilation instructions see README.1st in the root directory of the tar-ball or
http://gpitrsvn.gpi.uni-karlsruhe.de:8000/TFSoftware/wiki/docs/installation
This software is part of the project Seitosh. See README.1st in the root
directory of the collection or https://git.scc.kit.edu/Seitosh/Seitosh for
general installation instructions.
The AFF (Array of Friederich and Forbriger) is a lightweight class library. It
offers a simple and easy to use container for numbers as is necessary in
......
this is <COPYING>
============================================================================
libfourier
----------
$Id$
============================================================================
Copyright (c) 2002, 2013 by Thomas Forbriger
......
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id$
#
# Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
#
......@@ -76,6 +75,8 @@ $(LOCLIBDIR)/%: install-include %
#----------------------------------------------------------------------
DOXYTXT=$(wildcard doxygen*txt)
flist: Makefile $(wildcard *.f *.inc *.h *.cc *.c README *.cfg) COPYING \
$(DOXYTXT)
echo $^ | tr ' ' '\n' | sort > $@
......
this is <README>
============================================================================
FOURIER --- Fourier coefficients library
-------------------------------
$Id$
============================================================================
For compilation instructions see README.1st in the root directory of the tar-ball or
http://gpitrsvn.gpi.uni-karlsruhe.de:8000/TFSoftware/wiki/docs/installation
This software is part of the project Seitosh. See README.1st in the root
directory of the collection or https://git.scc.kit.edu/Seitosh/Seitosh for
general installation instructions.
The library provides modules to perform Fourier transformations and to
operate on Fourier transforms.An interface to FFTW is available.Signal
......
......@@ -3,8 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
*
* ----
* libfourier is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
......@@ -41,8 +39,6 @@
*/
#define TF_CFFTWTEST_C_VERSION \
"TF_CFFTWTEST_C V1.0 "
#define TF_CFFTWTEST_C_CVSID \
"$Id$"
#include<drfftw.h>
#include<stdio.h>
......@@ -105,7 +101,6 @@ void process(int n, int m)
int main()
{
printf("%s\n", TF_CFFTWTEST_C_VERSION);
printf("%s\n", TF_CFFTWTEST_C_CVSID);
/*
* the code provided in the tutorial for rfftw_one produces segmentation
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 16/05/2011
*
......@@ -37,8 +36,6 @@
*/
#define CXXFFTWARTEST_VERSION \
"CXXFFTWARTEST V1.2 test fftw3 array engine"
#define CXXFFTWARTEST_CVSID \
"$Id$"
// use the input/output facilities from the standard library
#include <iostream>
......@@ -140,7 +137,6 @@ int main(int iargc, char* argv[])
// define full help text
char help_text[]=
{
CXXFFTWARTEST_CVSID
"\n"
"input input file name\n"
"output output file name\n"
......@@ -218,7 +214,7 @@ int main(int iargc, char* argv[])
// report program version if in verbose mode
if (opt.verbose)
{ cout << CXXFFTWARTEST_VERSION << endl << CXXFFTWARTEST_CVSID << endl; }
{ cout << CXXFFTWARTEST_VERSION << endl; }
// extract commandline arguments
// here the rest of the command line is parsed; i.e. the names of input
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 12/09/2007
*
......@@ -35,8 +34,6 @@
*/
#define CXXFFTWTEST_VERSION \
"CXXFFTWTEST V1.1 a small test program for fftw when called from C++"
#define CXXFFTWTEST_CVSID \
"$Id$"
#include <iostream>
#include <tfxx/commandline.h>
......@@ -123,7 +120,6 @@ int main(int iargc, char* argv[])
// define full help text
char help_text[]=
{
CXXFFTWTEST_CVSID
"\n"
"-n n set number of samples to n\n"
"-m f set frequency to f\n"
......
......@@ -2,8 +2,6 @@
* \brief C++ library to operate on Fourier transforms (libfourierxx)
*
* ----------------------------------------------------------------------------
*
* $Id$
*
* Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
*
......@@ -34,8 +32,6 @@
/*! \mainpage
\author Thomas Forbriger
\since 2002
\date $Date$
\version $Id$
The library provides modules to perform Fourier transformations and to operate
on Fourier transforms.
......@@ -128,6 +124,14 @@ Further stuff will be collected in this library in the future.
* integral transform, where
* \f$\Delta\omega=2\pi/T\f$ and \f$T=N\Delta t\f$.
*
* The library provides functions which return the appropriate scaling factors
* such that application of a inverse transformation in sequence with a
* forward transformation is an identity-operation.
*
* \sa fourier::fft::DRFFTWAFFArrayEngine::scale_series,
* fourier::fft::DRFFTWAFFArrayEngine::scale_spectrum,
* fourier::fft::DRFFTWAFF::scale_series,
* fourier::fft::DRFFTWAFF::scale_spectrum
*/
// ----- END OF doxygen.txt -----
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 11/09/2007
*
......@@ -35,8 +34,6 @@
*/
#define FOURIER_ERROR_CC_VERSION \
"FOURIER_ERROR_CC V1.1"
#define FOURIER_ERROR_CC_CVSID \
"$Id$"
#include <iostream>
#include <fourier/error.h>
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 11/09/2007
*
......@@ -39,8 +38,6 @@
#define FOURIER_ERROR_H_VERSION \
"FOURIER_ERROR_H V1.0"
#define FOURIER_ERROR_H_CVSID \
"$Id$"
namespace fourier {
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 05/01/2003
*
......@@ -35,8 +34,6 @@
*/
#define TF_FCOMMAND_CC_VERSION \
"TF_FCOMMAND_CC V1.0 "
#define TF_FCOMMAND_CC_CVSID \
"$Id$"
#include <fourier/fcommand.h>
#include <fstream>
......@@ -79,11 +76,8 @@ namespace fourier {
{
os << endl;
os << TF_FCOMMAND_H_VERSION << endl;
os << TF_FCOMMAND_H_CVSID << endl;
os << TF_FILTERS_H_VERSION << endl;
os << TF_FILTERS_H_CVSID << endl;
os << TF_POLESNZEROES_H_VERSION << endl;
os << TF_POLESNZEROES_H_CVSID << endl;
}
}
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 05/01/2003
*
......@@ -39,8 +38,6 @@
#define TF_FCOMMAND_H_VERSION \
"TF_FCOMMAND_H V1.0 "
#define TF_FCOMMAND_H_CVSID \
"$Id$"
#include<iostream>
#include<fourier/filters.h>
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 11/07/2006
*
......@@ -42,8 +41,6 @@
*/
#define TF_FFTWAFF_CC_VERSION \
"TF_FFTWAFF_CC V1.4"
#define TF_FFTWAFF_CC_CVSID \
"$Id$"
#include <iostream>
#include <fourier/fftwaff.h>
......@@ -308,7 +305,15 @@ namespace fourier {
/*----------------------------------------------------------------------*/
//! Return appropriate scaling factor for sampling interval dt.
/*! \brief Return appropriate scaling factor for sampling interval dt.
*
* Factor to be applied when transforming to time domain.
*
* \param[in] dt sampling interval
* \return scalar factor to be applied to all samples
*
* \sa \ref sec_fftw3_integral_transform
*/
DRFFTWAFF::Tsample DRFFTWAFF::scale_series(const Tsample& dt) const
{
return(1./(Msize*dt));
......@@ -316,7 +321,15 @@ namespace fourier {
/*----------------------------------------------------------------------*/
//! Return appropriate scaling factor for sampling interval dt.
/*! \brief Return appropriate scaling factor for sampling interval dt.
*
* Factor to be applied when transforming to Fourier domain.
*
* \param[in] dt sampling interval
* \return scalar factor to be applied to all samples
*
* \sa \ref sec_fftw3_integral_transform
*/
DRFFTWAFF::Tsample DRFFTWAFF::scale_spectrum(const Tsample& dt) const
{
return(dt);
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 11/07/2006
*
......@@ -75,8 +74,6 @@
#define TF_FFTWAFF_H_VERSION \
"TF_FFTWAFF_H V1.3"
#define TF_FFTWAFF_H_CVSID \
"$Id$"
#include<complex>
#ifdef FFTWFALLBACK
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 13/05/2011
*
......@@ -37,8 +36,6 @@
*/
#define TF_FFTWAFFAR_CC_VERSION \
"TF_FFTWAFFAR_CC V1.1"
#define TF_FFTWAFFAR_CC_CVSID \
"$Id$"
#include <iostream>
#include <fourier/fftwaffar.h>
......@@ -253,7 +250,15 @@ namespace fourier {
/*----------------------------------------------------------------------*/
//! Return appropriate scaling factor for sampling interval dt.
/*! \brief Return appropriate scaling factor for sampling interval dt.
*
* Factor to be applied when transforming to time domain.
*
* \param[in] dt sampling interval
* \return scalar factor to be applied to all samples
*
* \sa \ref sec_fftw3_integral_transform
*/
DRFFTWAFFArrayEngine::Tsample
DRFFTWAFFArrayEngine::scale_series(const Tsample& dt) const
{
......@@ -262,7 +267,15 @@ namespace fourier {
/*----------------------------------------------------------------------*/
//! Return appropriate scaling factor for sampling interval dt.
/*! \brief Return appropriate scaling factor for sampling interval dt.
*
* Factor to be applied when transforming to Fourier domain.
*
* \param[in] dt sampling interval
* \return scalar factor to be applied to all samples
*
* \sa \ref sec_fftw3_integral_transform
*/
DRFFTWAFFArrayEngine::Tsample
DRFFTWAFFArrayEngine::scale_spectrum(const Tsample& dt) const
{
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 13/05/2011
*
......@@ -44,8 +43,6 @@
#define TF_FFTWAFFAR_H_VERSION \
"TF_FFTWAFFAR_H V1.2"
#define TF_FFTWAFFAR_H_CVSID \
"$Id$"
#include<complex>
#include<fftw3.h>
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 05/01/2003
*
......@@ -34,16 +33,15 @@
* ============================================================================
*/
#define TF_FILTER_CC_VERSION \
"TF_FILTER_CC V1.0 "
#define TF_FILTER_CC_CVSID \
"$Id$"
"TF_FILTER_CC V1.0"
#include <fourier/filters.h>
#include <fourier/error.h>
#include <cmath>
namespace fourier {
const double Filter::pi;
const double Filter::pi=M_PI;
const Filter::Tcvalue Filter::ime=Tcvalue(0.,1.);
/*----------------------------------------------------------------------*/
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 05/01/2003
*
......@@ -40,9 +39,7 @@
#ifndef TF_FILTERS_H_VERSION
#define TF_FILTERS_H_VERSION \
"TF_FILTERS_H V1.0 "
#define TF_FILTERS_H_CVSID \
"$Id$"
"TF_FILTERS_H V1.0"
#include<fourier/polesnzeroes.h>
......@@ -53,7 +50,7 @@ namespace fourier {
typedef PolesNZeroes Tbase;
typedef Tbase::Tcvalue Tcvalue;
static const double pi=3.141592653589793;
static const double pi;
static const Tcvalue ime;
Filter(): Tbase(), Mfrequency(false) { }
......
c this is <filters.inc>
c ----------------------------------------------------------------------------
c ($Id$)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 05/01/2003
*
......@@ -34,8 +33,6 @@
*/
#define TF_POLESNZEROES_CC_VERSION \
"TF_POLESNZEROES_CC V1.0 "
#define TF_POLESNZEROES_CC_CVSID \
"$Id$"
#include <fourier/polesnzeroes.h>
......
......@@ -3,7 +3,6 @@
*
* ----------------------------------------------------------------------------
*
* $Id$
* \author Thomas Forbriger
* \date 05/01/2003
*
......@@ -38,8 +37,6 @@
#define TF_POLESNZEROES_H_VERSION \
"TF_POLESNZEROES_H V1.0 "
#define TF_POLESNZEROES_H_CVSID \
"$Id$"
#include<cmath>
#include<complex>
......
c this is <polesnzeros.inc>
cS
c ----------------------------------------------------------------------------
c ($Id$)
c
c Copyright (c) 2002 by Thomas Forbriger (IMG Frankfurt)
c
......
this is <README>
============================================================================
SEIFE --- seismic waveform filters
----------------------------------
$Id$
============================================================================
For compilation instructions see README.1st in the root directory of the tar-ball or
http://gpitrsvn.gpi.uni-karlsruhe.de:8000/TFSoftware/wiki/docs/installation
For compilation instructions see README.1st in the root directory.
libseife is a collection of Fortran 77 functions and subroutines for time
series analysis and digital filters. It is a full-grown signal processing
......@@ -37,5 +34,15 @@ Dependencies
external 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 -----
......@@ -70,7 +70,7 @@ 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
approrpiate water level by application of the L-curve criterion (Groos, 2013,
Appendix G, page 148).
......@@ -87,6 +87,9 @@ tests/onlinehelp Provides access to these texts. Just issue
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
provided on the OpenTOAST web-page:
http://www.opentoast.de/Data_analysis_code_soutifu_and_libstfinv.php
......
......@@ -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 -----
......@@ -94,9 +94,9 @@
Programs using this library will require the following libraries in
addition:
- libfourierxx (available from TFSoftware)
- libfourierxx (available from Seitosh)
- libfftw3
- libaff (available from TFSoftware)
- libaff (available from Seitosh)
C programs will further be required to link against
......
......@@ -29,11 +29,12 @@
*
* REVISIONS and CHANGES
* - 06/05/2011 V1.0 Thomas Forbriger
* - 18/01/2016 V1.1 rename function to secomtospace
*
* ============================================================================
*/
#define STFINV_PARAMETERHANDLER_CC_VERSION \
"STFINV_PARAMETERHANDLER_CC V1.0"
"STFINV_PARAMETERHANDLER_CC V1.1"
#include <algorithm>
#include <stfinv/parameterhandler.h>
......@@ -50,7 +51,7 @@ std::string stfinv::tools::clipstring(std::string& s, const std::string& delim)
s.erase();
}
return(result);
} // std::string stfinv::tools::clipstring^
} // std::string stfinv::tools::clipstring()
/*----------------------------------------------------------------------*/
......@@ -71,11 +72,12 @@ stfinv::tools::Tparamap stfinv::tools::makeparamap(const std::string& p,
/*----------------------------------------------------------------------*/
std::string stfinv::tools::commatospace(std::string s)
std::string stfinv::tools::secomtospace(std::string s)
{
std::replace(s.begin(), s.end(), ',', ' ');
std::replace(s.begin(), s.end(), ';', ' ');
return(s);
} // std::string stfinv::tools::commatospace(const std::string& s)
} // std::string stfinv::tools::secomtospace(const std::string& s)
/*----------------------------------------------------------------------*/
......
......@@ -29,6 +29,7 @@
*
* REVISIONS and CHANGES
* - 06/05/2011 V1.0 Thomas Forbriger
* - 18/01/2016 V1.1 rename function to secomtospace
*
* ============================================================================
*/
......@@ -37,7 +38,7 @@
#ifndef STFINV_PARAMETERHANDLER_H_VERSION
#define STFINV_PARAMETERHANDLER_H_VERSION \
"STFINV_PARAMETERHANDLER_H V1.0"
"STFINV_PARAMETERHANDLER_H V1.1"
#include <string>
#include <map>
......@@ -45,12 +46,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 +66,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
......@@ -85,18 +86,19 @@ namespace stfinv {
/*----------------------------------------------------------------------*/
/*! replace comma by whitespace
* \ingroup tools
/*! replace commas and semicolons by whitespace
* \ingroup group_tools
*
* \param s input string
* \return input string with all commas replaced by whitespace
* \return input string with all commas replaced by whitespace and
* all semicolons replaced by whitespace
*/
std::string commatospace(std::string s);
std::string secomtospace(std::string s);
/*----------------------------------------------------------------------*/
/*! \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:
......
This diff is collapsed.
......@@ -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
*/
......
......@@ -19,6 +19,10 @@ remainder may consist of several control parameters being separated by colons
or may come along with a parameter value. The value is separated from the
parameter by an equal sign (=).
Where several values in an argument to a parameter must be separated (like in
the 'irtap' option of the Fourier domain procedures) white space ( ), commas
(,), and semicolons (;) are allowed as field delimiters, at your convenience.
Examples:
- To select Fourier domain least squares and shift the returned source
correction filter wavelet by 0.4s and switch on verbose mode, pass the
......
......@@ -4,17 +4,31 @@
# 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
These options define the number of samples N used for the FFT (Fast Fourier
Transform). This number N should be larger than the number of samples M in the
original input time series to avoid wrap-around. If fpow2 is set, N will be
the next power of 2 larger than M*f. Else if fdiv is set, N will be the next
integer multiple of d larger than M*f.
integer multiple of d larger than M*f. If fdiv is not set explicitely, a
default value for d (commonly 100) is used. If option fpad ist used, N will be
f times larger than without padding. Without explicitely setting fpad, a
default value for f is used (which commonly equals 1.5). When defining the
number of samples N, first padding is considered (fpad), then the either
selection of a power of two (pow2) or the divisor criterion (fdiv) is applied.
The latter is only applied, if pow2 ist not selected.
Input time series with M samples will be padded with (N-M) zeros to create the
time series which actually will be transformed to the Fourier domain. Upon
inverse FFT the additional (N-M) samples of the resulting time series will be
discarded before returning the M remaining samples to the caller. Note, that
this is a form of implicite taper. In particular the caller will not obtain
exactly the filter response, which was used for convolution internally.
The derived correction filter in some cases can contain acausal components.
This means that the impulse response is non-zero for negative time values.
......@@ -22,4 +36,24 @@ Since by definition, the impulse response is output for the time interval of
the input data, these acausal components can remain unnoticed. The option
tshift can be used to shift the impulse response as obtained by inverse FFT in
order to expose acausal components.
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
Time values are given in the same unit in which the sampling interval is given
in the input time series. I.e. if sampling interval is specified as a fraction
of seconds (which is standard) then all time values passed as parameters are
also given as fractions or multiples of seconds.
# ----- END OF stfinvfourier_description_usage.txt -----
......@@ -4,8 +4,9 @@
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)
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)
irtap=t1,t2,t3,t4 taper impulse response of correction filter
# ----- END OF stfinvfourier_summary_usage.txt -----
This diff is collapsed.
# VS VP RHO PRO F_FILT F_HIGH F_LOW WT J J EPRE EPSI
1 1 0 0 0.01 0 0 0 2 0 0 0 0.005
2 1 0 0 0.01 0 0 0 2 0 0 0 0.005
3 1 0 0 0.01 0 0 0 2 0 0 0 0.005
4 1 0 0 0.01 0 0 0 2 0 0 0 0.005
# VS VP RHO PRO F_FILT F_HIGH F_LOW WT J J EPRE EPSI GAMMA
1 1 0 0 0.01 0 0 0 2 0 0 0 0.005 20
2 1 0 0 0.01 0 0 0 2 0 0 0 0.005 10
3 1 0 0 0.01 0 0 0 2 0 0 0 0.005 5
4 1 0 0 0.01 0 0 0 2 0 0 0 0.005 0
This diff is collapsed.
......@@ -38,6 +38,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
extern float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
extern int EPRECOND;
extern float EPSILON_WE;
extern float GAMMA;
extern int GRAD_METHOD;
extern int WORKFLOW_STAGE;
......@@ -126,10 +127,12 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
if(EPRECOND==0 && workflow[WORKFLOW_STAGE][12]!=0){
if(MYID==0) printf(" WARNING: EPRECOND have to be set >0 in JSON (if so, ignore this message)");
}
EPRECOND=workflow[WORKFLOW_STAGE][12];
EPSILON_WE=workflow[WORKFLOW_STAGE][13];
GAMMA=workflow[WORKFLOW_STAGE][14];
if(*LBFGS_iter_start==*iter && GRAD_METHOD==2){
if(MYID==0)printf("\n L-BFGS will be used from iteration %d on.",*LBFGS_iter_start+1);
}
}
\ No newline at end of file
}
......@@ -29,7 +29,7 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
extern int TRKILL, NORMALIZE, F_LOW_PASS, TIMEWIN;
extern char TRKILL_FILE[STRING_SIZE];
extern int VELOCITY;
extern int WRITE_FILTERED_DATA;
int i,j;
float l2;
int h;
......@@ -202,6 +202,15 @@ double calc_misfit(float **sectiondata, float **section, int ntr, int ns, int LN
}
}
if(WRITE_FILTERED_DATA==2){
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectiondata[i][j]=intseis_sectiondata[i][j];
section[i][j]=intseis_section[i][j];
}
}
}
l2=L2;
/* printf("\n MYID = %i IN CALC_MISFIT: L2 = %10.12f \n",MYID,l2); */
......
......@@ -30,7 +30,7 @@ void exchange_par(void){
extern float XREC1, XREC2, YREC1, YREC2, FPML;
extern float REC_ARRAY_DEPTH, REC_ARRAY_DIST, MUN, EPSILON, EPSILON_u, EPSILON_rho;
extern int SEISMO, NDT, NGEOPH, SEIS_FORMAT, FREE_SURF, READMOD, READREC, SRCREC;
extern int BOUNDARY, REC_ARRAY, DRX, FW;
extern int BOUNDARY, REC_ARRAY, DRX, FW, STF_FULL;
extern int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
extern float TSNAP1, TSNAP2, TSNAPINC, REFREC[4];
extern char MFILE[STRING_SIZE], SIGNAL_FILE[STRING_SIZE],SIGNAL_FILE_SH[STRING_SIZE], LOG_FILE[STRING_SIZE];
......@@ -377,7 +377,7 @@ void exchange_par(void){
idum[116]=TRKILL_STF_OFFSET_INVERT;
idum[117]=JOINT_EQUAL_WEIGHTING;
idum[118]=STF_FULL;
} /** if (MYID == 0) **/
MPI_Barrier(MPI_COMM_WORLD);
......@@ -664,7 +664,7 @@ void exchange_par(void){
TRKILL_STF_OFFSET_INVERT=idum[116];
JOINT_EQUAL_WEIGHTING=idum[117];
STF_FULL=idum[118];
if ( MYID!=0 && L>0 ) {
FL=vector(1,L);
}
......
This diff is collapsed.
......@@ -14,7 +14,7 @@ float XREC1, XREC2, YREC1, YREC2;
float REC_ARRAY_DEPTH, REC_ARRAY_DIST;
float REFREC[4]={0.0, 0.0, 0.0, 0.0}, FPML;
int SEISMO, NDT, NGEOPH, NSRC=1, SEIS_FORMAT, FREE_SURF, READMOD, READREC, SRCREC, FW=0;
int NX, NY, NT, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, REC_ARRAY, RUN_MULTIPLE_SHOTS, NTRG;
int NX, NY, NT, SOURCE_SHAPE,SOURCE_SHAPE_SH, SOURCE_TYPE, SNAP, SNAP_FORMAT, REC_ARRAY, RUN_MULTIPLE_SHOTS, NTRG,STF_FULL;
int L, BOUNDARY, DC, DRX, NXG, NYG, IDX, IDY, FDORDER, MAXRELERROR;
char SNAP_FILE[STRING_SIZE], SOURCE_FILE[STRING_SIZE], SIGNAL_FILE[STRING_SIZE], SIGNAL_FILE_SH[STRING_SIZE];
char MFILE[STRING_SIZE], REC_FILE[STRING_SIZE];
......@@ -147,4 +147,4 @@ int JOINT_EQUAL_WEIGHTING;
float JOINT_INVERSION_PSV_SH_ALPHA_VS;
float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
\ No newline at end of file
int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
......@@ -21,7 +21,7 @@
#include "fd.h"
void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]){
void write_matrix_disk(float ** local_matrix,char path_name[STRING_SIZE]){
char joint[225];
FILE *FPjoint;
extern int POS[3],MYID;
......@@ -32,7 +32,7 @@ void write_matrix_disk(float ** gradient,char path_name[STRING_SIZE]){
for (i=1;i<=NX;i=i+IDX){
for (j=1;j<=NY;j=j+IDY){
fwrite(&gradient[j][i],sizeof(float),1,FPjoint);
fwrite(&local_matrix[j][i],sizeof(float),1,FPjoint);
}
}
......@@ -122,6 +122,70 @@ float matrix_product(float ** matrix1, float **matrix2) {
return global_sum;
}
float ** get_global_from_local_matrix(float ** local_matrix) {
extern int NXG, NYG;
extern int NX,NY;
extern int POS[3];
float ** global_matrix=NULL,** global_matrix_temp=NULL;
int i=0,j=0;
int ii=0, jj=0;
/* Allocate global matrix temp */
global_matrix_temp=matrix(1,NYG,1,NXG);
if(global_matrix_temp==NULL) {
declare_error("Allocation of global_matrix_temp in get_global_from_local_matrix failed!");
}
/* Allocate global matrix */
/* You have to deallocate this matrix on our own */
global_matrix=matrix(1,NYG,1,NXG);
if(global_matrix==NULL) {
declare_error("Allocation of global_matrix in get_global_from_local_matrix failed!");
}
/* Store local matrix in global matrix */
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
if ( (POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY)) ) {
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
global_matrix_temp[j][i]=local_matrix[jj][ii];
}
}
}
MPI_Allreduce(&global_matrix_temp[1][1],&global_matrix[1][1],NXG*NYG,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
free_matrix(global_matrix_temp,1,NYG,1,NXG);
return global_matrix;
}
void get_local_from_global_matrix(float ** global_matrix,float ** local_matrix) {
extern int NXG, NYG;
extern int NX,NY;
extern int POS[3];
int i=0,j=0;
int ii=0, jj=0;
/* Store local matrix in global matrix */
for (i=1;i<=NXG;i++){
for (j=1;j<=NYG;j++){
if ( (POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY)) ) {
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
local_matrix[jj][ii]=global_matrix[j][i];
}
}
}
}
......@@ -25,6 +25,7 @@
float *rd_sour(int *nts,FILE* fp_source){
extern int VERBOSE;
/* local variables */
float *psource;
int i, c;
......@@ -35,7 +36,7 @@ float *rd_sour(int *nts,FILE* fp_source){
while ((c=fgetc(fp_source)) != EOF)
if (c=='\n') ++(*nts);
rewind(fp_source);
printf(" Number of samples (nts) in source file: %i\n",*nts);
if (VERBOSE==1) printf(" Number of samples (nts) in source file: %i\n",*nts);
psource=vector(1,*nts);
for (i=1;i<=*nts;i++) fscanf(fp_source,"%e\n",&psource[i]);
......
......@@ -128,6 +128,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern float WOLFE_C1_SL;
extern float WOLFE_C2_SL;
extern int STF_FULL;
/* definition of local variables */
int number_readobjects=0,fserr=0;
......@@ -462,7 +463,9 @@ void read_par_json(FILE *fp, char *fileinp){
declare_error("Variable FORWARD_ONLY could not be retrieved from the json input file!");
else {
if (FORWARD_ONLY==0) { /* FWI is calculated */
/* Overwrite IDX/IDY option from forward modeling (used for snapshots), interpolation for FWI not yet implemented*/
IDX=1;
IDY=1;
/* General inversion parameters */
if (get_int_from_objectlist("ITERMAX",number_readobjects,&ITERMAX,varname_list, value_list))
declare_error("Variable ITERMAX could not be retrieved from the json input file!");
......@@ -763,6 +766,9 @@ void read_par_json(FILE *fp, char *fileinp){
TRKILL_STF=0;
fprintf(fp,"Variable TRKILL_STF is set to default value %d.\n",TRKILL_STF);}
else {
if (get_int_from_objectlist("STF_FULL",number_readobjects,&STF_FULL,varname_list, value_list)){
STF_FULL=0;
fprintf(fp,"Variable STF_FULL is set to default value %d.\n",STF_FULL);}
if (TRKILL_STF==1) {
if (get_int_from_objectlist("TRKILL_STF_OFFSET",number_readobjects,&TRKILL_STF_OFFSET,varname_list, value_list)){
TRKILL_STF_OFFSET=0;
......
......@@ -104,7 +104,7 @@ void readmod_elastic(float ** rho, float ** pi, float ** u){
}
if(feof(fp_vs) && feof(fp_rho)){
if(feof(fp_vs) || feof(fp_rho)){
declare_error("Model file VS or RHO is to small. Check dimensions NX*NY of file.");
}
......@@ -141,7 +141,7 @@ void readmod_elastic(float ** rho, float ** pi, float ** u){
fread(&vs, sizeof(float), 1, fp_vs);
fread(&rho, sizeof(float), 1, fp_rho);
if(!feof(fp_vs) && !feof(fp_rho)){
if(!feof(fp_vs) || !feof(fp_rho)){
declare_error("Model file VS or RHO is to big. Check dimensions NX*NY of file.");
}
fclose(fp_vs);
......
This diff is collapsed.
This diff is collapsed.
......@@ -28,7 +28,7 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
/* declaration of global variables */
extern float DT, DH;
extern int SEIS_FORMAT, MYID, NT, SOURCE_SHAPE, TIME_FILT, TIMEWIN, TAPER_STF, ORDER;
extern int SEIS_FORMAT, MYID, NT, SOURCE_SHAPE, TIME_FILT, TIMEWIN, TAPER_STF, ORDER, STF_FULL;
extern char PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
extern int TRKILL_STF, NORMALIZE, USE_WORKFLOW, WORKFLOW_STAGE;
extern char TRKILL_FILE_STF[STRING_SIZE];
......@@ -160,10 +160,10 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
}
/* trace killing ends here */
if(TIMEWIN==1){
if((TIMEWIN==1)&&(STF_FULL==0)){
time_window_glob(sectionvy, iter, ntr_glob, ns, ishot);
time_window_glob(sectionvy_obs, iter, ntr_glob, ns, ishot);
}
}
/* NORMALIZE TRACES */
if(NORMALIZE==1){
......@@ -358,4 +358,4 @@ void stf(FILE *fp, float **sectionvy, float ** sectionvy_obs, float ** sectionvy
free_imatrix(kill_tmp,1,ntr_glob,1,nshots);
free_ivector(kill_vector,1,ntr_glob);
}
}
\ No newline at end of file
}
......@@ -139,10 +139,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc,
time = (float)(j * DT);
dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
taper = exp(-GAMMA*dumpa*dumpa);
taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]);
taper1 = exp(-GAMMA*dumpb*dumpb);
taper1 = exp(-GAMMA*dumpb);
if(time>=picked_times_m[1][i]+picked_times_m[2][i]){
sectiondata[i][j] = sectiondata[i][j] * taper;}
......@@ -160,10 +160,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc,
time = (float)(j * DT);
dumpa = (time-picked_times_m[1][i]-picked_times_m[2][i]);
taper = exp(-GAMMA*dumpa*dumpa);
taper = exp(-GAMMA*dumpa);
dumpb = (time-picked_times_m[1][i]+picked_times_m[3][i]);
taper1 = exp(-GAMMA*dumpb*dumpb);
taper1 = exp(-GAMMA*dumpb);
dummysection[i][j] = sectiondata[i][j];
......@@ -174,10 +174,10 @@ void time_window(float **sectiondata, int iter, int ntr_glob, int **recpos_loc,
dummysection[i][j] = sectiondata[i][j] * taper1;}
dumpc = (time-picked_times_m[4][i]-picked_times_m[5][i]);
taper2 = exp(-GAMMA*dumpc*dumpc);
taper2 = exp(-GAMMA*dumpc);
dumpd = (time-picked_times_m[4][i]+picked_times_m[6][i]);
taper3 = exp(-GAMMA*dumpd*dumpd);
taper3 = exp(-GAMMA*dumpd);
<